text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Testing

function get_model()
    X = Dict{Int, Vector{Float64}}(
        1 => [2.0, 3.0, 4.0],
        2 => [6.0, 3.0, 2.0],
        3 => [2.0, 5.0, 1.0],
        4 => [4.0, 3.0, 6.0])
    u = Dict{Int, Vector{Float64}}(
        1 => [0.0, 0.0, 0.0],
        2 => [0.0, 0.0, 0.0],
        3 => [0.0, 0.0, 0.0],
        4 => [0.25, 0.25, 0.25])
    e1 = Element(Tet4, [1, 2, 3, 4])
    e2 = Element(Tri3, [1, 2, 3])
    update!([e1, e2], "geometry", X)
    update!([e1, e2], "displacement", 0.0 => u)
    update!(e1, "youngs modulus" => 96.0)
    update!(e1, "poissons ratio" => 1.0/3.0)
    update!(e1, "density" => 420.0)
    update!(e2, "displacement 1" => 0.0)
    update!(e2, "displacement 2" => 0.0)
    update!(e2, "displacement 3" => 0.0)
    p1 = Problem(Elasticity, 3)
    p1.properties.finite_strain = false
    p1.properties.geometric_stiffness = false
    p2 = Problem(Dirichlet, p1)
    push!(p1, e1)
    push!(p2, e2)
    solver = Solver(Modal)
    solver.properties.which = :LM
    push!(solver, p1, p2)
    return solver
end

@testset "test eigenvalues for single tet4 element" begin
    solver = get_model()
    solver(; debug=true)
    @test isapprox(solver.properties.eigvals, [4/3, 1/3])
end

@testset "test eigenvalues for single tet4 element, with geometric stiffness" begin
    solver = get_model()
    problem = first(solver.problems)
#   problem.properties.finite_strain = true
    problem.properties.geometric_stiffness = true
    solver.properties.geometric_stiffness = true
    solver(; debug=true)
    @test isapprox(solver.properties.eigvals, [5/3, 2/3])
end

@testset "test poisson problem modal analysis without tie" begin
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [1.0, 3.0],
        4 => [0.0, 3.0],
        5 => [0.0, 3.0],
        6 => [1.0, 3.0],
        7 => [1.0, 9.0],
        8 => [0.0, 9.0])
    el1 = Element(Quad4, [1, 2, 3, 4])
    el2 = Element(Quad4, [4, 3, 7, 8])
    el3 = Element(Seg2, [1, 2])
    el4 = Element(Seg2, [7, 8])
    update!([el1, el2, el3, el4], "geometry", X)
    update!([el1, el2], "density", 6.0)
    update!([el1, el2], "temperature thermal conductivity", 36.0)
    update!([el3, el4], "temperature 1", 0.0)
    p1 = Problem(Heat, "combined body", 1)
    p1.properties.formulation = "2D"
    p2 = Problem(Dirichlet, "fixed ends", 1, "temperature")
    push!(p1, el1, el2)
    push!(p2, el3, el4)
    solver = Solver(Modal)
    push!(solver, p1, p2)
    solver()
    @test isapprox(solver.properties.eigvals[1], 1.0)
end

@testset "test poisson modal problem with mesh tie" begin
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [1.0, 3.0],
        4 => [0.0, 3.0],
        5 => [0.0, 3.0],
        6 => [1.0, 3.0],
        7 => [1.0, 9.0],
        8 => [0.0, 9.0])
    el1 = Element(Quad4, [1, 2, 3, 4])
    el2 = Element(Quad4, [5, 6, 7, 8])
    el3 = Element(Seg2, [1, 2])
    el4 = Element(Seg2, [7, 8])
    el5 = Element(Seg2, [3, 4])
    el6 = Element(Seg2, [5, 6])
    update!([el1, el2, el3, el4, el5, el6], "geometry", X)
    update!([el1, el2], "density", 6.0)
    update!([el1, el2], "temperature thermal conductivity", 36.0)
    update!([el3, el4], "temperature 1", 0.0)
    update!(el5, "master elements", [el6])
    p1 = Problem(Heat, "body 1", 1)
    p2 = Problem(Heat, "body 2", 1)
    p1.properties.formulation = "2D"
    p2.properties.formulation = "2D"
    p3 = Problem(Dirichlet, "fixed ends", 1, "temperature")
    p4 = Problem(Mortar, "interface between bodies", 1, "temperature")
    p4.properties.dimension = 1
    push!(p1, el1)
    push!(p2, el2)
    push!(p3, el3, el4)
    push!(p4, el5, el6)
    solver = Solver(Modal)
    push!(solver, p1, p2, p3, p4)
    solver()
    @test isapprox(solver.properties.eigvals[1], 1.0)
end


text:
module TestDataTable
    using Base.Test
    using DataTables, Compat
    import Compat.String

    #
    # Equality
    #

    @test isequal(DataTable(a=[1, 2, 3], b=[4, 5, 6]), DataTable(a=[1, 2, 3], b=[4, 5, 6]))
    @test !isequal(DataTable(a=[1, 2], b=[4, 5]), DataTable(a=[1, 2, 3], b=[4, 5, 6]))
    @test !isequal(DataTable(a=[1, 2, 3], b=[4, 5, 6]), DataTable(a=[1, 2, 3]))
    @test !isequal(DataTable(a=[1, 2, 3], b=[4, 5, 6]), DataTable(a=[1, 2, 3], c=[4, 5, 6]))
    @test !isequal(DataTable(a=[1, 2, 3], b=[4, 5, 6]), DataTable(b=[4, 5, 6], a=[1, 2, 3]))
    @test !isequal(DataTable(a=[1, 2, 2], b=[4, 5, 6]), DataTable(a=[1, 2, 3], b=[4, 5, 6]))
    @test isequal(DataTable(a=Nullable{Int}[1, 2, Nullable()], b=[4, 5, 6]),
                  DataTable(a=Nullable{Int}[1, 2, Nullable()], b=[4, 5, 6]))

    # FIXME: equality operators won't work until JuliaStats/NullableArrays#84 is merged
    #@test get(DataTable(a=[1, 2, 3], b=[4, 5, 6]) == DataTable(a=[1, 2, 3], b=[4, 5, 6]))
    #@test get(DataTable(a=[1, 2], b=[4, 5]) != DataTable(a=[1, 2, 3], b=[4, 5, 6]))
    #@test get(DataTable(a=[1, 2, 3], b=[4, 5, 6]) != DataTable(a=[1, 2, 3]))
    #@test get(DataTable(a=[1, 2, 3], b=[4, 5, 6]) != DataTable(a=[1, 2, 3], c=[4, 5, 6]))
    #@test get(DataTable(a=[1, 2, 3], b=[4, 5, 6]) != DataTable(b=[4, 5, 6], a=[1, 2, 3]))
    #@test get(DataTable(a=[1, 2, 2], b=[4, 5, 6]) != DataTable(a=[1, 2, 3], b=[4, 5, 6]))
    #@test get(DataTable(a=Nullable{Int}[1, 3, Nullable()], b=[4, 5, 6]) !=
    #          DataTable(a=Nullable{Int}[1, 2, Nullable()], b=[4, 5, 6]))
    #@test isnull(DataTable(a=Nullable{Int}[1, 2, Nullable()], b=[4, 5, 6]) ==
    #             DataTable(a=Nullable{Int}[1, 2, Nullable()], b=[4, 5, 6]))
    #@test isnull(DataTable(a=Nullable{Int}[1, 2, Nullable()], b=[4, 5, 6]) ==
    #             DataTable(a=Nullable{Int}[1, 2, 3], b=[4, 5, 6]))

    #
    # Copying
    #

    dt = DataTable(a = NullableArray([2, 3]),
                   b = NullableArray([DataTable(c = 1), DataTable(d = 2)]))
    dtc = copy(dt)
    dtdc = deepcopy(dt)

    dt[1, :a] = 4
    get(dt[1, :b])[:e] = 5
    names!(dt, [:f, :g])

    @test names(dtc) == [:a, :b]
    @test names(dtdc) == [:a, :b]

    @test get(dtc[1, :a]) === 4
    @test get(dtdc[1, :a]) === 2

    @test names(get(dtc[1, :b])) == [:c, :e]
    @test names(get(dtdc[1, :b])) == [:c]

    #

    x = DataTable(a = [1, 2, 3], b = [4, 5, 6])
    v = DataTable(a = [5, 6, 7], b = [8, 9, 10])

    z = vcat(v, x)

    z2 = z[:, [1, 1, 2]]
    @test names(z2) == [:a, :a_1, :b]

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

    # Insert single value
    x[:d] = Nullable(3)
    @test isequal(x[:d], NullableArray([3, 3, 3]))

    x0[:d] = 3
    @test x0[:d] == Int[]

    # similar / nulls
    dt = DataTable(a = NullableArray([1]),
                   b = NullableArray(["b"]),
                   c = NullableCategoricalArray([3.3]))
    nulldt = DataTable(a = NullableArray{Int}(2),
                       b = NullableArray{String}(2),
                       c = NullableCategoricalArray{Float64}(2))
    @test isequal(nulldt, similar(dt, 2))

    # Associative methods

    dt = DataTable(a=[1, 2], b=[3., 4.])
    @test haskey(dt, :a)
    @test !haskey(dt, :c)
    @test get(dt, :a, -1) === dt.columns[1]
    @test get(dt, :c, -1) == -1
    @test !isempty(dt)

    @test empty!(dt) === dt
    @test isempty(dt.columns)
    @test isempty(dt)

    dt = DataTable(a=NullableArray([1, 2]), b=NullableArray([3., 4.]))
    @test_throws BoundsError insert!(dt, 5, ["a", "b"], :newcol)
    @test_throws ErrorException insert!(dt, 1, ["a"], :newcol)
    @test isequal(insert!(dt, 1, ["a", "b"], :newcol), dt)
    @test names(dt) == [:newcol, :a, :b]
    @test isequal(dt[:a], NullableArray([1, 2]))
    @test isequal(dt[:b], NullableArray([3., 4.]))
    @test isequal(dt[:newcol], ["a", "b"])

    dt = DataTable(a=[1, 2], b=[3., 4.])
    dt2 = DataTable(b=["a", "b"], c=[:c, :d])
    @test isequal(merge!(dt, dt2), dt)
    @test isequal(dt, DataTable(a=[1, 2], b=["a", "b"], c=[:c, :d]))

    #test_group("Empty DataTable constructors")
    dt = DataTable(Nullable{Int}, 10, 3)
    @test size(dt, 1) == 10
    @test size(dt, 2) == 3
    @test typeof(dt[:, 1]) == NullableVector{Int}
    @test typeof(dt[:, 2]) == NullableVector{Int}
    @test typeof(dt[:, 3]) == NullableVector{Int}
    @test all(isnull, dt[:, 1])
    @test all(isnull, dt[:, 2])
    @test all(isnull, dt[:, 3])

    dt = DataTable([Nullable{Int}, Nullable{Float64}, Nullable{String}], 100)
    @test size(dt, 1) == 100
    @test size(dt, 2) == 3
    @test typeof(dt[:, 1]) == NullableVector{Int}
    @test typeof(dt[:, 2]) == NullableVector{Float64}
    @test typeof(dt[:, 3]) == NullableVector{String}
    @test all(isnull, dt[:, 1])
    @test all(isnull, dt[:, 2])
    @test all(isnull, dt[:, 3])

    dt = DataTable([Nullable{Int}, Nullable{Float64}, Nullable{String}],
                   [:A, :B, :C], 100)
    @test size(dt, 1) == 100
    @test size(dt, 2) == 3
    @test typeof(dt[:, 1]) == NullableVector{Int}
    @test typeof(dt[:, 2]) == NullableVector{Float64}
    @test typeof(dt[:, 3]) == NullableVector{String}
    @test all(isnull, dt[:, 1])
    @test all(isnull, dt[:, 2])
    @test all(isnull, dt[:, 3])

    dt = DataTable([Nullable{Int}, Nullable{Float64}, Nullable{String}],
                   [:A, :B, :C], [false, false, true],100)
    @test size(dt, 1) == 100
    @test size(dt, 2) == 3
    @test typeof(dt[:, 1]) == NullableVector{Int}
    @test typeof(dt[:, 2]) == NullableVector{Float64}
    @test typeof(dt[:, 3]) == NullableCategoricalVector{Compat.UTF8String,UInt32}
    @test all(isnull, dt[:, 1])
    @test all(isnull, dt[:, 2])
    @test all(isnull, dt[:, 3])

    dt = convert(DataTable, zeros(10, 5))
    @test size(dt, 1) == 10
    @test size(dt, 2) == 5
    @test typeof(dt[:, 1]) == Vector{Float64}

    dt = convert(DataTable, ones(10, 5))
    @test size(dt, 1) == 10
    @test size(dt, 2) == 5
    @test typeof(dt[:, 1]) == Vector{Float64}

    dt = convert(DataTable, eye(10, 5))
    @test size(dt, 1) == 10
    @test size(dt, 2) == 5
    @test typeof(dt[:, 1]) == Vector{Float64}

    @test DataTable(NullableArray[[1,2,3],[2.5,4.5,6.5]], [:A, :B]) ==
        DataTable(A = NullableArray([1,2,3]), B = NullableArray([2.5,4.5,6.5]))

    # This assignment was missing before
    dt = DataTable(Column = [:A])
    dt[1, :Column] = "Testing"

    # zero-row datatable and subdatatable test
    dt = DataTable(x=[], y=[])
    @test nrow(dt) == 0
    dt = DataTable(x=[1:3;], y=[3:5;])
    sdt = view(dt, dt[:x] .== 4)
    @test size(sdt, 1) == 0

    @test hash(convert(DataTable, [1 2; 3 4])) == hash(convert(DataTable, [1 2; 3 4]))
    @test hash(convert(DataTable, [1 2; 3 4])) != hash(convert(DataTable, [1 3; 2 4]))

    # push!(dt, row)
    dt=DataTable( first=[1,2,3], second=["apple","orange","pear"] )

    dtb= DataTable( first=[1,2], second=["apple","orange"] )
    push!(dtb, Any[3,"pear"])
    @test isequal(dt, dtb)

    dtb= DataTable( first=[1,2], second=["apple","orange"] )
    push!(dtb, (3,"pear"))
    @test isequal(dt, dtb)

    dtb= DataTable( first=[1,2], second=["apple","orange"] )
    @test_throws ArgumentError push!(dtb, (33.33,"pear"))

    dtb= DataTable( first=[1,2], second=["apple","orange"] )
    @test_throws ArgumentError push!(dtb, ("coconut",22))

    dtb= DataTable( first=[1,2], second=["apple","orange"] )
    push!(dtb, @compat(Dict(:first=>3, :second=>"pear")))
    @test isequal(dt, dtb)

    dt=DataTable( first=[1,2,3], second=["apple","orange","banana"] )
    dtb= DataTable( first=[1,2], second=["apple","orange"] )
    push!(dtb, @compat(Dict("first"=>3, "second"=>"banana")))
    @test isequal(dt, dtb)

    dt0= DataTable( first=[1,2], second=["apple","orange"] )
    dtb= DataTable( first=[1,2], second=["apple","orange"] )
    @test_throws ArgumentError push!(dtb, @compat(Dict(:first=>true, :second=>false)))
    @test isequal(dt0, dtb)

    dt0= DataTable( first=[1,2], second=["apple","orange"] )
    dtb= DataTable( first=[1,2], second=["apple","orange"] )
    @test_throws ArgumentError push!(dtb, @compat(Dict("first"=>"chicken", "second"=>"stuff")))
    @test isequal(dt0, dtb)

    # delete!
    dt = DataTable(a=1, b=2, c=3, d=4, e=5)
    @test_throws ArgumentError delete!(dt, 0)
    @test_throws ArgumentError delete!(dt, 6)
    @test_throws KeyError delete!(dt, :f)

    d = copy(dt)
    delete!(d, [:a, :e, :c])
    @test names(d) == [:b, :d]
    delete!(d, :b)
    @test isequal(d, dt[[:d]])

    d = copy(dt)
    delete!(d, [2, 5, 3])
    @test names(d) == [:a, :d]
    delete!(d, 2)
    @test isequal(d, dt[[:a]])

    # deleterows!
    dt = DataTable(a=[1, 2], b=[3., 4.])
    @test deleterows!(dt, 1) === dt
    @test isequal(dt, DataTable(a=[2], b=[4.]))

    dt = DataTable(a=[1, 2], b=[3., 4.])
    @test deleterows!(dt, 2) === dt
    @test isequal(dt, DataTable(a=[1], b=[3.]))

    dt = DataTable(a=[1, 2, 3], b=[3., 4., 5.])
    @test deleterows!(dt, 2:3) === dt
    @test isequal(dt, DataTable(a=[1], b=[3.]))

    dt = DataTable(a=[1, 2, 3], b=[3., 4., 5.])
    @test deleterows!(dt, [2, 3]) === dt
    @test isequal(dt, DataTable(a=[1], b=[3.]))

    dt = DataTable(a=NullableArray([1, 2]), b=NullableArray([3., 4.]))
    @test deleterows!(dt, 1) === dt
    @test isequal(dt, DataTable(a=NullableArray([2]), b=NullableArray([4.])))

    dt = DataTable(a=NullableArray([1, 2]), b=NullableArray([3., 4.]))
    @test deleterows!(dt, 2) === dt
    @test isequal(dt, DataTable(a=NullableArray([1]), b=NullableArray([3.])))

    dt = DataTable(a=NullableArray([1, 2, 3]), b=NullableArray([3., 4., 5.]))
    @test deleterows!(dt, 2:3) === dt
    @test isequal(dt, DataTable(a=NullableArray([1]), b=NullableArray([3.])))

    dt = DataTable(a=NullableArray([1, 2, 3]), b=NullableArray([3., 4., 5.]))
    @test deleterows!(dt, [2, 3]) === dt
    @test isequal(dt, DataTable(a=NullableArray([1]), b=NullableArray([3.])))

    # describe
    #suppress output and test that describe() does not throw
    devnull = is_unix() ? "/dev/null" : "nul"
    open(devnull, "w") do f
        @test nothing == describe(f, DataTable(a=[1, 2], b=Any["3", Nullable()]))
        @test nothing ==
              describe(f, DataTable(a=NullableArray([1, 2]),
                                    b=NullableArray(Nullable{String}["3", Nullable()])))
        @test nothing ==
              describe(f, DataTable(a=CategoricalArray([1, 2]),
                                    b=NullableCategoricalArray(Nullable{String}["3", Nullable()])))
        @test nothing == describe(f, [1, 2, 3])
        @test nothing == describe(f, NullableArray([1, 2, 3]))
        @test nothing == describe(f, CategoricalArray([1, 2, 3]))
        @test nothing == describe(f, Any["1", "2", Nullable()])
        @test nothing == describe(f, NullableArray(Nullable{String}["1", "2", Nullable()]))
        @test nothing == describe(f, NullableCategoricalArray(Nullable{String}["1", "2", Nullable()]))
    end

    #Check the output of unstack
    dt = DataTable(Fish = NullableCategoricalArray(["Bob", "Bob", "Batman", "Batman"]),
                   Key = Nullable{String}["Mass", "Color", "Mass", "Color"],
                   Value = Nullable{String}["12 g", "Red", "18 g", "Grey"])
    # Check that reordering levels does not confuse unstack
    levels!(dt[1], ["XXX", "Bob", "Batman"])
    #Unstack specifying a row column
    dt2 = unstack(dt,:Fish, :Key, :Value)
    #Unstack without specifying a row column
    dt3 = unstack(dt,:Key, :Value)
    #The expected output
    dt4 = DataTable(Fish = Nullable{String}["XXX", "Bob", "Batman"],
                    Color = Nullable{String}[Nullable(), "Red", "Grey"],
                    Mass = Nullable{String}[Nullable(), "12 g", "18 g"])
    @test isequal(dt2, dt4)
    @test typeof(dt2[:Fish]) <: NullableCategoricalArray{String,1,UInt32}
    # first column stays as CategoricalArray in dt3
    @test isequal(dt3[:, 2:3], dt4[2:3, 2:3])
    #Make sure unstack works with NULLs at the start of the value column
    dt[1,:Value] = Nullable()
    dt2 = unstack(dt,:Fish, :Key, :Value)
    #This changes the expected result
    dt4[2,:Mass] = Nullable()
    @test isequal(dt2, dt4)

    dt = DataTable(A = 1:10, B = 'A':'J')
    @test !(dt[:,:] === dt)

    @test append!(DataTable(A = 1:2, B = 1:2), DataTable(A = 3:4, B = 3:4)) == DataTable(A=1:4, B = 1:4)
    dt = DataTable(A = NullableArray(1:3), B = NullableArray(4:6))
    @test all(c -> isa(c, NullableArray), categorical!(deepcopy(dt)).columns)
    @test all(c -> isa(c, NullableCategoricalArray), categorical!(deepcopy(dt), [1,2]).columns)
    @test all(c -> isa(c, NullableCategoricalArray), categorical!(deepcopy(dt), [:A,:B]).columns)
    @test find(c -> isa(c, NullableCategoricalArray), categorical!(deepcopy(dt), [:A]).columns) == [1]
    @test find(c -> isa(c, NullableCategoricalArray), categorical!(deepcopy(dt), :A).columns) == [1]
    @test find(c -> isa(c, NullableCategoricalArray), categorical!(deepcopy(dt), [1]).columns) == [1]
    @test find(c -> isa(c, NullableCategoricalArray), categorical!(deepcopy(dt), 1).columns) == [1]

    @testset "unstack nullable promotion" begin
        dt = DataTable(Any[repeat(1:2, inner=4), repeat('a':'d', outer=2), collect(1:8)],
                       [:id, :variable, :value])
        udt = unstack(dt)
        @test udt == unstack(dt, :variable, :value) == unstack(dt, :id, :variable, :value)
        @test udt == DataTable(Any[Nullable[1, 2], Nullable[1, 5], Nullable[2, 6],
                                   Nullable[3, 7], Nullable[4, 8]], [:id, :a, :b, :c, :d])
        @test all(typeof.(udt.columns) .== NullableVector{Int})
        dt = DataTable(Any[categorical(repeat(1:2, inner=4)),
                           categorical(repeat('a':'d', outer=2)), categorical(1:8)],
                       [:id, :variable, :value])
        udt = unstack(dt)
        @test udt == unstack(dt, :variable, :value) == unstack(dt, :id, :variable, :value)
        @test udt == DataTable(Any[Nullable[1, 2], Nullable[1, 5], Nullable[2, 6],
                                   Nullable[3, 7], Nullable[4, 8]], [:id, :a, :b, :c, :d])
        @test all(typeof.(udt.columns) .== NullableCategoricalVector{Int, UInt32})
    end

    @testset "duplicate entries in unstack warnings" begin
        dt = DataTable(id=NullableArray([1, 2, 1, 2]), variable=["a", "b", "a", "b"], value=[3, 4, 5, 6])
        @static if VERSION >= v"0.6.0-dev.1980"
            @test_warn "Duplicate entries in unstack." unstack(dt, :id, :variable, :value)
            @test_warn "Duplicate entries in unstack at row 3." unstack(dt, :variable, :value)
        end
        a = unstack(dt, :id, :variable, :value)
        b = unstack(dt, :variable, :value)
        @test a == b == DataTable(id = Nullable[1, 2], a = [5, Nullable()], b = [Nullable(), 6])

        dt = DataTable(id=NullableArray(1:2), variable=["a", "b"], value=3:4)
        @static if VERSION >= v"0.6.0-dev.1980"
            @test_nowarn unstack(dt, :id, :variable, :value)
            @test_nowarn unstack(dt, :variable, :value)
        end
        a = unstack(dt, :id, :variable, :value)
        b = unstack(dt, :variable, :value)
        @test a == b == DataTable(id = Nullable[1, 2], a = [3, Nullable()], b = [Nullable(), 4])
    end

    @testset "rename" begin
        dt = DataTable(A = 1:3, B = 'A':'C')
        @test names(rename(dt, :A, :A_1)) == [:A_1, :B]
        @test names(dt) == [:A, :B]
        @test names(rename!(dt, :A, :A_1)) == [:A_1, :B]
        @test names(dt) == [:A_1, :B]
    end

    @testset "size" begin
        dt = DataTable(A = 1:3, B = 'A':'C')
        @test_throws ArgumentError size(dt, 3)
        @test length(dt) == 2
        @test ndims(dt) == 2
    end

    @testset "description" begin
        dt = DataTable(A = 1:10)
        @test head(dt) == DataTable(A = 1:6)
        @test head(dt, 1) == DataTable(A = 1)
        @test tail(dt) == DataTable(A = 5:10)
        @test tail(dt, 1) == DataTable(A = 10)
    end

    @testset "misc" begin
        dt = DataTable(Any[collect('A':'C')])
        @test sprint(dump, dt) == """
                                  DataTables.DataTable  3 observations of 1 variables
                                    x1: Array{Char}((3,))
                                      1: Char A
                                      2: Char B
                                      3: Char C
                                  """
        dt = DataTable(A = 1:12, B = repeat('A':'C', inner=4))
        @test DataTables.without(dt, 1) == DataTable(B = repeat('A':'C', inner=4))
    end

    @testset "column conversions" begin
        dt = DataTable(Any[collect(1:10), collect(1:10)])
        @test !isa(dt[1], NullableArray)
        nullable!(dt, 1)
        @test isa(dt[1], NullableArray)
        @test !isa(dt[2], NullableArray)
        nullable!(dt, [1,2])
        @test isa(dt[1], NullableArray) && isa(dt[2], NullableArray)
    end
end

ext: jl


text:
function helper(f, meta, block, childsym)
    result = f(block)
    if result !== nothing
        return result
    end

    result = walk_and_replace(f, meta, getfield(block, childsym))
    if result !== nothing
        setfield!(block, childsym, result)
    end

    return nothing
end

walk_and_replace(f, meta, element) = f(element)

function walk_and_replace(f, meta, block::Vector)
    for (i, each) in enumerate(block)
        result = walk_and_replace(f, meta, each)
        if result !== nothing
            block[i] = result
        end
    end
end

const MDContentElements = Union{
    Markdown.BlockQuote,
    Markdown.Paragraph,
    Markdown.MD,
}
walk_and_replace(f, meta, block::MDContentElements) = helper(f, meta, block, :content)

const MDTextElements = Union{
    Markdown.Bold,
    Markdown.Header,
    Markdown.Italic,
}
walk_and_replace(f, meta, block::MDTextElements)      = helper(f, meta, block, :text)
walk_and_replace(f, meta, block::Markdown.Footnote)   = helper(f, meta, block, :text)
walk_and_replace(f, meta, block::Markdown.Admonition) = helper(f, meta, block, :content)
walk_and_replace(f, meta, block::Markdown.Image)      = helper(f, meta, block, :alt)
walk_and_replace(f, meta, block::Markdown.Table)      = helper(f, meta, block, :rows)
walk_and_replace(f, meta, block::Markdown.List)       = helper(f, meta, block, :items)
walk_and_replace(f, meta, block::Markdown.Link)       = helper(f, meta, block, :text)
walk_and_replace(f, meta, block::Documents.RawHTML)   = nothing
walk_and_replace(f, meta, block::Documents.DocsNodes) = helper(f, meta, block, :nodes)
walk_and_replace(f, meta, block::Documents.DocsNode)  = helper(f, meta, block, :docstr)
walk_and_replace(f, meta, block::Documents.EvalNode)  = helper(f, meta, block, :result)
walk_and_replace(f, meta, block::Documents.MetaNode)  = (merge!(meta, block.dict); nothing)
walk_and_replace(f, meta, block::Anchors.Anchor)      = helper(f, meta, block, :object)


text:
using NeXLCore

function Base.parse(
    ::Type{Material},
    expr::AbstractString;
    name::Union{AbstractString,Missing} = missing,
    properties::Dict{Symbol,Any} = Dict{Symbol,Any}(),
    atomicweights::Dict{Element,V} = Dict{Element,Float64}(),
    density::Union{Missing,AbstractFloat} = missing,
    description::Union{Missing,AbstractString} = missing,
    pedigree::Union{Missing,AbstractString} = missing,
    conductivity::Union{Missing,Symbol} = missing, # :Conductor, :Semiconductor, :Insulator
    lookup::Function = s->nothing
)::Material where {V<:AbstractFloat}
    mf, pex = _mp_level1(expr, lookup=lookup)
    pex = replace(pex, "*"=>"⋅")
    result = material(
        ismissing(name) ? pex : name,
        mf;
        properties = properties,
        atomicweights = atomicweights,
        density = density,
        description = description,
        conductivity = conductivity,
        pedigree = pedigree,
    )
    result[:Formula]=expr
    return result
end

macro mat_str(str)
    parse(Material, str)
end

"""
Parses "XXX+YYY+ZZZ" into add("XXX", "YYY", "ZZZ")
"""
function _mp_level1(expr::AbstractString; lookup::Function=s->nothing)
    plusuv(a::AbstractFloat, b::AbstractFloat) = a + b
    plusuv(a::UncertainValue, b::AbstractFloat) = NeXLUncertainties.sum(a,uv(b))
    plusuv(a::AbstractFloat, b::UncertainValue) = NeXLUncertainties.sum(uv(a),b)
    plusuv(a::UncertainValue, b::UncertainValue) = NeXLUncertainties.sum(a,b)
    function plus(d1,d2)
        d3 = Dict( elm => plusuv(get(d1[1],elm,0.0), get(d2[1],elm,0.0)) for elm in union(keys(d1[1]), keys(d2[1])))
        return (d3, "$(d1[2])+$(d2[2])")
    end
    function f(s) # Lookup while maintaining the raw representation
        # Julia 1.6 fails when this is reduced to a single call to replace
        t = reduce((a,b)->replace(a, b), ( "₀"=>"0", "₁"=>"1", "₂"=>"2", "₃"=>"3", "₄"=>"4", "₅"=>"5", "₆"=>"6", "₇"=>"7", "₈"=>"8", "₉"=>"9" ), init=string(s))
        return ( _mp_level2(t, s, lookup), s) # Finally parse expression
    end
    return mapreduce(f, plus, split(expr, "+"))
end

"""
Parses "#EXPR" or "#*EXPR" into times(#, "EXPR")
"""
function _mp_level2(expr::AbstractString, raw::AbstractString, lookup::Function)
    mult(n,d) = Dict(elm => n*q for (elm, q) in d)
    rfp = r"^(\d+(?:[.]\d*)?|[.]\d+)\s*[\*|⋅]?\s*(.*)$"
    ruv = r"^(\(\s*(?:\d+(?:[.]\d*)?|[.]\d+)\s*[±|+-|-+]\s*(?:\d+(?:[.]\d*)?|[.]\d+)\s*\))\s*[\*|⋅]\s*(.*)$"
    function luor3(s, raw)
        r=lookup(s)  # Look up is in mass fractions
        isnothing(r) && (r = lookup(raw))
        if isnothing(r)
            rt = _mp_level3(s)
            # Convert to mass fractions
            n = sum(a(elm) * value(v) for (elm, v) in rt)
            r = Dict( elm => (a(elm) / n) * v for (elm, v) in rt)
        end
        return r
    end
    muv = match(ruv, strip(expr))
    # Is there an UncertainValue multiplier
    if !isnothing(muv)
        k = parse(UncertainValue, strip(muv[1][2:end-1]))
        r = luor3(muv[2], raw)
    else
        # Is there an float multiplier
        mfp = match(rfp, expr)
        if !isnothing(mfp)
            k = parse(Float64, mfp[1])
            r = luor3(mfp[2], raw)
        else
            # Nope just parse the expression
            k = 1.0
            r = luor3(expr, raw)
        end
    end
    return mult(k, r)
end

"""
Parses "something⋅#H2O" or "something⋅#OH" or "something.PO4" or "something"
"""
function _mp_level3(expr::AbstractString)::Dict{Element, Int}
    plus(d1,d2) = Dict{Element,Int}( elm => get(d1,elm,0.0) + get(d2,elm,0.0) for elm in union(keys(d1), keys(d2)))
    mult(n,d) = Dict{Element,Int}(elm => n*q for (elm, q) in d)
    # Check for waters and other items postpended with an ⋅, · or .
    rjoined=r"^(.*)(⋅|·|\.)(\d*)(.*)$" # 1=>pre 2=>[⋅|.|·] 3=># 4=>post
    m = match(rjoined, expr)
    if !isnothing(m)
        res1=_mp_level4(m[1])
        n = length(m[3]) > 0 ? parse(Int,m[3]) : 1
        res2=_mp_level3(m[4])
        return plus(res1, mult(n, res2))
    end
    return _mp_level4(expr)
end

"""
Match parenthesis like "pre(inner)[n](post)"
"""
function _mp_level4(expr::AbstractString)
    mult(n,d) = Dict{Element,Int}(elm => n*q for (elm, q) in d)
    plus(d1,d2) = Dict{Element,Int}( elm => get(d1,elm,0.0) + get(d2,elm,0.0) for elm in union(keys(d1), keys(d2)))
    # Match parenthesis
    cx, start, stop, idx = 0, 0, 0, 0
    for _ in 1:length(expr)
        idx = nextind(expr, idx, 1)
        if expr[idx]=='('
            (start==0) && (start=idx)
            cx+=1
        end
        (expr[idx]==')') && (cx-=1)
        (cx<0) && error("Spurious ')` in $expr")
        (start!=0) && (cx==0) && (stop=idx)
        if (start!=0) && (stop!=0)
            inner = _mp_level4(expr[start+1:stop-1])
            pre = _mp_level5(expr[1:start-1])
            m = match(r"(\d*)(.*)", expr[stop+1:end])
            if !isnothing(m)
                n, post = length(m[1])>0 ? parse(Int,m[1]) : 1, _mp_level4(m[2])
            else
                n, post = 1, Dict{Element,Int}()
            end
            return plus(plus(pre, mult(n, inner)), post)
        end
    end
    (cx!=0) && error("Mismatched parenthesis in $expr")
    return _mp_level5(expr)
end

"""
Evaluate simple expressions like Al2O3 or SiO2 etc.
"""
function _mp_level5(expr::AbstractString)::Dict{Element,Int}
    function eval(res, s)
        relm=r"^(He|Li|Be|Ne|Na|Mg|Al|Si|Cl|Ar|Ca|Sc|Ti|Cr|Mn|Fe|Co|Ni|Cu|Zn|Ga|Ge|As|Se|Br|Kr|Rb|Sr|Zr|Nb|Mo|Tc|Ru|Rh|Pd|Ag|Cd|In|Sn|Sb|Te|Xe|Cs|Ba|La|Ce|Pr|Nd|Pm|Sm|Eu|Gd|Tb|Dy|Ho|Er|Tm|Yb|Lu|Hf|Ta|Re|Os|Ir|Pt|Au|Hg|Tl|Pb|Bi|Po|At|Rn|Fr|Ra|Ac|Th|Pa|Np|Pu|Am|Cm|Bk|Cf|Es|Fm|H|B|C|N|O|F|P|S|K|V|Y|I|W|U)(\d*)(.*)"
        m = match(relm,s)
        isnothing(m) && error("Unable to parse $expr as a simple formula.")
        elm=parse(Element,m[1])
        res[elm] = get(res,elm,0) + (length(m[2])==0 ? 1 : parse(Int,m[2]))
        return length(m[3])>0 ? eval(res, m[3]) : res
    end
    res = Dict{Element,Int}()
    return length(expr) > 0 ? eval(res, expr) : res
end


text:
@testset "reconstruction_matrix.jl" begin
    a = 0.0
    b = 2π
    N = 100
    M = N
    Δx = (b - a) / N

    # Domains
    closed_interval = ClosedIntervalDomain(a, b)
    periodic_interval = PeriodicIntervalDomain(a, b)
    unknown_domain = UnknownDomain()

    ## Identity filter
    f = IdentityFilter()

    # ClosedIntervalDomain
    R = reconstruction_matrix(f, closed_interval, M, N)
    @test R isa SparseMatrixCSC
    @test size(R) == (N + 1, M + 1)
    @test all(sum(R, dims = 2) .≈ 1)

    # PeriodicIntervalDomain
    R = reconstruction_matrix(f, periodic_interval, M, N)
    @test R isa SparseMatrixCSC
    @test size(R) == (N, M)
    @test all(sum(R, dims = 2) .≈ 1)

    # Unknown domain
    @test_throws Exception reconstruction_matrix(f, unknown_domain, M, N)

    ## Top hat filter
    h₀ = 0.05
    h = x -> h₀ * (1 - 1 / 2 * cos(x))
    f = TopHatFilter(h)

    # ClosedIntervalDomain
    R = reconstruction_matrix(f, closed_interval, M, N)
    @test R isa SparseMatrixCSC
    @test size(R) == (N + 1, M + 1)
    @test all(sum(R, dims = 2) .≈ 1)

    # PeriodicIntervalDomain
    R = reconstruction_matrix(f, periodic_interval, M, N)
    @test R isa SparseMatrixCSC
    @test size(R) == (N, M)
    @test all(sum(R, dims = 2) .≈ 1)

    # Unknown domain
    @test_throws Exception reconstruction_matrix(f, unknown_domain, M, N)

    ## Convolutional filter
    σ = Δx
    g = GaussianFilter(σ)

    # ClosedIntervalDomain
    R = reconstruction_matrix(g, closed_interval, M, N)
    @test R isa SparseMatrixCSC
    @test size(R) == (N + 1, M + 1)
    @test all(sum(R, dims = 2) .≈ 1)

    # PeriodicIntervalDomain
    R = reconstruction_matrix(g, periodic_interval, M, N)
    @test R isa SparseMatrixCSC
    @test size(R) == (N, M)
    @test all(sum(R, dims = 2) .≈ 1)

    # Unknown domain
    @test_throws Exception reconstruction_matrix(f, unknown_domain, M, N)

    ## Constant mesh-wide top hat filter
    h = x -> Δx / 2
    f = TopHatFilter(h)

    # ClosedIntervalDomain
    R = reconstruction_matrix_meshwidth(f, closed_interval, N)
    @test R isa SparseMatrixCSC
    @test size(R) == (N + 1, N + 1)
    @test all(sum(R, dims = 2) .≈ 1)

    # PeriodicIntervalDomain
    R = reconstruction_matrix_meshwidth(f, periodic_interval, N)
    @test R isa SparseMatrixCSC
    @test size(R) == (N, N)
    @test all(sum(R, dims = 2) .≈ 1)

    # Unknown domain
    @test_throws Exception reconstruction_matrix_meshwidth(f, unknown_domain, N)
end


text:
using HTTP
using JSON
using Base.Test

push!(LOAD_PATH, @__DIR__)
using Check

const json_headers = Dict{String, String}([
    "Accept" => "application/vnd.travis-ci.2+json",
    "Travis-API-Version" = "2",
    ])
const textplain_headers = Dict{String, String}([
    "Accept" => "text/plain",
    "Travis-API-Version" = "2",
    ])
const uriroot = "https://api.travis-ci.org"
const readtimeout = 60
const repo = "/repos/JuliaLang/julia"
const TravisDateFmt = dateformat"yyyy-mm-dd\THH:MM:SS\Z"

do_request(path) = HTTP.get(uriroot * path, headers = json_headers, readtimeout = readtimeout)

test_connection() = @testset "Sanity check" begin
    r = do_request("/")
    @test r.status== 200
    @test HTTP.headers(r)["Content-Type"] == "application/json"
    @test String(r.body) == """{"hello":"world"}"""
    r
end

struct IdIndex
    builds::Dict{Int, Dict{String, Any}}
    commits::Dict{Int, Dict{String, Any}}
    jobs::Dict{Int, Dict{String, Any}}
end
IdIndex() = IdIndex(
    Dict{Int, Dict{String, Any}}(),
    Dict{Int, Dict{String, Any}}(),
    Dict{Int, Dict{String, Any}}())
IdIndex(builds, other::IdIndex) = IdIndex(builds, other.commits, other.jobs)

function make_date(dict, elem)
    datestr = dict[elem]
    if isa(datestr, String) # it's Nullable!
        dict[elem] = parse(DateTime, datestr, TravisDateFmt)
    end
    nothing
end

function get_branch_builds!(state::IdIndex, after::String)
    if isempty(after)
        r = do_request("$repo/builds?event_type=push")
    else
        r = do_request("$repo/builds?event_type=push&after_number=$after")
    end
    @check r.status == 200

    data = JSON.parse(String(r.body))
    commits = data["commits"]
    for commit in commits
        commit_id = commit["id"]::Int
        make_date(commit, "committed_at")
        state.commits[commit_id] = commit
    end
    builds = data["builds"]
    build_id = nothing
    for build in builds
        build_id = build["id"]::Int
        make_date(build, "started_at")
        make_date(build, "finished_at")
        state.builds[build_id] = build
    end
    return build_id
end

function get_branch_builds_until!(state::IdIndex, before::DateTime)
    after = ""
    while true
        after_id = get_branch_builds!(state, after)
        if after_id === nothing || state.builds[after_id]["started_at"] < before
            return
        end
        after = state.builds[after_id]["number"]
    end
end

function filter_branch(state::IdIndex, branch::String)
    builds = filter(state.builds) do (id, build)
        return state.commits[build["commit_id"]]["branch"] == branch
    end
    return IdIndex(builds, state)
end

function map_build(f, state::IdIndex)
    newdict = map(state.builds) do build
        return build.first => f(build.second)
    end
    return newdict
end

function foreach_build(f, state::IdIndex)
    foreach(state.builds) do build
        f(build.second)
    end
end

function fetch_jobs!(state)
    foreach_build(state) do build
        foreach(build["job_ids"]) do job_id
            get!(state.jobs, job_id) do
                job = do_request("/jobs/$job_id")
                @check job.status == 200
                data = JSON.parse(String(job.body))
                return data["job"]
            end
        end
    end
end

const blacklist = Set([338558293, 339157435]) # job_id that are broken (406) on Travis
function fetch_sysimg_size!(state)
    foreach_build(state) do build
        foreach(build["job_ids"]) do job_id
            if !haskey(build, "bytes-linux-x86_64")
                job = state.jobs[job_id]
                job_config = job["config"]
                if job_config["os"] == "linux" && contains(job_config["env"], "x86_64")
                    if job["state"] == "passed" && !in(job_id, blacklist)
                        raw_log = HTTP.get("$uriroot/jobs/$job_id/log", headers = textplain_headers)
                        @check raw_log.status == 200
                        s = String(raw_log.body)
                        m = match(r"\.data\s+(\d+)", s)
                        @check isa(m, RegexMatch)
                        bytes = parse(Int, m[1])
                    else
                        bytes = 0
                    end
                    build["bytes-linux-x86_64"] = bytes
                end
            end
        end
    end
end

function build_stats_bytes(state)
    min_bytes = typemax(Int)
    max_bytes = 0
    foreach_build(state) do build
        bytes = build["bytes-linux-x86_64"]::Int
        if bytes > 0
            min_bytes = min(bytes, min_bytes)
            max_bytes = max(bytes, max_bytes)
        end
    end
    return (min_bytes, max_bytes)
end

if false
    # ulimit -n 8000
    io = STDOUT
    cols = displaysize(io)[2]
    state = IdIndex()
    until = now(Dates.UTC) - Dates.Day(30) # last 30 days
    get_branch_builds_until!(state, until)
    state = filter_branch(state, "master")
    fetch_jobs!(state)
    fetch_sysimg_size!(state)
    min_bytes, max_bytes = build_stats_bytes(state)
    for id in sort(collect(keys(state.builds)))
        bytes = state.builds[id]["bytes-linux-x86_64"]
        commit = state.commits[state.builds[id]["commit_id"]]
        sha = commit["sha"][1:10]
        msg = split(commit["message"], '\n', limit=2)[1][1:min(end, 40)]
        spark = "#" ^ (bytes > 0 ? 1 + cols ÷ 2 * (bytes - min_bytes) ÷ (max_bytes - min_bytes) : 0)
        println("$sha $bytes bytes  $spark    \t$msg")
    end
end

max_stars_repo_name: vtjnash/PerfPlots

max_issues_repo_name: vtjnash/PerfPlots

max_forks_repo_name: vtjnash/PerfPlots


text:
_maybe_apply(f::Function, x, p, t) = f(x, p, t)
_maybe_apply(f, x, p, t) = f

"""
    f = apply_inputs(func; kwargs...)

Creates functions for use in ODEProblems that allow functions of type ```(x,p,t)->{...}```
    to be passed in through keyword arguments and applied before passing their result into
    the main ODE function. Creates methods for both in-place and out-of-place use.

# Examples

```jldoctest
julia> using ControlSimTools, OrdinaryDiffEq

julia> function eom!(dx, x, m, t; F=0) # <--Note the input F
       dx[1] = x[2]
       dx[2] = F/m
       return nothing
       end
eom! (generic function with 1 method)

julia> prob = ODEProblem(apply_inputs(eom!, F=(x,p,t)->sin(t)), [0.0, 0.0], (0.0, 10.0), 100.0);

julia> sol = solve(prob, Tsit5());

julia> sol.u[end]
2-element Array{Float64,1}:
 0.105440211435822
 0.01839060375124582
```
"""
function apply_inputs(func; kwargs...)
    simfun(dx, x, p, t) = func(dx, x, p, t; map(f->_maybe_apply(f, x, p, t), (;kwargs...))...)
    simfun(x, p, t) = func(x, p, t; map(f->_maybe_apply(f, x, p, t), (;kwargs...))...)
    return simfun
end


text:
using Images, FileIO, Netpbm

rgbimg = load("data/bitmapInputTest.ppm")
greyimg = Gray.(rgbimg)
save("data/bitmapOutputTest.ppm", greyimg)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Test of some schemes

using Test
using StaticArrays
using FiniteVolumes
using FiniteVolumes.CartesianMeshes: Half, dx

@testset "Schemes" begin

@testset "Upwind linear advection" begin
    scheme = Upwind()
    @test scheme isa FiniteVolumes.Scheme

    # 1D
    mesh = CartesianMesh(0.0, 1.0, 2)
    left_flux = LinearAdvectionFlux(1.0)
    right_flux = LinearAdvectionFlux(-1.0)

    middle_face = (Half(3),)
    @test scheme(left_flux, mesh, [1.0, 2.0], middle_face) == 1.0
    @test scheme(right_flux, mesh, [1.0, 2.0], middle_face) == -2.0

    # Upwind scheme should ignore the time step passed as 5th argument
    @test scheme(left_flux, mesh, [1.0, 2.0], middle_face, 0.2) == 1.0
    @test scheme(left_flux, mesh, [1.0, 2.0], middle_face, nothing) == 1.0

    # Other types from StaticArrays
    @test scheme(left_flux, mesh, [Scalar(1.0), Scalar(2.0)], middle_face) == Scalar(1.0)
    @test scheme(left_flux, mesh, [SVector(1.0), SVector(2.0)], middle_face) == SVector(1.0)
    @test scheme(left_flux, mesh, [SVector(1.0, 2.0), SVector(2.0, 2.0)], middle_face) == SVector(1.0, 2.0)
    @test scheme(right_flux, mesh, [Scalar(1.0), Scalar(2.0)], middle_face) == -Scalar(2.0)
    @test scheme(right_flux, mesh, [SVector(1.0), SVector(2.0)], middle_face) == -SVector(2.0)
    @test scheme(right_flux, mesh, [SVector(1.0, 2.0), SVector(2.0, 2.0)], middle_face) == -SVector(2.0, 2.0)

    # 2D horizontal
    hmesh_2d = PeriodicCartesianMesh(2, 1)
    left_flux_2d = LinearAdvectionFlux([1.0, 0.0])
    right_flux_2d = LinearAdvectionFlux([-1.0, 0.0])

    middle_face = (Half(3), Half(2))
    @test scheme(left_flux_2d, hmesh_2d, [1.0, 2.0], middle_face) == 1.0
    @test scheme(left_flux_2d, hmesh_2d, [SVector(1.0), SVector(2.0)], middle_face) == SVector(1.0)
    @test scheme(left_flux_2d, hmesh_2d, [SVector(1.0, 1.0), SVector(2.0, 2.0)], middle_face) == SVector(1.0, 1.0)
    @test scheme(right_flux_2d, hmesh_2d, [1.0, 2.0], middle_face) == -2.0
    @test scheme(right_flux_2d, hmesh_2d, [SVector(1.0), SVector(2.0)], middle_face) == -SVector(2.0)
    @test scheme(right_flux_2d, hmesh_2d, [SVector(1.0, 1.0), SVector(2.0, 2.0)], middle_face) == -SVector(2.0, 2.0)

    # 2D vertical
    vmesh_2d = PeriodicCartesianMesh(1, 2)
    bottom_flux = LinearAdvectionFlux([0.0, 1.0])
    top_flux = LinearAdvectionFlux([0.0, -1.0])

    middle_face = (Half(2), Half(3))
    @test scheme(bottom_flux, vmesh_2d, [1.0 2.0], middle_face) == 1.0
    @test scheme(top_flux, vmesh_2d, [1.0 2.0], middle_face) == -2.0
end

@testset "Courant" begin
    mesh = CartesianMesh(2)
    left_flux = LinearAdvectionFlux(1.0)
    faster_left_flux = LinearAdvectionFlux(10.0)
    right_flux = LinearAdvectionFlux(-1.0)

    @test FiniteVolumes.courant(1.0, left_flux, mesh, [1.0, 1.0]) == 2.0
    @test FiniteVolumes.courant(1.0, faster_left_flux, mesh, [1.0, 1.0]) == 20.0
    @test FiniteVolumes.courant(1.0, right_flux, mesh, [1.0, 1.0]) == 2.0
end

@testset "Div" begin
    # 1D
    flux = LinearAdvectionFlux(2.0)
    mesh = PeriodicCartesianMesh(45)
    dt = 0.4*dx(mesh)[1]
    w0 = rand(FiniteVolumes.nb_cells(mesh))
    D = FiniteVolumes.inner_faces_to_cells_matrix(mesh)
    Δw1 = D * FiniteVolumes.numerical_fluxes(flux, mesh, w0, Upwind(), dt)
    Δw2 = -FiniteVolumes.div(flux, mesh, w0, Upwind(), dt)
    @test Δw1 ≈ Δw2

    # 2D
    flux = LinearAdvectionFlux([9.0, 4.0])
    mesh = PeriodicCartesianMesh(5, 8)
    dt = 0.4*dx(mesh)[1]
    w0 = rand(size(FiniteVolumes.cell_centers(mesh))...)
    D = FiniteVolumes.inner_faces_to_cells_matrix(mesh)
    Δw1 = D * FiniteVolumes.numerical_fluxes(flux, mesh, w0, Upwind(), dt)
    Δw2 = -FiniteVolumes.div(flux, mesh, w0, Upwind(), dt)
    @test reshape(Δw1, Tuple(mesh.nb_cells)) ≈ Δw2
end

end


text:
# Feynman/utils/graph_wrapper.jl
#
# MIT License
#
# Copyright (c) 2019 Brandon Gomes
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

import LightGraphs: AbstractGraph, is_directed, edgetype, inneighbors, outneighbors,
                    ne, nv, edges, vertices, has_edge, has_vertex

"""
    GraphWrapper{I}
        <: AbstractGraph{I}

Abstract Base Type for objects which contain and behave like a graph.
"""
abstract type GraphWrapper{I} <: AbstractGraph{I} end

"""
    graph(g::AbstractGraph)

Returns internal graph representation of graph.
"""
graph(g::AbstractGraph) = g

"""
    is_directed(g::GraphWrapper)
        === is_directed(graph(g))
"""
is_directed(g::GraphWrapper) = is_directed(graph(g))

"""
    edgetype(g::GraphWrapper)
        === edgetype(graph(g))
"""
edgetype(g::GraphWrapper) = edgetype(graph(g))

"""
    inneighbors(g::GraphWrapper, v)
        === inneighbors(graph(g), v)
"""
inneighbors(g::GraphWrapper, v) = inneighbors(graph(g), v)

"""
    outneighbors(g::GraphWrapper, v)
        === outneighbors(graph(g), v)
"""
outneighbors(g::GraphWrapper, v) = outneighbors(graph(g), v)

"""
    ne(g::GraphWrapper)
        === ne(graph(g))
"""
ne(g::GraphWrapper) = ne(graph(g))

"""
    nv(g::GraphWrapper)
        === nv(graph(g))
"""
nv(g::GraphWrapper) = nv(graph(g))

"""
    edges(g::GraphWrapper)
        === edges(graph(g))
"""
edges(g::GraphWrapper) = edges(graph(g))

"""
    vertices(g::GraphWrapper)
        === vertices(graph(g))
"""
vertices(g::GraphWrapper) = vertices(graph(g))

"""
    has_edge(g::GraphWrapper, args...)
        === has_edge(graph(g), args...)
"""
has_edge(g::GraphWrapper, args...) = has_edge(graph(g), args...)

"""
    has_vertex(g::GraphWrapper, args...)
        === has_vertex(graph(g), args...)
"""
has_vertex(g::GraphWrapper, args...) = has_vertex(graph(g), args...)


text:
struct GaussianProcess{M,K,L,O} <: AbstractGaussianProcess
    mean::M
    kernel::K
    lognoise::L
    optimize::Bool
    optimizer::O
end

"""
    GaussianProcess(mean, kernel;
        lognoise=-2,
        optimize=false,
        optimizer=Optim.LBFGS())

Wrapper around [GaussianProcesses](http://stor-i.github.io/GaussianProcesses.jl/latest/)
gives the possibility of training every model in the library using Gaussian Regression.
GaussianProcesses.jl needs to be called in order to use these models. The use of Gaussian
Regression for ESNs has first been explored in [1].

[1] Chatzis, Sotirios P., and Yiannis Demiris. "_Echo state Gaussian process._"
IEEE Transactions on Neural Networks 22.9 (2011): 1435-1445.
"""
function GaussianProcess(
    mean,
    kernel;
    lognoise=-2,
    optimize=false,
    optimizer=Optim.LBFGS()
)
    return GaussianProcess(mean, kernel, lognoise, optimize, optimizer)
end

function _train(states, target_data, gp::GaussianProcess)
    out_size = size(target_data, 1)
    output_matrix = []

    for i in 1:out_size
        out_size == 1 ? target = vec(target_data) : target = target_data[i,:]
        push!(output_matrix, GP(states, target, gp.mean, gp.kernel, gp.lognoise))
        gp.optimize ? optimize!(output_matrix[i]; method=gp.optimizer) : nothing
    end

    return OutputLayer(gp, output_matrix, out_size, target_data[:,end])
end

max_stars_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_issues_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_forks_repo_name: MartinuzziFrancesco/EchoStateNetwork


text:
using QuantEcon
using Statistics
using LinearAlgebra

include(string(pwd(),"/src/gensys.jl"))
include(string(pwd(),"/hmc/diffs_estrut_v3.jl"))
include(string(pwd(),"/hmc/diff_kalman.jl"))

function solve_lyapunov_vec(A,Q)
    aux = (I - kron(A,A))
    vec_ans = \(aux,vec(Q))
    res = reshape(vec_ans,size(Q))
    return res
end

function log_like_dsge(par,data;kalman_tol = 1e-10,reduce=true)
    #order to par
    #alfa
    #beta
    #epsilon
    #theta
    #sig
    #sigma: this is the std dev of the innovation
    #phi
    #phi_pi
    #phi_y
    #rho_v

    #data will have t x p dimension: lines are periods p are variables

    alfa = par[1]
    bet = par[2]
    epsilon = par[3]
    theta = par[4]
    sig = par[5]
    #par 6 is coded bellow see line 56
    phi = par[7]
    phi_pi = par[8]
    phi_y = par[9]
    rho_v = par[10]

    THETA = (1-alfa)/(1-alfa+alfa*epsilon)
    lamb = (1-theta)*(1-bet*theta)/theta*THETA
    kappa = lamb*(sig+(phi+alfa)/(1-alfa))

    nobs = size(data,1)
    l = size(data,2)

    GAMMA_0 = [bet  0    0  0;
               1    sig  0  0;
               0    0    0  0;
               0    0    0  1]

    GAMMA_1 = [ 1       -kappa  0   0;
                0        sig    1   0;
               -phi_pi  -phi_y  1  -1;
                0        0      0   rho_v]

    PSI = [0; 0; 0; par[6]]

    PI = [bet  0;
          1    sig;
          0    0;
          0    0]

    p = size(GAMMA_1,1) #number of endogenous vars

    sol = gensys(GAMMA_0,GAMMA_1,PSI,PI; verbose = false)

    if sum(sol.eu) != 2
        return -Inf, repeat([0],length(par))
    end

    #Sig = zeros(p,p)
    #Sig[4,4] = par[6]

    G = zeros(1,p)
    G[1,2] = 1

    A = sol.Theta1
    B = sol.Theta2
    A[abs.(A) .< 10*eps()] .= 0
    B[abs.(B) .< 10*eps()] .= 0
    R = [0] .+ 1e-2
    Q = sol.Theta2*sol.Theta2'
    #Q_eig = eigen(Q)
    #Q = Q_eig.vectors*diagm(repeat([mean(Q_eig.values)],p))*Q_eig.vectors'

    kalman_res = Kalman(A,G,Q,R) #create a Kalman filter instance

    y_mean = mean(data,dims=1)
    y_var = var(data,dims=1)
    y_var = reshape(y_var,size(y_var,2))
    y_var = repeat([1],p)

    x_hat = repeat(y_mean,p) #initial mean of the state
    x_var = I(4) #diagm(repeat(y_var,p))#variance initial of state
    #x_var = x_var*x_var'
    set_state!(kalman_res,x_hat,x_var)

    fit = zeros(nobs,4)

    dA,dB = diff_mod(par)
    #dB = dB[1:4,:]
    dA = dA'*I
    dG = zeros(10,4)
    dQ = (kron(sol.Theta2,I(p)) + kron(I(p),sol.Theta2)*commutation_matrix(4,1))*dB
    dQ = dQ'
    dR = zeros(size(par,1),1)

    llh = zeros(nobs)
    dll = zeros(nobs,size(par,1))

    dx0 = zeros(p)
    S0 = zeros(p,p)

    g_x_hat_l = zeros(size(par,1),p)
    P = kalman_res.G*S0*kalman_res.G' + kalman_res.R
    d_Sl = zeros(10,p^2)#transpose(diff_S0(A,S0,dA,dQ)) #vec, not vech

    for j in 1:nobs
        med = kalman_res.cur_x_hat
        fit[j,:] = med
        varian = kalman_res.cur_sigma
        #println(det(varian))
        eta = data[j,:] - kalman_res.G*med #mean loglike
        P = kalman_res.G*varian*kalman_res.G' + kalman_res.R
        teste_cond = 1/cond(P)

        gain = (A*varian*G')*inv(P)

        d_p = transpose(diff_P(G,varian,dG,d_Sl,dR))
        d_gain = transpose(diff_K(A,P,G,varian,dA,dG,d_p,d_Sl))
        g_eta = transpose(grad_y(med,g_x_hat_l,G,dG))
        g_x_hat = transpose(grad_x_hat(med,dA,A,g_x_hat_l,eta,g_eta, gain,d_gain))

        p_inv = inv(P)

        if teste_cond < kalman_tol
            llh[j] = -500
            dll[j,:] .= 0
        #elseif det(P) <= 0
        #    llh[j] = -500
        else
            llh[j] = -(p*log(2*pi) + logdet(P) .+ eta'*p_inv*eta)/2
            dll[j,:] = -1/2*p_inv*d_p' - (eta'*p_inv)*g_eta' + 1/2*kron(eta'*p_inv,eta'*p_inv)*d_p'
            QuantEcon.update!(kalman_res,data[j,:]) #updating the kalman estimates
        end #end if
        #println(kalman_res.cur_sigma)
        #println(det(varian))
        #println(P)
        g_x_hat_l = g_x_hat
        d_Sl = transpose(diff_S(A,P,gain,varian,dA,d_p,d_gain,d_Sl,dQ))
    end #end for
    llh = llh[10:length(llh)]
    dll = dll[10:length(llh),:]
    if reduce == true
        return sum(llh),sum(dll,dims=1)
    else
        return llh,dll
    end
end

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _rhotheta{T<:AbstractFloat}(period::T, periastron::T, eccentricity::T,
                                     semimajor_axis::T, inclination::T,
                                     omega::T, omega2::T, epoch::T)
    rho = theta = -one(period)
    # See chapter 55.
    n = 360.0/period
    M = deg2rad(n*(epoch - periastron))
    E = kepler_solver(M, eccentricity)
    r  = semimajor_axis*(1.0 - eccentricity*cos(E))
    nu = trueanom(E, eccentricity)
    # Convert variables in radians.
    omega2      = deg2rad(omega2)
    inclination = deg2rad(inclination)
    omega       = deg2rad(omega)
    theta = omega + atan2(sin(nu + omega2)*cos(inclination), cos(nu + omega2))
    rho   = r*cos(nu + omega2)/cos(theta - omega)
    # Convert theta to degrees and for it to be in [0, 360) range.
    theta = cirrange(rad2deg(theta))
    return rho, theta
end

"""
    rhotheta(period, periastron, eccentricity, semimajor_axis, inclination, omega, omega2, epoch) -> rho, theta

### Purpose ###

Calculate the separation and position angle of a binary star.

### Explanation ###

This function will return the separation \$\\rho\$ and position angle
\$\\theta\$ of a visual binary star derived from its orbital elements.  The
algorithms described in the following book will be used: Meeus J., 1992,
Astronomische Algorithmen, Barth.  Compared to the examples given at page 400
and no discrepancy found.

### Arguments ###

* `period`: period [year]
* `periastro`: time of periastron passage [year]
* `eccentricity`: eccentricity of the orbit
* `semimajor_axis`: semi-major axis [arc second]
* `inclination`: inclination angle [degree]
* `omega`: node [degree]
* `omega2`: longitude of periastron [degree]
* `epoch`: epoch of observation [year]

All input parameters have to be scalars.

### Output ###

The 2-tuple \$(\\rho, \\theta)\$, where

* \$\\rho\$: separation [arc second]
* \$\\theta\$: position angle [degree]

### Example ###

Find the position of Eta Coronae Borealis at the epoch 2016

``` julia
ρ, θ = rhotheta(41.623, 1934.008, 0.2763, 0.907, 59.025, 23.717, 219.907, 2016)
# => (0.6351167848228113,214.42513388052114)
```

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
rhotheta(period::Real, periastron::Real, eccentricity::Real,
         semimajor_axis::Real, inclination::Real,
         omega::Real, omega2::Real, epoch::Real) =
             _rhotheta(promote(float(period), float(periastron),
                               float(eccentricity), float(semimajor_axis),
                               float(inclination), float(omega),
                               float(omega2), float(epoch))...)


text:
module TestSKPreprocessing

using Random
using TSML
using TSMLextra
using PyCall
using Test

const IRIS = getiris()
const X = IRIS[:,1:4] |> DataFrame
const Y = IRIS[:,5] |> Vector

# "KernelCenterer","MissingIndicator","KBinsDiscretizer","OneHotEncoder",
const preprocessors = [
     "FastICA", "IncrementalPCA",
     "LatentDirichletAllocation",
     #"KernelPCA",
     #"DictionaryLearning",
     #"FactorAnalysis",
     #"MiniBatchDictionaryLearning",
     #"MiniBatchSparsePCA",
     #"NMF",
     #"PCA",
     #"TruncatedSVD",
     "VarianceThreshold",
     "SimpleImputer",
     "Binarizer", "FunctionTransformer",
     "MultiLabelBinarizer", "MaxAbsScaler", "MinMaxScaler", "Normalizer",
     "OrdinalEncoder", "PolynomialFeatures", "PowerTransformer",
     "QuantileTransformer", "RobustScaler", "StandardScaler"
 ]

function fit_test(preproc::String,in::DataFrame,out::Vector)
	_preproc=SKPreprocessor(Dict(:preprocessor=>preproc))
	fit!(_preproc,in,out)
	@test _preproc.model != nothing
	return _preproc
end

function transform_test(preproc::String,in::DataFrame,out::Vector)
	_preproc=SKPreprocessor(Dict(:preprocessor=>preproc))
	fit!(_preproc,in,out)
	res = transform!(_preproc,in)
	@test size(res)[1] == size(out)[1]
end

@testset "scikit preprocessors fit test" begin
    Random.seed!(123)
    for cl in preprocessors
	println(cl)
	fit_test(cl,X,Y)
    end
end

@testset "scikit preprocessors transform test" begin
    Random.seed!(123)
    for cl in preprocessors
	println(cl)
	transform_test(cl,X,Y)
    end
end

function skptest()
    features = X
    labels = Y

    pca = SKPreprocessor(Dict(:preprocessor=>"IncrementalPCA",:impl_args=>Dict(:n_components=>3)))
    fit!(pca,features)
    @test transform!(pca,features) |> x->size(x,2) == 3

    norml = SKPreprocessor(Dict(:preprocessor=>"Normalizer"))
    fit!(norml,features)
    @test transform!(norml,features) |> x->size(x,2) == 4

    ica = SKPreprocessor(Dict(:preprocessor=>"FastICA",:impl_args=>Dict(:n_components=>2)))
    fit!(ica,features)
    @test transform!(ica,features) |> x->size(x,2) == 2

    stdsc = SKPreprocessor(Dict(:preprocessor=>"StandardScaler",:impl_args=>Dict()))
    fit!(stdsc,features)
    @test abs(mean(transform!(stdsc,features) |> Matrix)) < 0.00001

    minmax = SKPreprocessor(Dict(:preprocessor=>"MinMaxScaler",:impl_args=>Dict()))
    fit!(minmax,features)
    @test mean(transform!(minmax,features) |> Matrix) ≈ 0.4486931104833648
end
@testset "scikit preprocessor fit/transform test with real data" begin
    Random.seed!(123)
    skptest()
end

end

max_stars_repo_name: ppalmes/TSMLextra

max_issues_repo_name: ppalmes/TSMLextra

max_forks_repo_name: ppalmes/TSMLextra


text:
# This file generates builtins.jl.

# Look up the expected number of arguments in Core.Compiler.tfunc data
function generate_fcall(f, table, id)
    if id !== nothing
        minarg, maxarg, tfunc = table[id]
    else
        minarg = 0
        maxarg = typemax(Int)
    end
    # The tfunc tables are wrong for fptoui and fptosi
    if f == "Base.fptoui" || f == "Base.fptosi"
        minarg = 2
    end
    # Generate a separate call for each number of arguments
    if maxarg < typemax(Int)
        wrapper = minarg == maxarg ? "" : "if nargs == "
        for nargs = minarg:maxarg
            if minarg < maxarg
                wrapper *= "$nargs\n            "
            end
            argcall = ""
            for i = 1:nargs
                argcall *= "@lookup(frame, args[$(i+1)])"
                if i < nargs
                    argcall *= ", "
                end
            end
            wrapper *= "return Some{Any}($f($argcall))"
            if nargs < maxarg
                wrapper *= "\n        elseif nargs == "
            end
        end
        if minarg < maxarg
            wrapper *= "\n        end"
        end
        return wrapper
    end
    # A built-in with arbitrary or unknown number of arguments.
    # This will (unfortunately) use dynamic dispatch.
    return "return Some{Any}($f(getargs(args, frame)...))"
end

# `io` is for the generated source file
# `intrinsicsfile` is the path to Julia's `src/intrinsics.h` file
function generate_builtins(io::IO)
    pat = r"(ADD_I|ALIAS)\((\w*),"
    print(io,
"""
# This file is generated by `generate_builtins.jl`. Do not edit by hand.

function getargs(args, frame)
    nargs = length(args)-1  # skip f
    callargs = resize!(frame.callargs, nargs)
    for i = 1:nargs
        callargs[i] = @lookup(frame, args[i+1])
    end
    return callargs
end

\"\"\"
    ret = maybe_evaluate_builtin(frame, call_expr)

If `call_expr` is to a builtin function, evaluate it, returning the result inside a `Some` wrapper.
Otherwise, return `call_expr`.
\"\"\"
function maybe_evaluate_builtin(frame, call_expr)
    # By having each call appearing statically in the "switch" block below,
    # each gets call-site optimized.

    args = call_expr.args
    nargs = length(args) - 1
    fex = args[1]
    if isa(fex, QuoteNode)
        f = fex.value
    else
        f = @lookup(frame, fex)
    end
    # Builtins and intrinsics have empty method tables. We can circumvent
    # a long "switch" check by looking for this.
    mt = typeof(f).name.mt
    if isa(mt, Core.MethodTable)
        isempty(mt) || return call_expr
    end
    # Builtins
""")
    firstcall = true
    for ft in subtypes(Core.Builtin)
        ft === Core.IntrinsicFunction && continue
        ft === getfield(Core, Symbol("#kw##invoke")) && continue  # handle this one later
        head = firstcall ? "if" : "elseif"
        firstcall = false
        f = ft.instance
        # Tuple is common, especially for returned values from calls. It's worth avoiding
        # dynamic dispatch through a call to `ntuple`.
        if f === tuple
            print(io,
"""
    $head f === $f
        return Some{Any}(ntuple(i->@lookup(frame, args[i+1]), length(args)-1))
""")
            continue
        end
        id = findfirst(isequal(f), Core.Compiler.T_FFUNC_KEY)
        fcall = generate_fcall(f, Core.Compiler.T_FFUNC_VAL, id)
        print(io,
"""
    $head f === $f
        $fcall
""")
        firstcall = false
    end
    print(io,
"""
    # Intrinsics
""")
    for fsym in names(Core.Intrinsics)
        fsym == :Intrinsics && continue
        isdefined(Base, fsym) || (println("skipping ", fname); continue)
        f = getfield(Base, fsym)
        f isa Core.IntrinsicFunction || error("not an intrinsic")
        id = reinterpret(Int32, f) + 1
        f = isdefined(Base, fsym) ? "Base.$fsym" :
            isdefined(Core, fsym) ? "Core.$fsym" : error("whoops on $f")
        fcall = generate_fcall(f, Core.Compiler.T_IFUNC, id)
        print(io,
"""
    elseif f === $f
        $fcall
""")
    end
    print(io,
"""
    end
    if isa(f, getfield(Core, Symbol("#kw##invoke")))
        return Some{Any}(getfield(Core, Symbol("#kw##invoke"))(getargs(args, frame)...))
    end
    return call_expr
end
""")
end


text:
tests = ["basicfuns", "rmath"]

for t in tests
    to_test = true
    # if VERSION < v"0.4.0-dev"
    #     to_test = false
    # end

    if to_test
        fp = "$t.jl"
        println("* running $fp")
        include(fp)
    end
end


text:
import Base.getindex,
	   LinQu.ITensors.getindex,
	   LinQu.ITensors.IndexVal

function Rotate(θ,nx,ny,nz)
	len = norm([nx,ny,nz])
	nx/=len
	ny/=len
	nz/=len
	return cos(θ/2)*I_DATA-1im*sin(θ/2)*(nx*X_DATA+ ny*Y_DATA+nz*Z_DATA)
end

U1 = VarGate(Rotate, [-0.75, 1,1,2])
U2 = VarGate(Rotate, [2.14, 1,2,1])

function test_circuit(N::Int)
	circuit= QCircuit(N)
	add!(circuit, U2(4))
	for i=1:4
		add!(circuit, U1(i))
	end
	add!(circuit, CNOT(1,2), CNOT(3,4))
	for i =1:4
		add!(circuit, U1(i))
	end
	add!(circuit, control(Z)(2,3))
	return circuit
end

function getindex(T::ITensor, inds::IndexSet, i::Vector{Int})
	if length(inds) != length(i)
		error("dimension of inds and i mismatch!")
	end
	ivs = Array{IndexVal,1}(undef, 0)
	for j = 1:length(i)
		push!(ivs, IndexVal(inds[j], i[j]))
	end
	return getindex(T, ivs...)
end

function print_prob(state::ExactState)
	N = order(state.site)
	inds = IndexSet()
	for i = 1:4
		push!(inds, findindex(state.site, "q=$(i)"))
	end
	for ii =0:2^N-1
	    index = [Int((ii&a)/a) for a∈ [1,2,4,8]]
	    amp = state.site[inds, (index.+1)]
	    print(index,"	>	", Real(amp*conj(amp)), "		", amp , "\n")
	end
end

function main()
	N=4
	mps_state = MPSState(N)
	exact_state = toExactState(mps_state)
	circuit = test_circuit(N)
	apply!(exact_state, circuit)
	apply!(mps_state, circuit)
	normalize!(mps_state)
	@test exact_state ≈ toExactState(mps_state)
	print_prob(exact_state)
end

main()


text:
module StressTest

"""
    dream(seconds)

Like Base.sleep() except maxes out the thread for a specified number of seconds. The minimum dream time is 1
millisecond or input of `0.001`.
"""
function dream(sec::Real)
    sec ≥ 0 || throw(ArgumentError("cannot dream for $sec seconds"))
    t = Timer(sec)
    while isopen(t)
        yield()
    end
    nothing
end

export dream

end # module


text:
using Test, TreesHeaps
import TreesHeaps: link!, cut!, getdir, opposite

test_show(x) = show(IOBuffer(), x)
test_showmime(x) = show(IOBuffer(), MIME{Symbol("text/plain")}(), x)

@testset "Constructors" begin
    @test SBN(1) == SBN(Int, 1.0)
    @test HBN(1) == HBN(Int, 1.0)
    @test RBN(1) == RBN(Int, 1.0)
    @test BST() == BST(Any)
    @test BST(0, 1) == insert!(BST(Int, 0.0), 1.0)
    @test BST(0, 1) == BST(Int, 0.0, 1.0)
    @test AVL() == AVL(Any)
    @test AVL(0, 1) == insert!(AVL(Int, 0.0), 1.0)
    @test Splay() == Splay(Any)
    @test Splay(0, 1) == insert!(Splay(Int, 0.0), 1.0)
    @test RBT() == RBT(Any)
    @test RBT(0, 1) == insert!(RBT(Int, 0.0), 1.0)
    @test RBT(0, 1) == RBT(Int, 0.0, 1.0)
end

@testset "Interface" begin
    s1 = BST(Float64, 1.0)
    insert!(s1, 0, 2)
    insert!(s1, -1, 3)
    nodes1 = [n for n in s1.root]
    s2 = BST(Float64, 1.0, height = true)
    insert!(s2, 0, 2)
    insert!(s2, -1, 3)
    nodes2 = [n for n in s2.root]
    @test eltype(s1) == eltype(s1.root)
    @test eltype(typeof(s1)) == eltype(typeof(s1.root))
    @test size(s2) == s2.size
    @test length(s2) == length(s2.root)
    @test height(s1.root) == 0
    @test first(nodes1) == s1.root.left
    @test last(nodes2) == s2.root.right
    x = collect(PreOrderDFS(s1))
    @test first(x) === s1.root
    @test last(x) === s1.root.right.right
    y = collect(PostOrderDFS(s2))
    @test first(y) === s2.root.left.left
    @test last(y) === s2.root
    z = [i * n.data for (i, n) in pairs(s1.root)]
    @test last(z) == 4
    @test s1.root == s2.root
    @test s1.root != NullNode()
    @test NullNode() != s2.root
    @test NullNode() == NullNode()
    @test s1 == s2
    @test findmax(NullNode()) == NullNode() == findmin(NullNode())
    test_show(s1)
    test_show(s2)
    test_showmime(s1)
    test_showmime(s2)
    test_show(s1.root)
    test_show(s2.root)
    test_showmime(s1.root)
    test_showmime(s2.root)
    test_show(NullNode())
    test_showmime(NullNode())
    plot_init()
end

@testset "Node operations" begin
    @test opposite(:left) == :right
    @test isnothing(opposite(nothing))
    @test !isnull(SBN(1))
    @test isnull(NullNode())
    @test height(SBN(0)) == 0
    @test height(HBN(0)) == 0
    @test height(NullNode()) == -1
    s = SBN(1)
    link!(s, SBN(2), :left)
    @test s.left.data == 2
    cut!(s, s.left)
    @test isnothing(s.left.data)
    link!(s, NullNode(), :left)
    @test isnothing(s.left.data)
    s = SBN(2)
    link!(NullNode(), s, :left)
    @test isnothing(s.parent.data)
    @test isnothing(cut!(NullNode(), s, :left))
    @test isnothing(cut!(s, NullNode(), :left))
    link!(NullNode(), s, nothing)
    @test isnothing(s.parent.data)
    @test isnothing(getdir(NullNode(), s))
    @test isnothing(getdir(s, NullNode()))
    @test isnothing(getdir(NullNode(), NullNode()))
    @test isnothing(link!(NullNode(), NullNode(), :left))
    @test isnothing(link!(NullNode(), NullNode(), nothing))
end

@testset "BinarySearchTree" begin
    s = BST(Float64)
    search(s, 1)
    search(s, 1, true)
    delete!(s, 1)
    delete!(s, 1, 2)
    search(s, 1.0)
    @test !first(search(s, 1.0, true))
    insert!(s, 1.0)
    insert!(s, 1.0, 0, 2)
    insert!(s, -1)
    insert!(s, 2.5)
    search(s, 1)
    search(s, 0.5)
    search(s, 1.5)
    @test first(search(s, 0, true))
    @test !first(search(s, 3, true))
    @test s.root.data == 1.0
    delete!(s, 0)
    delete!(s, -1.5)
    @test !first(search(s, 0, true))
    @test isnull(s.root.parent)
    x = collect(PostOrderDFS(s))
    @test last(x) === s.root
    @test first(x) === s.root.left
    @test size(s) == 4
    @test findmax(s).data == 2.5
    @test findmin(s).data == -1
    test_show(s)
    test_showmime(s)
end

@testset "BinarySearchTree with height" begin
    s = BST(1.0, height = true)
    search(s, 1)
    insert!(s, 0, 2)
    insert!(s, -1)
    insert!(s, 2.5)
    @test first(search(s, 0, true))
    @test !first(search(s, 3, true))
    @test s.root.data == 1.0
    delete!(s, 0)
    delete!(s, -1.5)
    @test !first(search(s, 0, true))
    @test isnull(s.root.parent)
    @test size(s) == 4
    @test length(s) == 2
    delete!(s, 2, 2.5)
    delete!(s, 1)
    @test s.root.data == -1
    insert!(s, 0)
    delete!(s, -1)
    @test s.root.data == 0
    delete!(s, 0)
    @test isnull(s.root)
    test_show(s)
    test_showmime(s)
end

@testset "AVLTree" begin
    s = AVL(Float64)
    delete!(s, 1, 2)
    insert!(s, 1, 2, 0, 3)
    @test s.root.data == 1.0
    insert!(s, 1)
    insert!(s, 4)
    @test s.root.right.data == 3.0
    insert!(s, 0.5, -1, -2, 5)
    insert!(s, -3)
    @test s.root.left.left.data == -2.0
    insert!(s, 4.5)
    insert!(s, 6, 5.5, 7, 8)
    @test s.root.right.right.data == 5.5
    delete!(s, 5.5)
    @test s.root.right.right.data == 6.0
    subs =
    delete!(s, 0.5)
    @test s.root.left.data == -1.0
    delete!(s, -1)
    @test s.root.left.right.data == 3.0
    delete!(s, 1)
    @test s.root.left.data == 2.0
    delete!(s, 4, 2)
    @test s.root.left.data == 0.0
    delete!(s, 4.5)
    @test s.root.data == 5
    insert!(s, -1)
    delete!(s, -3, 3, 0, -2, -1)
    @test size(s) == 4
    @test length(s) == 2
    test_show(s)
    test_showmime(s)
end

@testset "SplayTree" begin
    s = Splay(Float64)
    delete!(s, 1)
    delete!(s, 1, 2)
    splay!(s, 1)
    insert!(s, 1, 2, 0, 3, -1)
    @test s.root.data == -1.0
    splay!(s, 3)
    @test s.root.data == 3.0
    @test s.root.left.left.data == -1.0
    insert!(s, 10, 5, -10, -7, -2, 6, -1.5)
    splay!(s, 3)
    @test s.root.right.data == 6.0
    @test s.root.left.data == -1.5
    splay!(s, 1)
    @test s.root.right.data == 3.0
    delete!(s, 0)
    @test !first(search(s, 0.0, true))
    @test s.root.data == -1.0
    delete!(s, -1)
    @test s.root.data == -1.5
    delete!(s, -10)
    delete!(s, -10, 6, 2, 5)
    @test s.root.data == 3.0
    @test s.root.left.data == 1.0
    @test size(s) == 6
    test_show(s)
    test_showmime(s)
end

@testset "SplayTree with topdown operations" begin
    s = Splay(Float64)
    topdowndelete!(s, 1)
    topdowndelete!(s, 1, 2)
    topdownsplay!(s, 1)
    topdowninsert!(s, 1, 2)
    insert!(s, 0, 3, -1)
    topdownsplay!(s, 2)
    @test s.root.data == 2.0
    @test s.root.right.data == 3.0
    topdownsplay!(s, 3)
    @test s.root.left.data == 2.0
    @test s.root.left.left.data == 0.0
    topdowninsert!(s, 1, 0.5)
    @test s.root.left.data == 0.0
    @test s.root.right.data == 1.0
    topdowninsert!(s, 1.5)
    @test s.root.left.data == 1.0
    @test s.root.right.data == 2.0
    topdowndelete!(s, 1.2)
    @test s.root.data == 1.0
    @test s.root.right.data == 1.5
    topdowndelete!(s, 2, 0)
    @test s.root.data == -1
    @test s.root.right.data == 1.0
    topdowninsert!(s, -2.5, -3, -2, -1.5)
    topdownsplay!(s, -1)
    topdownsplay!(s, 0.5)
    @test size(s) == 10
    nodes = [n for n in s.root]
    @test first(nodes) == s.root.left
    test_show(s)
    test_showmime(s)
end


text:
using JuLIP, ACE, JuLIP.Potentials, LinearAlgebra
using JuLIP.Potentials: SZList, ZList
using ACE.Testing: ToyModel, trainset, lsq

function get_basis(species; N = 3, maxdeg = 10, rcut = 7.0 )
   rcut = 7.0
   basis = ACE.Utils.rpi_basis(; species=species, N = N, r0 = 2.5,
   maxdeg = maxdeg, rcut = rcut,
   rin = rnn(:Fe) * 0.6,
   constants = false )
   return basis
end

function calc_rmse(train, Vfit)
   sqe = 0.0
   sqy = 0.0
   for (at, E, F) in train
      sqy += (E/length(at))^2 + norm(mat(F)[:])^2
      Efit = energy(Vfit, at)
      Ffit = forces(Vfit, at)
      sqe += ((E - Efit)/length(at))^2 +
             norm( mat(F)[:] - mat(Ffit)[:] )^2
   end
   return sqrt(sqe / sqy)
end

#---

species = [:Fe, :Al]
at = ACE.Testing.rand_config(ToyModel(species))

basis = get_basis(species; N = 3, maxdeg = 6, rcut = 6.0)
c = rand(length(basis)) .- 0.5
Vfit = JuLIP.MLIPs.combine(basis, c)
@show energy(Vfit, at) ≈ sum(c .* energy(basis, at))
@show forces(Vfit, at) ≈ sum(c .* forces(basis, at))
@show virial(Vfit, at) ≈ sum(c .* virial(basis, at))

#---

@info("Quick run + consistency test")
# species = [:Al, :Fe]
species = [:Fe, :Al]
Vref = ToyModel(species)
train = trainset(Vref, 10)
basis = get_basis(species; N = 3, maxdeg = 6, rcut = 6.0)
c, relrmse = lsq(train, basis)

Vfit = JuLIP.MLIPs.combine(basis, c)

norm(forces(Vfit, at) - sum(c .* forces(basis, at)))

relrmse2 = calc_rmse(train, Vfit)
@assert abs(relrmse - relrmse2) < 0.01

#--- simple convergence test

species = [:Al, :Fe]
Vref = ToyModel(species)
train = trainset(Vref, 300; rattle=0.1)

for maxdeg in [6, 8, 10, 12]
   B = get_basis(species; N = 6, maxdeg = maxdeg, rcut = 7.0)
   lsq(train, B)
end


text:
struct OptimizationResultCacheKey
    model::Symbol
    key::OptimizationContainerKey
end

struct CacheFlushRule
    keep_in_cache::Bool
end

CacheFlushRule() = CacheFlushRule(false)

const DEFAULT_SIMULATION_STORE_CACHE_SIZE_MiB = 1024
const DEFAULT_SIMULATION_STORE_CACHE_SIZE = DEFAULT_SIMULATION_STORE_CACHE_SIZE_MiB * MiB
const MIN_CACHE_FLUSH_SIZE_MiB = 1
const MIN_CACHE_FLUSH_SIZE = MIN_CACHE_FLUSH_SIZE_MiB * MiB

"""
Informs the flusher on what data to keep in cache.
"""
struct CacheFlushRules
    data::Dict{OptimizationResultCacheKey, CacheFlushRule}
    min_flush_size::Int
    max_size::Int
end

function CacheFlushRules(;
    max_size=DEFAULT_SIMULATION_STORE_CACHE_SIZE,
    min_flush_size=MIN_CACHE_FLUSH_SIZE,
)
    return CacheFlushRules(
        Dict{OptimizationResultCacheKey, CacheFlushRule}(),
        min_flush_size,
        max_size,
    )
end

function add_rule!(
    rules::CacheFlushRules,
    model_name,
    op_container_key,
    keep_in_cache::Bool,
)
    key = OptimizationResultCacheKey(model_name, op_container_key)
    rules.data[key] = CacheFlushRule(keep_in_cache)
    return
end

function get_rule(x::CacheFlushRules, model, op_container_key)
    return get_rule(x, OptimizationResultCacheKey(model, op_container_key))
end

get_rule(x::CacheFlushRules, key::OptimizationResultCacheKey) = x.data[key]

mutable struct CacheStats
    hits::Int
    misses::Int
end

CacheStats() = CacheStats(0, 0)

function get_cache_hit_percentage(x::CacheStats)
    total = x.hits + x.misses
    total == 0 && return 0.0
    return x.hits / (total) * 100
end


text:
"""
# SeisTau

SeisTau integrates the [TauPy](https://github.com/anowacki/TauPy.jl) and
[Seis](https://github.com/anowacki/Seis.jl) packages, allowing one to easily
use seismic travel time predictions for 1D Earth models with `Seis.Trace`s,
`Seis.Event`s and `Seis.Station`s.
"""
module SeisTau

import Seis
using Seis: AbstractTrace, GeogEvent, GeogStation, distance_deg, picks
import TauPy

export
    path,
    travel_time,
    travel_times

"""
    add_pick!(t, p::TauPy.Phase, name=p.name) -> (time, name)

Add a travel time pick to the `Trace` `t` from a `TauPy.Phase` arrival.

The pick name will be taken from the phase name by default.
"""
Seis.add_pick!(t::AbstractTrace, p::Union{TauPy.Phase,TauPy.PhaseGeog},
    name=p.name) = add_pick!(t, p.time, name)

"""
    add_picks!(t, phase; model="iasp91", exact=false, sphere=false)

Add travel time picks to the trace `t` for the 1D Earth `model`, for arrivals
with a name matching `phase`.

If `exact` is `true`, only phases which are an exact match for `phase` will
be added.  If `sphere` is `true`, then assume the Earth is spherical.

Available models are: $(TauPy.available_models()).
"""
function Seis.add_picks!(t::AbstractTrace, phase::AbstractString="ttall";
        model="iasp91", exact=false, sphere=false)
    _check_headers_taup(t)
    arrivals = travel_time(t, phase; model=model, sphere=sphere)
    for arr in arrivals
        exact && arr.name != phase && continue
        Seis.add_pick!(t, arr.time, arr.name)
    end
    t.picks
end

"""
    path(t, phase="all"; model="iasp91", exact=false, kwargs...) -> phases::Vector{TauPy.PhaseGeog}

Return ray paths for the geometry specified in the `Trace` t for the 1D Earth `model`,
for arrivals with a name matching `phase`.

If `exact` is `true`, then only return phases whose names are an exact match for `phase`.

    path(event, station, phase="ttall"; kwargs...)

Compute paths for individual `event`s and `station`s instead of taking
these from a `Trace`.

Available models are: $(TauPy.available_models()).

Note that TauPy only computes geographic paths for a spherical Earth.
"""
function path(evt::GeogEvent, sta::GeogStation, phase="ttall";
        model="iasp91", exact=false, kwargs...)
    _check_headers_taup(evt, sta)
    arrivals = TauPy.path(evt.lon, evt.lat, evt.dep, sta.lon, sta.lat, phase; model=model, kwargs...)
    exact && filter!(arr -> arr.name == phase, arrivals)
    arrivals
end

path(t::AbstractTrace, phase="ttall"; kwargs...) = path(t.evt, t.sta, phase; kwargs...)

"""
    travel_time(t, phase="ttall"; model="iasp91", sphere=false) -> phases::Vector{TauPy.Phase}

Return travel times for the geometry specified in the `Trace` `t` for the 1D
Earth `model`, for arrivals with a name matching `phase`.

If `exact` is `true`, then only return phases whose names are an exact match for `phase`.

    travel_time(event, station, phase="all"; kwargs...)

Compute travel times for individual `event`s and `station`s instead of taking
these from a `Trace`.

Available models are: $(TauPy.available_models()).
"""
travel_time(t::AbstractTrace, phase="ttall"; kwargs...) = travel_time(t.evt, t.sta, phase; kwargs...)

function travel_time(evt::GeogEvent, sta::GeogStation, phase::AbstractString="ttall";
        model="iasp91", exact=false, sphere=false, kwargs...)
    _check_headers_taup(evt, sta)
    arrivals = TauPy.travel_time(evt.dep, distance_deg(evt, sta, sphere=sphere), phase; model=model, kwargs...)
    exact && filter!(arr -> arr.name == phase, arrivals)
    arrivals
end

"""
    travel_times(t; model="iasp91", sphere=false) -> ::Vector{Vector{TauPy.Phase}}

Return the list of `TauPy.Phase`s associated with each travel time pick in `t`.
"""
travel_times(t::AbstractTrace; model::AbstractString="iasp91", sphere=false) =
    [travel_time(t, name; model=model, sphere=sphere)
     for (_, name) in picks(t) if name !== missing && occursin(r"^[psPS]", name)]

"Throw an error if a Trace doesn't contain the right headers to call TauPy.travel_time."
_check_headers_taup(evt::GeogEvent, sta::GeogStation) = any(ismissing,
    (evt.lon, evt.lat, evt.dep, sta.lon, sta.lat)) &&
    throw(ArgumentError("Insufficient information in trace to compute travel times." *
                        "  (Need: evt.lon, evt.lat, evt.dep, sta.lon and sta.lat)")) ||
    nothing

_check_headers_taup(t::AbstractTrace) = _check_headers_taup(t.evt, t.sta)

end # module


text:
include("daggen.jl")

using GraphViz
using Dagger
import Dagger: Thunk, Chunk, istask, dependents

# generate the dag that was executed here
const dag = gen_sort_dag(10^6, 40, 4, 40)
const execution_order = [
# fill in the execution order from actual run output
]
const colors = ["black", "red", "green", "blue", "cyan", "orange", "purple", "brown", "yellow"]

function label_attribs(t::Thunk)
    order = 1
    node = 1
    for idx in 1:length(execution_order)
        p = execution_order[idx]
        if p.first == t.id
            order = idx
            node = p.second
            break
        end
    end
    grayscale = 50 + floor(Int, order*49/length(execution_order))
    color = colors[node]

     """ color="$color" style="filled" fillcolor="gray$grayscale" penwidth="5" """
end

function node_label(io, t::Thunk, c)
    if isa(t.f, Function)
        println(io, t.id, " [label=\"", t.f, " - ", t.id, "\"", label_attribs(t), "]")
    else
        println(io, t.id, " [label=\"fn - ", t.id, "\"", label.attribs(t), "]")
    end
    c
end

function node_label(io, t,c)
    l = replace(string(t), "\"", "")
    println(io, dec(hash(t)), " [label=\"$l\"]")
end

global _part_labels = Dict()

function node_label(io, t::Chunk, c)
    _part_labels[t]="part_$c"
    c+1
end

function node_id(t::Thunk)
    t.id
end

function node_id(t)
    dec(hash(t))
end

function node_id(t::Chunk)
    _part_labels[t]
end

function write_dag(io, t)
    !istask(t) && return
    deps = dependents(t)
    c=1
    for k in keys(deps)
        c = node_label(io, k, c)
    end
    for (k, v) in deps
        for dep in v
            if isa(k, Union{Chunk, Thunk})
                println(io, "$(node_id(k)) -> $(node_id(dep))")
            end
        end
    end
end

function viz(t::Thunk)
    io = IOBuffer()
    write_dag(io, t)
    """digraph {
        graph [layout=dot, rankdir=TB];
        $(String(take!(io)))
    }"""
end

graph = viz(dag)
dotfile = joinpath(pwd(), "viz.dot")
htmlfile = joinpath(pwd(), "viz.html")
svgfile = joinpath(pwd(), "viz.svg")
pngfile = joinpath(pwd(), "viz.png")
open(dotfile, "w") do io
    println(io, graph)
end
open(htmlfile, "w") do io
    iob = IOBuffer()
    show(iob, MIME"image/svg+xml"(), Graph(graph))
    svgstr = String(take!(iob))
    svgstart = searchindex(svgstr, "<g id=\"graph0\"")
    svgend = searchindex(svgstr, "</svg>") - 1
    svgstr = svgstr[svgstart:svgend]

    print(io, """

    <input id="horzscale" type="range" min="-100" max="5" style="width: 800px" onchange="onrescale();">

    <input id="vertscale" type="range" min="-10" max="1" style="width: 800px" onchange="onrescale();">

    """)
end
open(svgfile, "w") do io
    show(io, MIME"image/svg+xml"(), Graph(graph))
end
#=
open(pngfile, "w") do io
    show(io, MIME"image/png"(), Graph(graph))
end
=#
#dotbin = joinpath(Pkg.dir("GraphViz"), "deps", "usr", "bin", "dot")
#run(`$dotbin $dotfile -Tsvg -o $pngfile`)


text:
const INTERIOR_COEFFICIENTS = Dict{Int, Vector{Float64}}(
2 => [ 0.482962913145 ; 0.836516303738 ; 0.224143868042 ; -0.129409522551 ]
,
3 => [ .332670552950 ; .806891509311 ; .459877502118 ; -.135011020010 ; -.085441273882 ; .035226291882 ]
,
4 => [ 0.045570345896 ; -0.0178247014417 ; -0.140317624179 ; 0.421234534204 ; 1.13665824341 ; 0.703739068656 ; -0.0419109651251 ; -0.107148901418 ] / sqrt2
,
5 => [ 0.0276321529578 ; -0.0298424998687 ; -0.247951362613 ; 0.0234789231361 ; 0.89658164838 ; 1.02305296689 ; 0.281990696854 ; -0.0553441861166 ; 0.0417468644215 ; 0.0386547959548 ] / sqrt2
,
6 => [ -0.0110318675094 ; 0.00249992209279 ; 0.06325056266 ; -0.0297837512985  ; -0.102724969862 ; 0.477904371333 ; 1.11389278393 ; 0.694457972958 ; -0.0683231215866  ; -0.166863215412 ; 0.00493661237185 ; 0.0217847003266 ] / sqrt2
,
7 => [ 0.014521394762 ; 0.00567134268574 ; -0.152463871896 ; -0.198056706807 ; 0.408183939725 ; 1.08578270981 ; 0.758162601964  ; 0.0246656594886  ; -0.070078291222 ; 0.0960147679355 ; 0.043155452582 ; -0.0178704316511 ; -0.0014812259146 ; 0.0037926585342 ] / sqrt2
,
8 => [ 0.00267279339281 ; -0.000428394300246 ; -0.0211456865284 ; 0.00538638875377 ; 0.0694904659113 ; -0.0384935212634 ; -0.0734625087609 ; 0.515398670374 ; 1.09910663054 ; 0.68074534719 ; -0.0866536154058 ; -0.202648655286 ; 0.0107586117505 ; 0.0448236230437 ; -0.000766690896228 ; -0.0047834585115 ] / sqrt2
)


text:
#=
  Coupon collecter's problem, or card collecter's problem.

  There are N different collecter's cards hidden in a package, but we don't
  know which card there is in the package we buy.
  We want to collect all of them, how many packages must one buy to collect
  all the different cards?

  See https://en.wikipedia.org/wiki/Coupon_collector%27s_problem
  """
  In probability theory, the coupon collector's problem describes 'collect all coupons and win'
  contests. It asks the following question: If each box of a brand of cereals contains a
  coupon, and there are n different types of coupons, what is the probability that more
  than t boxes need to be bought to collect all n coupons?

  An alternative statement is: Given n coupons, how many coupons do you expect you need
  to draw with replacement before having drawn each coupon at least once? The mathematical
  analysis of the problem reveals that the expected number of trials needed grows as
  Θ(n log(n).
  For example, when n = 50 it takes about 225[b] trials on average to collect all 50 coupons.

  ...

  [b]: E(50) = 50(1 + 1/2 + 1/3 + ... + 1/50) = 224.9603, the expected number of trials to
  collect all 50 coupons.

   The approximation n*log(n) + γ*n + 1/2 for this expected number gives in this case
   ≈ 195.6011 + 28.8608 + 0.5 ≈ 224.9619.  [log is the natural logarithm]
  """

  Note: The assumptions are that the slots are filled using a Poisson process
  and the parameter will influence the result much.

  p (mean): 5.7496

  Distributions of variable p (num:0)
  7.00000 =>    3172  (0.317200)
  6.00000 =>    2920  (0.292000)
  5.00000 =>    2349  (0.234900)
  4.00000 =>    1317  (0.131700)
  3.00000 =>     242  (0.024200)

  Distributions of variable total_filled (num:0)
  50.00000 =>   10000  (1.000000)

  Distributions of variable theoretical (num:0)
  224.96027 =>   10000  (1.000000)

  mean sumx: 288.7327

  See ~/blog/coupon_collector.blog
      ~/webppl/cookie_problem.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function coupon_collector()
    n = 50
    m = round(Int,sqrt(n))
    #  "Exact" probability from https:# en.wikipedia.org/wiki/Coupon_collector%27s_problem (footnote [b])
    theoretical ~ Dirac(n*sum([1/i for i in 1:n])) #  N*Math.log(N,Math.e);

    #  p = beta(14,14);
    p ~ DiscreteUniform(1,m)

    #  Fill this slot with cards
    fill ~ filldist(Poisson(p),n)

    #  How many slots are filled (i.e. > 0)?
    total_filled ~ Dirac(sum([fill[c] > 0 for c in 1:n]))

    true ~ Dirac(total_filled == n)

    sumx ~ Dirac(sum(fill))
    sum0s ~ Dirac(sum([fill[c] == 0 for c in 1:n] ))

end

model = coupon_collector()
num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns, :p)
show_var_dist_pct(chns, :total_filled)
show_var_dist_pct(chns, :theoretical)
println(["mean sumx: ", mean(chns[:sumx].data)])

show_var_dist_pct(chns, Symbol("fill[1]"))
show_var_dist_pct(chns, Symbol("fill[2]"))
show_var_dist_pct(chns, Symbol("fill[3]"))
show_var_dist_pct(chns, Symbol("fill[4]"))
show_var_dist_pct(chns, Symbol("fill[5]"))
show_var_dist_pct(chns, Symbol("fill[6]"))
show_var_dist_pct(chns, Symbol("fill[7]"))
show_var_dist_pct(chns, Symbol("fill[8]"))
show_var_dist_pct(chns, Symbol("fill[9]"))
show_var_dist_pct(chns, Symbol("fill[10]"))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# Test of Depth-first visit

using Graphs, SparseArrays
using Test

mutable struct GraphTest
    graph_edges::Array{Tuple{Int,Int},1}
    dfs_path::Array{Int,1}
    is_cyclic::Bool
    topo_sort::Array{Int,1}

    GraphTest(gedges, dfspath, iscyclic, toposort) =
        new(gedges, dfspath, iscyclic, toposort)
end

global dir_acyclic = GraphTest(
    [(1,2), (1,3), (1,6), (2,4), (2,5), (3,5), (3,6)],
    [1, 2, 4, 5, 3, 6],
    false,
    [1,3,6,2,5,4])
global undir_acyclic = GraphTest(
    [(1,2), (1,3), (1,6), (2,4), (2,5)],
    [],
    false,
    [1,6,3,2,5,4])
global cyclic  = GraphTest(
    [(1,2), (1,3), (1,6), (2,4), (2,5), (3,5), (3,6), (5,1)],
    [1, 2, 4, 5, 3, 6],
    true,
    [])

global testsets = [
    (true, [dir_acyclic, cyclic]),
    (false, [undir_acyclic, cyclic])]

for tset in testsets
    global (is_dir, graphtests) = tset

    for gtest in graphtests

        global g = simple_inclist(6, is_directed = is_dir)
        map((edg) -> add_edge!(g, edg[1], edg[2]), gtest.graph_edges)

        global gEx = graph(ExVertex[], ExEdge{ExVertex}[], is_directed = is_dir)
        map((x) -> add_vertex!(gEx, "edge:" * string(x)), 1:6)
        global VV = vertices(gEx)
        map((edg) -> add_edge!(gEx, VV[edg[1]], VV[edg[2]]), gtest.graph_edges)

        # DFS traversal
        if !isempty(gtest.dfs_path)
            global vs1 = visited_vertices(g, DepthFirst(), 1)
            @assert vs1 == gtest.dfs_path

            global vs2 = visited_vertices(gEx, DepthFirst(), VV[1])
            @assert vs2 == collect(map((x) -> gEx.vertices[x], gtest.dfs_path))
        end

        # Cyclic test
        @assert test_cyclic_by_dfs(g) == gtest.is_cyclic
        @assert test_cyclic_by_dfs(gEx) == gtest.is_cyclic

        # Topological sort
        if gtest.is_cyclic
            @test_throws ArgumentError topological_sort_by_dfs(g)  # g2 contains a loop
            @test_throws ArgumentError topological_sort_by_dfs(gEx)  # g2 contains a loop

        elseif !isempty(gtest.topo_sort)
            ts = topological_sort_by_dfs(g)
            @assert ts == gtest.topo_sort

            ts = topological_sort_by_dfs(gEx)
            @assert [vertex_index(e,gEx) for e in ts] == gtest.topo_sort
        end
    end
end


text:
module Reflect

  export anonymous_type, anonymous_module, define_function

  # Stores the ID of the next anonymous type.
  global __counter = zero(Int)
  global __counter_mod = zero(Int)

  # Generates and returns the name of the next anonymous type.
  fresh() = "Anon$(global __counter += 1)"
  fresh_mod() = "AnonMod$(global __counter_mod += 1)"

  # Generates an empty anonymous type.
  anonymous_type(m::Module) = anonymous_type(m, "type;end")

  # Generates an anonymous module from a provided definition.
  function anonymous_module(d::AbstractString)
    name = fresh_mod()
    d = "module $(name)\n$(d)\nend"
    eval(Base.parse(d))
    return eval(Base.parse(name))
  end

  # Generates an anonymous type from a provided definition.
  function anonymous_type(m::Module, d::AbstractString)
    # Insert the name of the type after the first modifier.
    name = fresh()
    d = replace(d, r"(?<=\w)\b", " $(name)", 1)

    # Format the constructors using the name of the type.
    d = replace(d, r"constructor(?=\()", name)

    eval(m, Base.parse(d))
    return eval(m, Base.parse(name))
  end
end


text:
module TestDocs

"""
This is a test function
    """
foo(x) = 1

end # module