text:
using Test, Tables, TableTraits, DataValues, QueryOperators, IteratorInterfaceExtensions

@testset "utils.jl" begin

    NT = NamedTuple{(), Tuple{}}
    @test Tables.names(NT) === ()
    @test Tables.types(NT) === Tuple{}
    @test isempty(Tables.runlength(Tables.types(NT)))
    @test Tables.columnindex(Tables.names(NT), :i) == 0
    @test Tables.columntype(Tables.names(NT), Tables.types(NT), :i) == Union{}

    NT = NamedTuple{(:a, :b, :c), NTuple{3, Int64}}
    @test Tables.names(NT) === (:a, :b, :c)
    @test Tables.types(NT) === Tuple{Int64, Int64, Int64}
    @test Tables.runlength(Tables.types(NT)) == [(Int64, 3)]
    @test Tables.columnindex(Tables.names(NT), :a) == 1
    @test Tables.columnindex(Tables.names(NT), :i) == 0
    @test Tables.columntype(Tables.names(NT), Tables.types(NT), :a) == Int64
    @test Tables.columntype(Tables.names(NT), Tables.types(NT), :i) == Union{}

    NT = NamedTuple{Tuple(Symbol("a$i") for i = 1:20), Tuple{vcat(fill(Int, 10), fill(String, 10))...}}
    @test Tables.names(NT) === Tuple(Symbol("a$i") for i = 1:20)
    @test Tables.types(NT) === Tuple{vcat(fill(Int, 10), fill(String, 10))...}
    @test Tables.runlength(Tables.types(NT)) == [(Int, 10), (String, 10)]
    @test Tables.columnindex(Tables.names(NT), :a20) == 20
    @test Tables.columnindex(Tables.names(NT), :i) == 0
    @test Tables.columntype(Tables.names(NT), Tables.types(NT), :a20) == String
    @test Tables.columntype(Tables.names(NT), Tables.types(NT), :i) == Union{}

    nt = (a=1, b=2, c=3)
    @test getproperty(nt, Int, 1, :a) === 1

    NT = typeof(nt)
    output = [0, 0, 0]
    Tables.eachcolumn(Tables.Schema(Tables.names(NT), Tables.types(NT)), nt, output) do val, col, nm, out
        out[col] = val
    end
    @test output == [1, 2, 3]

    nt = NamedTuple{Tuple(Symbol("a$i") for i = 1:101)}(Tuple(i for i = 1:101))
    NT = typeof(nt)
    @test Tables.runlength(Tables.types(NT)) == [(Int, 101)]
    output = zeros(Int, 101)
    Tables.eachcolumn(Tables.Schema(Tables.names(NT), Tables.types(NT)), nt, output) do val, col, nm, out
        out[col] = val
    end
    @test output == [i for i = 1:101]

    nt = NamedTuple{Tuple(Symbol("a$i") for i = 1:101)}(Tuple(i % 2 == 0 ? i : "$i" for i = 1:101))
    NT = typeof(nt)
    @test Tables.runlength(Tables.types(NT)) == [i % 2 == 0 ? (Int, 1) : (String, 1) for i = 1:101]
    output = Vector{Any}(undef, 101)
    Tables.eachcolumn(Tables.Schema(Tables.names(NT), Tables.types(NT)), nt, output) do val, col, nm, out
        out[col] = val
    end
    @test output == [i % 2 == 0 ? i : "$i" for i = 1:101]

    nt = (a=Ref(0), b=Ref(0))
    Tables.eachcolumn(Tables.Schema((:a, :b), nothing), (a=1, b=2)) do val, col, nm
        nt[nm][] = val
    end
    @test nt.a[] == 1
    @test nt.b[] == 2

    nt = (a=[1,2,3], b=[4,5,6])
    @test collect(Tables.eachcolumn(nt)) == [[1,2,3], [4,5,6]]
    @test Tables.columnindex(nt, :i) == 0
    @test Tables.columnindex(nt, :a) == 1
    @test Tables.columntype(nt, :a) == Int
    @test Tables.columntype(nt, :i) == Union{}

    rows = Tables.rows(nt)
    @test eltype(rows) == Tables.ColumnsRow{typeof(nt)}
    @test Tables.schema(rows) == Tables.Schema((:a, :b), (Int, Int))
    row = first(rows)
    @test row.a == 1

    @test Tables.sym(1) === 1
    @test Tables.sym("hey") == :hey

    @test propertynames(Tables.Schema((:a, :b), nothing)) == (:names, :types)

    v = Tables.EmptyVector(1)
    @test_throws UndefRefError v[1]
    @test Base.IndexStyle(typeof(v)) == Base.IndexLinear()

    @test Tables.istable(Tables.CopiedColumns)
    @test Tables.columnaccess(Tables.CopiedColumns)
    c = Tables.CopiedColumns(nt)
    @test Tables.columns(c) === c
    @test Tables.materializer(c) == Tables.materializer(nt)

    @test_throws ArgumentError Tables.columntable([1,2,3])

    tt = [(1,2,3), (4,5,6)]
    r = Tables.nondatavaluerows(tt)
    row = first(r)
    @test getproperty(row, 1) == 1
    @test Tables.columntable(tt) == NamedTuple{(Symbol("1"), Symbol("2"), Symbol("3"))}(([1, 4], [2, 5], [3, 6]))

    @test Tables.getarray([1,2,3]) == [1,2,3]
    @test Tables.getarray((1,2,3)) == [1,2,3]
end

@testset "namedtuples.jl" begin

    nt = (a=1, b=2, c=3)
    rt = [nt, nt, nt]
    @test Tables.rows(rt) === rt
    @test Tables.schema(rt).names == Tables.names(typeof(nt))
    @test Tables.namedtupleiterator(eltype(rt), rt) === rt

    rt = [(a=1, b=4.0, c="7"), (a=2, b=5.0, c="8"), (a=3, b=6.0, c="9")]
    nt = (a=[1,2,3], b=[4.0, 5.0, 6.0], c=["7", "8", "9"])
    @test Tables.rowcount(nt) == 3
    @test Tables.schema(nt) == Tables.Schema((:a, :b, :c), Tuple{Int, Float64, String})
    @test Tables.istable(typeof(nt))
    @test Tables.columnaccess(typeof(nt))
    @test Tables.columns(nt) === nt
    @test rowtable(nt) == rt
    @test columntable(rt) == nt
    @test rt == (rt |> columntable |> rowtable)
    @test nt == (nt |> rowtable |> columntable)

    @test Tables.buildcolumns(nothing, rt) == nt
    @test Tables.columntable(nothing, nt) == nt

    # test push!
    rtf = Iterators.Filter(x->x.a >= 1, rt)
    @test Tables.columntable(rtf) == nt
    @test Tables.buildcolumns(nothing, rtf) == nt

    # append
    nt2 = columntable(nt, rt)
    @test Tables.rowcount(nt2) == 6
    @test Tables.schema(nt2) == Tables.Schema((:a, :b, :c), Tuple{Int, Float64, String})
    @test nt2 == (a = [1, 2, 3, 1, 2, 3], b = [4.0, 5.0, 6.0, 4.0, 5.0, 6.0], c = ["7", "8", "9", "7", "8", "9"])
    rt2 = rowtable(rt, nt)
    @test length(rt2) == 9

    rt = [(a=1, b=4.0, c="7"), (a=2.0, b=missing, c="8"), (a=3, b=6.0, c="9")]
    @test Tables.istable(typeof(rt))
    @test Tables.rowaccess(typeof(rt))
    tt = Tables.buildcolumns(nothing, rt)
    @test isequal(tt, (a = [1.0, 2.0, 3.0], b = Union{Missing, Float64}[4.0, missing, 6.0], c = ["7", "8", "9"]))
    @test tt.a[1] === 1.0
    @test tt.a[2] === 2.0
    @test tt.a[3] === 3.0

    nti = Tables.NamedTupleIterator{Nothing, typeof(rt)}(rt)
    @test Base.IteratorEltype(typeof(nti)) == Base.EltypeUnknown()
    @test Base.IteratorSize(typeof(nti)) == Base.HasShape{1}()
    @test length(nti) == 3
    nti2 = collect(nti)
    @test isequal(rt, nti2)
    nti = Tables.NamedTupleIterator{typeof(Tables.Schema((:a, :b, :c), (Union{Int, Float64}, Union{Float64, Missing}, String))), typeof(rt)}(rt)
    @test eltype(typeof(nti)) == NamedTuple{(:a, :b, :c),Tuple{Union{Float64, Int},Union{Missing, Float64},String}}

    # test really wide tables
    nms = Tuple(Symbol("i", i) for i = 1:101)
    vals = Tuple(rand(Int, 3) for i = 1:101)
    nt = NamedTuple{nms}(vals)
    rt = Tables.rowtable(nt)
    @test length(rt) == 3
    @test length(rt[1]) == 101
    @test eltype(rt).parameters[1] == nms
    @test Tables.columntable(rt) == nt
    @test Tables.buildcolumns(nothing, rt) == nt
end

@testset "Materializer" begin
    rt = [(a=1, b=4.0, c="7"), (a=2, b=5.0, c="8"), (a=3, b=6.0, c="9")]
    nt = (a=[1,2,3], b=[4.0, 5.0, 6.0], c=["7", "8", "9"])

    @test nt == Tables.materializer(nt)(Tables.columns(nt))
    @test nt == Tables.materializer(nt)(Tables.columns(rt))
    @test nt == Tables.materializer(nt)(rt)
    @test rt == Tables.materializer(rt)(nt)

    function select(table, cols::Symbol...)
        Tables.istable(table) || throw(ArgumentError("select requires a table input"))
        nt = Tables.columntable(table)  # columntable(t) creates a NamedTuple of AbstractVectors
        newcols = NamedTuple{cols}(nt)
        Tables.materializer(table)(newcols)
    end

    @test select(nt, :a, :b, :c) == nt
    @test select(nt, :c, :a) == NamedTuple{(:c, :a)}(nt)
    @test select(rt, :a) == [(a=1,), (a=2,), (a=3,)]

    @test Tables.materializer(1) === Tables.columntable
end

@testset "Matrix integration" begin
    rt = [(a=1, b=4.0, c="7"), (a=2, b=5.0, c="8"), (a=3, b=6.0, c="9")]
    nt = (a=[1,2,3], b=[4.0, 5.0, 6.0])

    mat = Tables.matrix(rt)
    @test nt.a == mat[:, 1]
    @test size(mat) == (3, 3)
    @test eltype(mat) == Any
    @test_throws ArgumentError Tables.rows(mat)
    @test_throws ArgumentError Tables.columns(mat)
    mat2 = Tables.matrix(nt)
    @test eltype(mat2) == Float64
    @test mat2[:, 1] == nt.a
    @test !Tables.istable(mat2)
    @test !Tables.istable(typeof(mat2))
    mat3 = Tables.matrix(nt; transpose=true)
    @test size(mat3) == (2, 3)
    @test mat3[1, :] == nt.a
    @test mat3[2, :] == nt.b

    tbl = Tables.table(mat) |> columntable
    @test keys(tbl) == (:Column1, :Column2, :Column3)
    @test tbl.Column1 == [1, 2, 3]
    tbl2 = Tables.table(mat2) |> rowtable
    @test length(tbl2) == 3
    @test map(x->x.Column1, tbl2) == [1.0, 2.0, 3.0]

    mattbl = Tables.table(mat)
    @test Tables.istable(typeof(mattbl))
    @test Tables.rowaccess(typeof(mattbl))
    @test Tables.rows(mattbl) === mattbl
    @test Tables.columnaccess(typeof(mattbl))
    @test Tables.columns(mattbl) === mattbl
    @test mattbl.Column1 == [1,2,3]
    matrow = first(mattbl)
    @test eltype(mattbl) == typeof(matrow)
    @test matrow.Column1 == 1
    @test propertynames(mattbl) == propertynames(matrow) == [:Column1, :Column2, :Column3]
end

import Base: ==
struct GenericRow
    a::Int
    b::Float64
    c::String
end
==(a::GenericRow, b::GenericRow) = a.a == b.a && a.b == b.b && a.c == b.c

struct GenericRowTable
    data::Vector{GenericRow}
end
==(a::GenericRowTable, b::GenericRowTable) = all(a.data .== b.data)

Base.eltype(g::GenericRowTable) = GenericRow
Base.length(g::GenericRowTable) = length(g.data)
Base.size(g::GenericRowTable) = (length(g.data),)
Tables.istable(::Type{GenericRowTable}) = true
Tables.rowaccess(::Type{GenericRowTable}) = true
Tables.rows(x::GenericRowTable) = x
Tables.schema(x::GenericRowTable) = Tables.Schema((:a, :b, :c), Tuple{Int, Float64, String})

function Base.iterate(g::GenericRowTable, st=1)
    st > length(g.data) && return nothing
    return g.data[st], st + 1
end

genericrowtable(x) = GenericRowTable(collect(map(x->GenericRow(x.a, x.b, x.c), Tables.rows(x))))

struct GenericColumn{T} <: AbstractVector{T}
    data::Vector{T}
end
Base.eltype(g::GenericColumn{T}) where {T} = T
Base.length(g::GenericColumn) = length(g.data)
==(a::GenericColumn, b::GenericColumn) = a.data == b.data
Base.getindex(g::GenericColumn, i::Int) = g.data[i]

struct GenericColumnTable
    names::Dict{Symbol, Int}
    data::Vector{GenericColumn}
end

Tables.istable(::Type{GenericColumnTable}) = true
Tables.columnaccess(::Type{GenericColumnTable}) = true
Tables.columns(x::GenericColumnTable) = x
Tables.schema(g::GenericColumnTable) = Tables.Schema(Tuple(keys(getfield(g, 1))), Tuple{(eltype(x) for x in getfield(g, 2))...})
Base.getproperty(g::GenericColumnTable, nm::Symbol) = getfield(g, 2)[getfield(g, 1)[nm]]
Base.propertynames(g::GenericColumnTable) = Tuple(keys(getfield(g, 1)))

function genericcolumntable(x)
    cols = Tables.columns(x)
    sch = Tables.schema(x)
    data = [GenericColumn(getproperty(cols, nm)) for nm in sch.names]
    return GenericColumnTable(Dict(nm=>i for (i, nm) in enumerate(sch.names)), data)
end
==(a::GenericColumnTable, b::GenericColumnTable) = getfield(a, 1) == getfield(b, 1) && getfield(a, 2) == getfield(b, 2)

@testset "Tables.jl interface" begin

    @test !Tables.istable(1)
    @test !Tables.istable(Int)
    @test !Tables.rowaccess(1)
    @test !Tables.rowaccess(Int)
    @test !Tables.columnaccess(1)
    @test !Tables.columnaccess(Int)
    @test Tables.schema(1) === nothing

    sch = Tables.Schema{(:a, :b), Tuple{Int64, Float64}}()
    @test Tables.Schema((:a, :b), Tuple{Int64, Float64}) === sch
    @test Tables.Schema(NamedTuple{(:a, :b), Tuple{Int64, Float64}}) === sch
    @test Tables.Schema((:a, :b), nothing) === Tables.Schema{(:a, :b), nothing}()
    @test Tables.Schema([:a, :b], [Int64, Float64]) === sch
    show(sch)
    @test sch.names == (:a, :b)
    @test sch.types == (Int64, Float64)
    @test_throws ArgumentError sch.foobar

    gr = GenericRowTable([GenericRow(1, 4.0, "7"), GenericRow(2, 5.0, "8"), GenericRow(3, 6.0, "9")])
    gc = GenericColumnTable(Dict(:a=>1, :b=>2, :c=>3), [GenericColumn([1,2,3]), GenericColumn([4.0, 5.0, 6.0]), GenericColumn(["7", "8", "9"])])
    @test gc == (gr |> genericcolumntable)
    @test gr == (gc |> genericrowtable)
    @test gr == (gr |> genericrowtable)

    @test_throws ArgumentError Tables.columns(Int64)
    @test_throws ArgumentError Tables.rows(Int64)
end

@testset "isless" begin
    t = (x = [1, 1, 0, 2], y = [-1, 1, 3, 2])
    a,b,c,d = Tables.rows(t)
    @test isless(a, b)
    @test isless(c, d)
    @test !isless(d, a)
    @test !isequal(a, b)
    @test isequal(a, a)
    @test sortperm([a, b, c, d]) == [3, 1, 2, 4]
end

struct ColumnSource
end

TableTraits.supports_get_columns_copy_using_missing(::ColumnSource) = true

function TableTraits.get_columns_copy_using_missing(x::ColumnSource)
    return (a=[1,2,3], b=[4.,5.,6.], c=["A", "B", "C"])
end

let x=ColumnSource()
    @test Tables.source(Tables.columns(x)) == Tables.source(Tables.CopiedColumns(TableTraits.get_columns_copy_using_missing(x)))
end

struct ColumnSource2
end

IteratorInterfaceExtensions.isiterable(x::ColumnSource2) = true
TableTraits.isiterabletable(::ColumnSource2) = true

IteratorInterfaceExtensions.getiterator(::ColumnSource2) =
    Tables.rows((a=[1,2,3], b=[4.,5.,6.], c=["A", "B", "C"]))

let x=ColumnSource2()
    @test Tables.source(Tables.columns(x)) == (a=[1,2,3], b=[4.,5.,6.], c=["A", "B", "C"])
end

@testset "operations.jl" begin
ctable = (A=[1, missing, 3], B=[1.0, 2.0, 3.0], C=["hey", "there", "sailor"])
rtable = Tables.rowtable(ctable)

## Tables.transform
tran = ctable |> Tables.transform(C=Symbol)
@test Tables.istable(typeof(tran))
@test !Tables.rowaccess(typeof(tran))
@test Tables.columnaccess(typeof(tran))
@test Tables.columns(tran) === tran
@test IteratorInterfaceExtensions.isiterable(tran)
@test typeof(IteratorInterfaceExtensions.getiterator(tran)) <: Tables.DataValueRowIterator

tran2 = rtable |> Tables.transform(C=Symbol)
@test Tables.istable(typeof(tran2))
@test Tables.rowaccess(typeof(tran2))
@test !Tables.columnaccess(typeof(tran2))
@test Tables.rows(tran2) === tran2
@test Base.IteratorSize(typeof(tran2)) == Base.HasShape{1}()
@test length(tran2) == 3
@test eltype(tran2) == Tables.TransformsRow{NamedTuple{(:A, :B, :C),Tuple{Union{Missing, Int},Float64,String}},NamedTuple{(:C,),Tuple{DataType}}}
trow = first(tran2)
@test trow.A === 1
@test trow.B === 1.0
@test trow.C == :hey
ctable2 = Tables.columntable(tran2)
@test isequal(ctable2.A, ctable.A)
@test ctable2.C == map(Symbol, ctable.C)

# test various ways of inputting Tables.transform functions
table = Tables.transform(ctable, Dict{String, Base.Callable}("C" => Symbol)) |> Tables.columntable
@test table.C == [:hey, :there, :sailor]

table = ctable |> Tables.transform(C=Symbol) |> Tables.columntable
@test table.C == [:hey, :there, :sailor]

table = Tables.transform(ctable, Dict{Symbol, Base.Callable}(:C => Symbol)) |> Tables.columntable
@test table.C == [:hey, :there, :sailor]

table = Tables.transform(ctable, Dict{Int, Base.Callable}(3 => Symbol)) |> Tables.columntable
@test table.C == [:hey, :there, :sailor]

# test simple Tables.transforms + return types
table = ctable |> Tables.transform(Dict("A"=>x->x+1)) |> Tables.columntable
@test isequal(table.A, [2, missing, 4])
@test typeof(table.A) == Vector{Union{Missing, Int}}

table = ctable |> Tables.transform(Dict("A"=>x->coalesce(x+1, 0))) |> Tables.columntable
@test table.A == [2, 0, 4]

table = ctable |> Tables.transform(Dict("A"=>x->coalesce(x+1, 0.0))) |> Tables.columntable
@test table.A == [2, 0.0, 4]

table = ctable |> Tables.transform(Dict(2=>x->x==2.0 ? missing : x)) |> Tables.columntable
@test isequal(table.B, [1.0, missing, 3.0])
@test typeof(table.B) == Vector{Union{Float64, Missing}}

# test row sinks
# test various ways of inputting Tables.transform functions
table = Tables.transform(ctable, Dict{String, Base.Callable}("C" => Symbol)) |> Tables.rowtable
@test table[1].C == :hey

table = ctable |> Tables.transform(C=Symbol) |> Tables.rowtable
@test table[1].C == :hey

table = Tables.transform(ctable, Dict{Symbol, Base.Callable}(:C => Symbol)) |> Tables.rowtable
@test table[1].C == :hey

table = Tables.transform(ctable, Dict{Int, Base.Callable}(3 => Symbol)) |> Tables.rowtable
@test table[1].C == :hey

# test simple transforms + return types
table = ctable |> Tables.transform(Dict("A"=>x->x+1)) |> Tables.rowtable
@test isequal(map(x->x.A, table), [2, missing, 4])
@test typeof(map(x->x.A, table)) == Vector{Union{Missing, Int}}

table = ctable |> Tables.transform(Dict("A"=>x->coalesce(x+1, 0))) |> Tables.rowtable
@test map(x->x.A, table) == [2, 0, 4]

table = ctable |> Tables.transform(Dict("A"=>x->coalesce(x+1, 0.0))) |> Tables.rowtable
@test map(x->x.A, table) == [2, 0.0, 4]

table = ctable |> Tables.transform(Dict(2=>x->x==2.0 ? missing : x)) |> Tables.rowtable
@test isequal(map(x->x.B, table), [1.0, missing, 3.0])
@test typeof(map(x->x.B, table)) == Vector{Union{Float64, Missing}}

## Tables.select

# 117
sel = Tables.select(ctable)
@test Tables.istable(typeof(sel))
@test Tables.schema(sel) == Tables.Schema((), ())
@test Tables.columnaccess(typeof(sel))
@test Tables.columns(sel) === sel
@test propertynames(sel) == ()
@test Tables.columntable(sel) == NamedTuple()
@test Tables.rowtable(sel) == NamedTuple{(), Tuple{}}[]

sel = ctable |> Tables.select(:A)
@test Tables.istable(typeof(sel))
@test IteratorInterfaceExtensions.isiterable(sel)
@test IteratorInterfaceExtensions.getiterator(sel) == Tables.datavaluerows(sel)
@test Tables.schema(sel) == Tables.Schema((:A,), (Union{Int, Missing},))
@test Tables.columnaccess(typeof(sel))
@test Tables.columns(sel) === sel
@test propertynames(sel) == (:A,)

sel = ctable |> Tables.select(1)
@test Tables.istable(typeof(sel))
@test IteratorInterfaceExtensions.isiterable(sel)
@test IteratorInterfaceExtensions.getiterator(sel) == Tables.datavaluerows(sel)
@test Tables.schema(sel) == Tables.Schema((:A,), (Union{Int, Missing},))
@test Tables.columnaccess(typeof(sel))
@test Tables.columns(sel) === sel
@test propertynames(sel) == (:A,)

sel = Tables.select(rtable)
@test Tables.rowaccess(typeof(sel))
@test Tables.rows(sel) === sel
@test Tables.schema(sel) == Tables.Schema((), ())
@test Base.IteratorSize(typeof(sel)) == Base.HasShape{1}()
@test length(sel) == 3
@test Base.IteratorEltype(typeof(sel)) == Base.HasEltype()
@test eltype(sel) == Tables.SelectRow{NamedTuple{(:A, :B, :C),Tuple{Union{Missing, Int},Float64,String}},()}
@test Tables.columntable(sel) == NamedTuple()
@test Tables.rowtable(sel) == [NamedTuple(), NamedTuple(), NamedTuple()]
srow = first(sel)
@test propertynames(srow) == ()

sel = rtable |> Tables.select(:A)
@test Tables.rowaccess(typeof(sel))
@test Tables.rows(sel) === sel
@test Tables.schema(sel) == Tables.Schema((:A,), (Union{Int, Missing},))
@test Base.IteratorSize(typeof(sel)) == Base.HasShape{1}()
@test length(sel) == 3
@test Base.IteratorEltype(typeof(sel)) == Base.HasEltype()
@test eltype(sel) == Tables.SelectRow{NamedTuple{(:A, :B, :C),Tuple{Union{Missing, Int},Float64,String}},(:A,)}
@test isequal(Tables.columntable(sel), (A = [1, missing, 3],))
@test isequal(Tables.rowtable(sel), [(A=1,), (A=missing,), (A=3,)])
srow = first(sel)
@test propertynames(srow) == (:A,)

sel = rtable |> Tables.select(1)
@test Tables.rowaccess(typeof(sel))
@test Tables.rows(sel) === sel
@test Tables.schema(sel) == Tables.Schema((:A,), (Union{Int, Missing},))
@test Base.IteratorSize(typeof(sel)) == Base.HasShape{1}()
@test length(sel) == 3
@test Base.IteratorEltype(typeof(sel)) == Base.HasEltype()
@test eltype(sel) == Tables.SelectRow{NamedTuple{(:A, :B, :C),Tuple{Union{Missing, Int},Float64,String}},(1,)}
@test isequal(Tables.columntable(sel), (A = [1, missing, 3],))
@test isequal(Tables.rowtable(sel), [(A=1,), (A=missing,), (A=3,)])
srow = first(sel)
@test propertynames(srow) == (:A,)

table = ctable |> Tables.select(:A) |> Tables.columntable
@test length(table) == 1
@test isequal(table.A, [1, missing, 3])

table = ctable |> Tables.select(1) |> Tables.columntable
@test length(table) == 1
@test isequal(table.A, [1, missing, 3])

table = ctable |> Tables.select("A") |> Tables.columntable
@test length(table) == 1
@test isequal(table.A, [1, missing, 3])

# column re-ordering
table = ctable |> Tables.select(:A, :C) |> Tables.columntable
@test length(table) == 2
@test isequal(table.A, [1, missing, 3])
@test isequal(table[2], ["hey", "there", "sailor"])

table = ctable |> Tables.select(1, 3) |> Tables.columntable
@test length(table) == 2
@test isequal(table.A, [1, missing, 3])
@test isequal(table[2], ["hey", "there", "sailor"])

table = ctable |> Tables.select(:C, :A) |> Tables.columntable
@test isequal(ctable.A, table.A)
@test isequal(ctable[1], table[2])

table = ctable |> Tables.select(3, 1) |> Tables.columntable
@test isequal(ctable.A, table.A)
@test isequal(ctable[1], table[2])

# row sink
table = ctable |> Tables.select(:A) |> Tables.rowtable
@test length(table[1]) == 1
@test isequal(map(x->x.A, table), [1, missing, 3])

table = ctable |> Tables.select(1) |> Tables.rowtable
@test length(table[1]) == 1
@test isequal(map(x->x.A, table), [1, missing, 3])

table = ctable |> Tables.select("A") |> Tables.rowtable
@test length(table[1]) == 1
@test isequal(map(x->x.A, table), [1, missing, 3])

# column re-ordering
table = ctable |> Tables.select(:A, :C) |> Tables.rowtable
@test length(table[1]) == 2
@test isequal(map(x->x.A, table), [1, missing, 3])
@test isequal(map(x->x[2], table), ["hey", "there", "sailor"])

table = ctable |> Tables.select(1, 3) |> Tables.rowtable
@test length(table[1]) == 2
@test isequal(map(x->x.A, table), [1, missing, 3])
@test isequal(map(x->x[2], table), ["hey", "there", "sailor"])

table = ctable |> Tables.select(:C, :A) |> Tables.rowtable
@test isequal(ctable.A, map(x->x.A, table))
@test isequal(ctable[1], map(x->x[2], table))

table = ctable |> Tables.select(3, 1) |> Tables.rowtable
@test isequal(ctable.A, map(x->x.A, table))
@test isequal(ctable[1], map(x->x[2], table))

end

@testset "TableTraits integration" begin
    rt = (a = Real[1, 2.0, 3], b = Union{Missing, Float64}[4.0, missing, 6.0], c = ["7", "8", "9"])

    dv = Tables.datavaluerows(rt)
    @test Base.IteratorSize(typeof(dv)) == Base.HasLength()
    @test eltype(dv) == NamedTuple{(:a, :b, :c),Tuple{Real,DataValue{Float64},String}}
    @test_throws MethodError size(dv)
    rt2 = collect(dv)
    @test rt2[1] == (a = 1, b = DataValue{Float64}(4.0), c = "7")

    ei = Tables.nondatavaluerows(QueryOperators.EnumerableIterable{eltype(dv), typeof(dv)}(dv))
    @test Tables.istable(typeof(ei))
    @test Tables.rowaccess(typeof(ei))
    @test Tables.rows(ei) === ei
    @test Base.IteratorEltype(typeof(ei)) == Base.HasEltype()
    @test Base.IteratorSize(typeof(ei)) == Base.HasLength()
    @test eltype(ei) == Tables.IteratorRow{NamedTuple{(:a, :b, :c),Tuple{Real,DataValue{Float64},String}}}
    @test eltype(typeof(ei)) == Tables.IteratorRow{NamedTuple{(:a, :b, :c),Tuple{Real,DataValue{Float64},String}}}
    @test_throws MethodError size(ei)
    nt = ei |> columntable
    @test isequal(rt, nt)
    rt3 = ei |> rowtable
    @test isequal(rt |> rowtable, rt3)

    # rt = [(a=1, b=4.0, c="7"), (a=2, b=5.0, c="8"), (a=3, b=6.0, c="9")]
    mt = Tables.nondatavaluerows(ei.x |> y->QueryOperators.map(y, x->(a=x.b, c=x.c), Expr(:block)))
    @inferred (mt |> columntable)
    @inferred (mt |> rowtable)

    # uninferrable case
    mt = Tables.nondatavaluerows(ei.x |> y->QueryOperators.map(y, x->(a=x.a, c=x.c), Expr(:block)))
    @test (mt |> columntable) == (a = Real[1, 2.0, 3], c = ["7", "8", "9"])
    @test length(mt |> rowtable) == 3

    rt = (a = Missing[missing, missing], b=[1,2])
    dv = Tables.datavaluerows(rt)
    @test eltype(dv) == NamedTuple{(:a, :b), Tuple{DataValue{Union{}}, Int}}

    # DataValue{Any}
    @test isequal(Tables.columntable(Tables.nondatavaluerows([(a=DataValue{Any}(), b=DataValue{Int}())])), (a = Any[missing], b = Union{Missing, Int64}[missing]))
end


text:
vents = map(split(read("input.txt", String), '\n')) do line
    parse.(Int, match(r"(\d*),(\d*) -> (\d*),(\d*)", line).captures)
end
ocean_floor = fill(0, (999, 999))

isvert(p) = p[1] == p[3] || p[2] == p[4]
function mark_vents!(ocean_floor, vents)
    for vent in vents
        x₁, y₁, x₂, y₂ = vent
        dx, dy = x₂-x₁, y₂-y₁
        for i in range(0, mapreduce(abs, max, (dx, dy)), step=1)
            ocean_floor[x₁ + i*sign(dx), y₁ + i*sign(dy)] += 1
        end
    end
end

# Part 1
mark_vents!(ocean_floor, filter(isvert, vents))
@info "Horizontal/vertical danger zones: $(count(>(1), ocean_floor))"

# Part 2
mark_vents!(ocean_floor, filter(!isvert, vents))
@info "All danger zones: $(count(>(1), ocean_floor))"

ext: jl


text:
import BioSequences: SkipmerFactory, MerIterResult

function char_complement(::Type{DNA}, c::Char)
    c === 'A' && return 'T'
    c === 'T' && return 'A'
    c === 'C' && return 'G'
    c === 'G' && return 'C'
    error("Oh god, wrong character!")
end

function char_complement(::Type{RNA}, c::Char)
    c === 'A' && return 'U'
    c === 'U' && return 'A'
    c === 'C' && return 'G'
    c === 'G' && return 'C'
    error("Oh god, wrong character!")
end

function string_reverse_complement(::Type{N}, seq::AbstractString) where {N<:NucleicAcid}
    return String([char_complement(N, c) for c in reverse(seq)])
end

function string_iscanonical(::Type{N}, seq::AbstractString) where {N<:NucleicAcid}
    i = 1
    j = length(seq)
    @inbounds while i <= j
        f = seq[i]
        r = char_complement(N, seq[j])
        f < r && return true
        r < f && return false
        i += 1
        j -= 1
    end
    return true
end

function string_canonical(::Type{N}, seq::AbstractString) where {N<:NucleicAcid}
    if string_iscanonical(N, seq)
        return seq
    else
        return string_reverse_complement(N, seq)
    end
end

@testset "Kmers" begin
    function string_eachkmer(::Type{N}, seq::AbstractString, k, step) where {N<:NucleicAcid}
        kmers = Tuple{String,String}[]
        i = 1
        for i in 1:step:length(seq) - k + 1
            subseq = seq[i:i + k - 1]
            if !in('N', subseq)
                push!(kmers, (subseq, string_reverse_complement(N, subseq)))
            end
        end
        return kmers
    end

    function test_eachkmer(S, seq::AbstractString, k, step)
        alph = BioSequences.minimal_alphabet(Alphabet(S))
        mertype = Mer{alph, k}
        iter = each(mertype, S(seq), step)

        @test eltype(iter) == MerIterResult{mertype}

        xs = [(String(x), String(y)) for (i, x, y) in iter]
        ys = string_eachkmer(eltype(S), seq, k, step)

        @test xs == ys
    end

    function test_iteratorlength(S, seq::AbstractString, k, step)
         @test length(each(Mer{BioSequences.minimal_alphabet(Alphabet(S)),k}, S(seq), step)) == length(string_eachkmer(eltype(S), seq, k, step))
    end

    for k in [1, 3, 16, 32], step in 1:3, len in [1, 2, 3, 5, 10, 100, 1000]
        test_eachkmer(LongSequence{DNAAlphabet{4}}, random_dna(len), k, step)
        test_eachkmer(LongSequence{RNAAlphabet{4}}, random_rna(len), k, step)
        test_eachkmer(ReferenceSequence, random_dna(len), k, step)

        probs = [0.25, 0.25, 0.25, 0.25, 0.00]

        test_eachkmer(LongSequence{DNAAlphabet{2}}, random_dna(len, probs), k, step)
        test_eachkmer(LongSequence{RNAAlphabet{2}}, random_rna(len, probs), k, step)

        test_iteratorlength(LongSequence{DNAAlphabet{2}}, random_dna(len, probs), k, step)
        test_iteratorlength(LongSequence{RNAAlphabet{2}}, random_rna(len, probs), k, step)
    end

    @test isempty(collect(each(DNAMer{1}, dna"")))
    @test isempty(collect(each(DNAMer{1}, dna"NNNNNNNNNN")))
    @test_throws Exception each(DNAMer{-1}, dna"ACGT")
    @test_throws Exception each(DNAMer{0}, dna"ACGT")
    @test_throws Exception each(DNAMer{33}, dna"ACGT")
    @test collect(each(DNAMer{3}, dna"AC-TGAG--TGC")) == [MerIterResult(4, DNACodon(DNA_T, DNA_G, DNA_A), DNACodon(DNA_T, DNA_C, DNA_A)),
                                                          MerIterResult(5, DNACodon(DNA_G, DNA_A, DNA_G), DNACodon(DNA_C, DNA_T, DNA_C)),
                                                          MerIterResult(10, DNACodon(DNA_T, DNA_G, DNA_C), DNACodon(DNA_G, DNA_C, DNA_A))]
end

@testset "SkipmerFactory" begin
    seq2 = LongSequence{DNAAlphabet{2}}("GAGGGGGATGCCCCTCTTTGAGCCCAAGG")
    seq4 = LongSequence{DNAAlphabet{4}}("GAGGGGGATGCCCCTCTTTGAGCCCAAGG")

    #function test_iter_traits(it::CanonicalSkipmers{ST,U,S}) where {ST,U,S}
    function test_iter_traits(it::SkipmerFactory{S,UInt64,K}) where {S,K}
        A = ifelse(eltype(S) == DNA, DNAAlphabet{2}, RNAAlphabet{2})
        @test eltype(it) == MerIterResult{Mer{A,K}}
        @test Base.IteratorSize(it) == Base.HasLength()
        @test Base.IteratorEltype(it) == Base.HasEltype()
    end

    function test_iter_traits(it::SkipmerFactory{S,UInt128,K}) where {S,K}
        A = ifelse(eltype(S) == DNA, DNAAlphabet{2}, RNAAlphabet{2})
        @test eltype(it) == MerIterResult{BigMer{A,K}}
        @test Base.IteratorSize(it) == Base.HasLength()
        @test Base.IteratorEltype(it) == Base.HasEltype()
    end

    # Test when Skipmers are based on UInt64, and sample nucletides using
    # M = 2, N = 3, K = 14 pattern.

    #"GAGGGGGATGCCCCTCTTTGAGCCCAAGG"
    #"GA GG GA GC CC CT TG GC CA GG"
    #" AG GG AT CC CT TT GA CC AA G"
    #"G GG GG TG CC TC TT AG CC AG "

    fwans = ["GAGGGAGCCCCTTG", "AGGGATCCCTTTGA", "GGGGTGCCTCTTAG",
             "GGGAGCCCCTTGGC", "GGATCCCTTTGACC", "GGTGCCTCTTAGCC",
             "GAGCCCCTTGGCCA", "ATCCCTTTGACCAA", "TGCCTCTTAGCCAG",
             "GCCCCTTGGCCAGG"]

    M = Mer{DNAAlphabet{2},14}
    BIGM = BigMer{DNAAlphabet{2},14}

    ans = collect(MerIterResult(x...) for x in zip(eachindex(fwans), M.(fwans), M.(string_reverse_complement.(DNA, fwans))))
    bigans = collect(MerIterResult(x...) for x in zip(eachindex(fwans), BIGM.(fwans), BIGM.(string_reverse_complement.(DNA, fwans))))

    #ans = ["CAAGGGGCTCCCTC", "AGGGATCCCTTTGA", "CTAAGAGGCACCCC",
    #       "GCCAAGGGGCTCCC", "GGATCCCTTTGACC", "GGCTAAGAGGCACC",
    #       "GAGCCCCTTGGCCA", "ATCCCTTTGACCAA", "CTGGCTAAGAGGCA",
    #       "CCTGGCCAAGGGGC"]

    #test_iter_traits(eachcanonical(ST, seq2))
    test_iter_traits(SkipmerFactory(M, seq2, 2, 3))
    #test_iter_traits(eachcanonical(ST, seq4))
    test_iter_traits(SkipmerFactory(M, seq4, 2, 3))
    #@test collect(eachcanonical(ST, seq2)) == ST.(ans)
    @test collect(SkipmerFactory(M, seq2, 2, 3)) == ans
    #@test collect(eachcanonical(ST, seq4)) == ST.(ans)
    @test collect(SkipmerFactory(M, seq4, 2, 3)) == ans

    # Test the same pattern, but with UInt128 for big capacity skipmers.
    #test_iter_traits(eachcanonical(BIGST, seq2))
    test_iter_traits(SkipmerFactory(BIGM, seq2, 2, 3))
    #test_iter_traits(eachcanonical(BIGST, seq4))
    test_iter_traits(SkipmerFactory(BIGM, seq4, 2, 3))
    #@test collect(eachcanonical(BIGST, seq2)) == BIGST.(ans)
    @test collect(SkipmerFactory(BIGM, seq2, 2, 3)) == bigans
    #@test collect(eachcanonical(BIGST, seq4)) == BIGST.(ans)
    @test collect(SkipmerFactory(BIGM, seq4, 2, 3)) == bigans

    # Test when Skipmers are based on UInt64, and sample nucleotides using
    # M = 1, N = 3, K = 7 pattern.
    M2 = Mer{DNAAlphabet{2},7}
    BIGM2 = BigMer{DNAAlphabet{2},7}

    fwans2 = ["GGGGCCT", "AGACCTG", "GGTCTTA", "GGGCCTG",
             "GACCTGC", "GTCTTAC", "GGCCTGC", "ACCTGCA",
             "TCTTACA", "GCCTGCG", "CCTGCAG"]

    ans2 = collect(MerIterResult(x...) for x in zip(eachindex(fwans2), M2.(fwans2), M2.(string_reverse_complement.(DNA, fwans2))))
    bigans2 = collect(MerIterResult(x...) for x in zip(eachindex(fwans2), BIGM2.(fwans2), BIGM2.(string_reverse_complement.(DNA, fwans2))))

    # Does iterator stop you if the span of your skipmer exceeds the len of the
    # sequence?
    @test_throws ArgumentError SkipmerFactory(Mer{DNAAlphabet{2},14}, seq2, 1, 3)

    test_iter_traits(SkipmerFactory(M2, seq2, 1, 3))
    test_iter_traits(SkipmerFactory(M2, seq4, 1, 3))
    @test collect(SkipmerFactory(M2, seq2, 1, 3)) == ans2
    @test collect(SkipmerFactory(M2, seq4, 1, 3)) == ans2

    # Test the same pattern, but with UInt128 for big capacity skipmers.
    test_iter_traits(SkipmerFactory(BIGM2, seq2, 1, 3))
    test_iter_traits(SkipmerFactory(BIGM2, seq4, 1, 3))
    @test collect(SkipmerFactory(BIGM2, seq2, 1, 3)) == bigans2
    @test collect(SkipmerFactory(BIGM2, seq4, 1, 3)) == bigans2

    # Test when Skipmers are generated using M = 2, N = 4, K = 8 pattern.
    M3 = Mer{DNAAlphabet{2},8}
    BIGM3 = BigMer{DNAAlphabet{2},8}

    fwans3 = ["GAGGTGCC", "AGGGGCCT", "GGGACCTC", "GGATCCCT",
              "GGTGCCTT", "GGGCCTTT", "GACCTCTG", "ATCCCTGA",
              "TGCCTTAG", "GCCTTTGC", "CCTCTGCC","CCCTGACC",
              "CCTTAGCA", "CTTTGCAA", "TCTGCCAG", "CTGACCGG"]

    ans3 = collect(MerIterResult(x...) for x in zip(eachindex(fwans3), M3.(fwans3), M3.(string_reverse_complement.(DNA, fwans3))))
    bigans3 = collect(MerIterResult(x...) for x in zip(eachindex(fwans3), BIGM3.(fwans3), BIGM3.(string_reverse_complement.(DNA, fwans3))))

    test_iter_traits(SkipmerFactory(M3, seq2, 2, 4))
    test_iter_traits(SkipmerFactory(M3, seq4, 2, 4))
    @test collect(SkipmerFactory(M3, seq2, 2, 4)) == ans3
    @test collect(SkipmerFactory(M3, seq4, 2, 4)) == ans3

    # Test the same pattern, but with UInt128 for big capacity skipmers.
    test_iter_traits(SkipmerFactory(BIGM3, seq2, 2, 4))
    test_iter_traits(SkipmerFactory(BIGM3, seq4, 2, 4))
    @test collect(SkipmerFactory(BIGM3, seq2, 2, 4)) == bigans3
    @test collect(SkipmerFactory(BIGM3, seq4, 2, 4)) == bigans3
end


text:
module TimeSeriesA3P1

#Use this to print
#=
weave(Pkg.dir("$(pwd())\\FamaMacbeth.jl"),
  informat="script",
  out_path = "$(pwd())\\FamaMacbethOUT.html",
  doctype = "md2html")
=#

if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

using DataFrames, Distributions, StatsBase, GZip, JLD,
  Gadfly, JuMP, NLopt, HCubature, StaticArrays, CTNumerical,
  ParallelAccelerator, Measures, Formatting

importall CTIO, CTReg, CTStat

const DATA_PATH = pwd() * "\\data" #path where we store the data
const OUTPUT_PATH = pwd() * "\\output"
const YAHOO_DATE_FORMAT = "m/d/yyyy"
const FRED_DATE_FORMAT = "m/d/yyyy"
const SP500_NAME = "GSPC"
const DBV_NAME = "DBV"
const FRED_NAME = "CPI"

const FOOTER_NAME = "footer.tex" #these are just headers for the tex tables
const HEADER_NAME = "header.tex"
const DECIMALS = 2
const NUM_FORMAT = "{:.$(DECIMALS)f}"
const PlotContainer = Union{Plot,Gadfly.Compose.Context} #holds graphs
const OVERRIDE_STAR_STRINGS = ["\\ostar","\\dstar","\\tstar"]
const COLORS = ["DodgerBlue", "OrangeRed", "Green",
  "Brown", "Black", "BlueViolet"]

#preprocess inflation data
function preProcessFredDF()::Void

  fredDF::DataFrame = readtable("$DATA_PATH\\$FRED_NAME.csv")

  #read and pre-process the data
  dateFormat::String = FRED_DATE_FORMAT
  if :Date ∈ names(fredDF)
    rename!(fredDF, :Date, :DATE)
  elseif :_Date ∈ names(fredDF)
    rename!(fredDF, :_Date, :DATE)
  end

  fredDF[:,:DATE] = ((s::String)->Date(s, dateFormat)::Date).(fredDF[:,:DATE])::DataVector{Date}
  sort!(fredDF, cols = [:DATE])

  fredDF[:,:lCPI] = log.(fredDF[:,:CPI])

  stream::IOStream = open("$DATA_PATH\\$FRED_NAME.jls", "w")
  serialize(stream, fredDF)
  close(stream)

  return nothing
end

#simple container function for the autocorrelation work
function A3BuildAutoρDF(fredDF::DataFrame; lagRange::Range = 1:36)::DataFrame

  #make the dataframe
  ρDF::DataFrame = DataFrame(lag=lagRange, actualρ=autoρ(fredDF[:,:inflation], lagRange))
  ρ1::Float64 = autoρ(fredDF[:,:inflation], 1)
  ρDF[:,:AR1ρ] = (ρ1).^(1:36)

  return ρDF

end

function estimateA3P5Once(sampleMoments::NTuple{4,Float64},
    lowerBound::NTuple{4,Float64}, upperBound::NTuple{4,Float64},
    initialValues::NTuple{4,Float64}, alg::Symbol, timeLimit::Float64)

  function obj(πₐ::A, σ²ᵢ::B, ρ::C, σ²ₓ::D) where {A<:Number, B<:Number, C<:Number, D<:Number}
    oneMinusρSq = (1. - ρ^2)
    SxSi = (σ²ᵢ*σ²ₓ)^0.5
    σ² = (σ²ₓ+σ²ᵢ*oneMinusρSq) / oneMinusρSq
    ρ1 = ((ρ*σ²ₓ+SxSi*oneMinusρSq)/oneMinusρSq)/σ²
    ρ2 = ρ*ρ1

    return (πₐ - sampleMoments[1])^2 + (σ² - sampleMoments[2])^2 +
      (ρ1 - sampleMoments[3])^2 + (ρ2 - sampleMoments[4])^2
  end

  #build the optimizer
  #Algorithms: :LD_LBFGS (fails), :LN_COBYLA (bad), :LN_BOBYQA (bad)
  #  :LD_SLSQP (fails), :LD_TNEWTON (fails), :LD_VAR2/:LD_VAR1 (fails),
  #  :LD_MMA (bad), :GN_ESCH (hang), :GN_DIRECT (questionable), GN_DIRECT_L (questionable),
  # GN_ISRES (bad), :LN_NELDERMEAD(bad), :LN_SBPLX(bad), :LD_LBFGS (fails), :GD_STOGO (fails)
  #GN_CRS2_LM (OK)

  #beginning to solve model
  mod = Model(solver=NLoptSolver(algorithm=alg, maxtime=timeLimit))

  #large boundaries
  @variable(mod, lowerBound[1] ≤ πₐ ≤ upperBound[1])
  @variable(mod, lowerBound[2] ≤ σ²ᵢ ≤ upperBound[2])
  @variable(mod, lowerBound[3] ≤ ρ ≤ upperBound[3])
  @variable(mod, lowerBound[4] ≤ σ²ₓ ≤ upperBound[4])

  setvalue(πₐ, initialValues[1])
  setvalue(σ²ᵢ, initialValues[2])
  setvalue(ρ, initialValues[3])
  setvalue(σ²ₓ, initialValues[4])

  JuMP.register(mod, :obj, 4, obj, autodiff=true)
  @NLobjective(mod, Min, obj(πₐ, σ²ᵢ, ρ, σ²ₓ))

  #println("Running Solve: ")
  output = solve(mod)
  πₐAns::Float64, σ²ᵢAns::Float64, ρAns::Float64, σ²ₓAns::Float64 =
  getvalue(πₐ), getvalue(σ²ᵢ), getvalue(ρ), getvalue(σ²ₓ)

  println("Dennis Obj: $(obj(.0287445, 1.271375E-03^2,.7237045,  2.919860E-03^2))")
  println(" Dennis Inf: $(((1.271375E-03^2+2.919860E-03^2*(1-.7237045^2)) / (1-.7237045^2))^0.5)")

  return (obj(πₐAns, σ²ᵢAns, ρAns, σ²ₓAns), πₐAns, σ²ᵢAns, ρAns, σ²ₓAns)
end

function estimateA3P5(inflation::Vector{Float64}; Δx::Float64 = 0.5,
    lowerBound::NTuple{4,Float64} = (0.0, 0.0, 0.0, 0.0),
    upperBound::NTuple{4,Float64} = (0.01, 0.01, 0.9999999, 0.01),
    initialValue::NTuple{4,Float64} = (upperBound .+ lowerBound)./2,
    hybridGridSearch::Bool = true, alg::Symbol = :LN_SBPLX, timeLimit::Float64 = 2.)

  #first calculate the sample moments
  μ::Float64 = mean(inflation)
  σ²::Float64 = var(inflation)
  ρ1::Float64 = autoρ(inflation, 1)
  ρ2::Float64 = autoρ(inflation, 2)
  sampleMoments::NTuple{4,Float64} = (μ, σ², ρ1, ρ2)
  bestObj::Float64 = 10.^ 10.
  runCtr::Int = 0
  failureCtr::Int = 0

  ΔxAbs::NTuple{4,Float64} = Δx .* (upperBound .- lowerBound)

  if hybridGridSearch
    #this does the grid serach
    for πₐInit::Float64 ∈ (lowerBound[1]+ΔxAbs[1]/2):ΔxAbs[1]:upperBound[1]-ΔxAbs[1]/2,
      σ²ᵢInit::Float64 ∈ (lowerBound[2]+ΔxAbs[2]/2):ΔxAbs[2]:upperBound[2]-ΔxAbs[2]/2,
      ρInit::Float64 ∈ (lowerBound[3]+ΔxAbs[3]/2):ΔxAbs[3]:upperBound[3]-ΔxAbs[3]/2,
      σ²ₓInit::Float64 ∈ (lowerBound[4]+ΔxAbs[4]/2):ΔxAbs[4]:upperBound[4]-ΔxAbs[4]/2

      runCtr += 1

      initialValues::NTuple{4, Float64} = (πₐInit, σ²ᵢInit, ρInit, σ²ₓInit)
      candidateObj::Float64 = 10^10.

      try
        candidateObj,  πₐAns::Float64, σ²ᵢAns::Float64, ρAns::Float64, σ²ₓAns::Float64 =
          estimateA3P5Once(sampleMoments, initialValues .- ΔxAbs./2, initialValues .+ ΔxAbs./2,
            initialValues, alg, timeLimit)
      catch
        failureCtr+=1
      end

      if candidateObj < bestObj
        bestObj=candidateObj
        bestValues::NTuple{4, Float64} = (πₐAns, σ²ᵢAns, ρAns, σ²ₓAns)
      end
    end

  else
    bestObj,  πₐAns, σ²ᵢAns, ρAns, σ²ₓAns =
      estimateA3P5Once(sampleMoments, lowerBound, upperBound, (lowerBound .+ upperBound) ./ 2,
        alg, timeLimit)
    bestValues = (πₐAns, σ²ᵢAns, ρAns, σ²ₓAns)
  end

  println("Simulation Results:")
  println("πₐ = $(bestValues[1]), σᵢ = $(bestValues[2]^0.5), ρ = $(bestValues[3]), σ²ₓ=$(bestValues[4]^0.5)")
  println("Objective: $(bestObj)")
  println("Implied Sigma: $(((bestValues[4]+bestValues[2]*(1-bestValues[3]^2))/(1-bestValues[3]^2))^0.5)")
  println("Sample moments: μ=$μ, σ=$(σ²^0.5), ρ1=$ρ1, ρ2=$ρ2)")
  println("$failureCtr failures out of $runCtr")

  return nothing
end

#problem specific graphs
function visualizeA3P1(fredDF::DataFrame, ρDF::DataFrame)::Void
  plots::Vector{PlotContainer} = Vector{PlotContainer}()
  plotNames::Vector{String} = Vector{String}()
  lagRange::Range = 1:size(ρDF,1)

  Gadfly.push_theme(:default)
  outDF::DataFrame = deepcopy(melt(fredDF[:,[:DATE,:inflationx12,:inflationYoY]], [:DATE]))
  completecases!(outDF)

  push!(plotNames, "Inflation Over Time")
  push!(plots, plot(outDF, x=:DATE, y=:value, color=:variable,  Geom.line,
    Guide.xlabel("Date"),Guide.ylabel("Inflation"),
    Guide.title("Inflation over Time")))

  push!(plotNames, "Inflation auto ρ")
  push!(plots, plot(
    layer(x=collect(lagRange), y=ρDF[:,:actualρ],  Geom.line,
      Theme(default_color=COLORS[1])),
    layer(x=collect(lagRange), y=ρDF[:,:AR1ρ],  Geom.line,
      Theme(default_color=COLORS[2])),
    Guide.manual_color_key("Legend", ["Actual", "AR1"], COLORS[1:2]),
    Guide.xlabel("Lag"),Guide.ylabel("Autocorrelation"),
    Guide.title("AutoCorrelation")))

  #write the graphs
  for i::Int ∈ 1:length(plots)
    draw(SVG("$OUTPUT_PATH\\$(plotNames[i]).svg", 9inch, 7inch), plots[i])
    #display(plots[i])
  end

  return nothing
end

#problem specific code
function A3P1Script(; refreshData::Bool = true)
  fredDF::DataFrame = getDF("$DATA_PATH\\$FRED_NAME.jls", refreshData, preProcessFredDF)

  #some inflation math
  fredDF[:,:inflation] = similar(fredDF[:,:CPI])
  fredDF[2:end,:inflation] = fredDF[2:end,:lCPI] .-  fredDF[1:(end-1),:lCPI]

  fredDF[:,:inflationx12] = fredDF[:,:inflation] .* 12.
  fredDF[:,:inflationYoY] =  similar(fredDF[:,:inflation])
  fredDF[13:end,:inflationYoY] .= fredDF[13:end, :lCPI] .- fredDF[1:(end-12), :lCPI]

  #build the autocorrelation dataframe
  ρDF::DataFrame = A3BuildAutoρDF(fredDF)

  inflation::Vector{Float64} = fredDF[2:end,:inflation]

  #perform the LB test
  TLBTest::Int = round(log(length(inflation)))
  LBTest::UnivariateTest = UnivariateTest(ljungBox, Chisq(TLBTest))

  Q::Float64, p::Float64 = test(LBTest, inflation)

  println("Ljung Box Test Completed with statistic of $Q and p value of $p")

  #problem 5 solve
  #Algorithms: :LD_LBFGS (fails), :LN_COBYLA (bad), :LN_BOBYQA (bad)
  #  :LD_SLSQP (fails), :LD_TNEWTON (fails), :LD_VAR2/:LD_VAR1 (fails),
  #  :LD_MMA (bad), :GN_ESCH (hang), :GN_DIRECT (questionable), GN_DIRECT_L (questionable),
  # GN_ISRES (bad), :LN_NELDERMEAD(bad), :LN_SBPLX(bad), :LD_LBFGS (fails), :GD_STOGO (fails)
  #GN_CRS2_LM (OK)
  #GN_DIRECT .060, LN_SBPLX .03, MMA .278 (Slow, 0.5), CRS2 0.3 (Slow, 0.5), SLSQP 0.08 (Rocky)

  estimateA3P5(inflation, hybridGridSearch=false, Δx = 0.3, alg=:LD_SLSQP, timeLimit = 2.)

  visualizeA3P1(fredDF, ρDF)

end

@time begin
  A3P1Script()
end

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using JLD2
using Statistics

mutable struct Emphasis <: LoggerKeyData
    emphasis::Array{Any}
    log_interval::Int

    function Emphasis(logger_init_info)
        num_logged_steps = fld(logger_init_info[LoggerInitKey.TOTAL_STEPS], logger_init_info[LoggerInitKey.INTERVAL])
        println("The Emphasis logger might be pretty inefficient and should probably only be used for Auto debug. If you are planning on running a sweep, it might be worth not logging step sizes.")
        new([], logger_init_info[LoggerInitKey.INTERVAL])
    end
end

function lg_step!(self::Emphasis, env, agent, s, a, s_next, r, is_terminal, cur_step_in_episode, cur_step_total)
    if rem(cur_step_total, self.log_interval) == 0
        ind = fld(cur_step_total, self.log_interval)
        learner = agent.demon_learner
        lu = update(learner)
        @assert lu isa ETB "Update is not ETB Please set it to ETB (or add my emphasis type here) or remove LoggerKey.EMPHASIS from logger keys"

        emphasis = lu.emphasis_logging
        push!(self.emphasis, copy(emphasis))
    end
end

function lg_episode_end!(self::Emphasis, cur_step_in_episode, cur_step_total)
end

function save_log(self::Emphasis, save_dict::Dict)
    save_dict[:emphasis] = cat(self.emphasis[2:end]..., dims=2)
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
function iscomposedof(cons::Constitution, elems::Vector{<:AbstractString})
    for latt in cons
        for spec in latt
            !iscomposedof(spec, elems) && return false
        end
    end
    return true
end

function isused(func::GFunction, phas::Phase)
    for param in phas.params
        any(str -> occursin(func.name, str), [param.funcstr, join(param.exprs)]) && return true
    end
    return false
end

function isused(func1::GFunction, func2::GFunction)
    return any(str -> occursin(func1.name, str), [func2.funcstr, join(func2.exprs)])
end

function constitutionstring(phas::Phase, s::Int)
    cons = phas.cons
    str = ""
    if length(cons[s]) > 1
        str *= '('
    end
    for spec in cons[s]
        str *= subscript(spec)
        if spec ≠ cons[s][end]
            str *= ','
        end
    end
    if length(cons[s]) > 1
        str *= ')'
    end
    n = phas.sites[s]
    if n ≠ one(n)
        str *= subscript(n)
    end
    return str
end

function constitutionstring(phas::Phase)
    S = length(phas.cons)
    str = ""
    for s in 1:S
        if phas.cons[s] ≠ []
            str *= constitutionstring(phas, s)
        end
    end
    return str
end

function print_func(db::Database)
    for func in db.funcs
        println(func.funcstr)
    end
end


text:
"""Définit la structure d'instance et fournit des méthodes associées."""

"""
Instance du problème.

Args:
    planes (Vector{Plane}): Vecteur des avions
    name (String): Nom de l'instance
    nb_kinds (Int): Nombre de classes d'avions
    nb_planes (Int): Nombre d'avions
    freeze_time (Int):
    sep_mat (Matrix{Int}):
"""
mutable struct Instance
    planes::Vector{Plane}
    name::String
    nb_kinds::Int
    nb_planes::Int
    freeze_time::Int
    sep_mat::Matrix{Int} # tableau d'éléments de type Int et de dimension 2

    function Instance()
        this = new()
        return this
    end

    """
    Lit une instance en fonction du paramètre format
    - Si format vaut alp, ampl ou orlib : ce format est utilisé
    - Si format vaut AUTO alors l'argument de la ligne de commande
      est utilisé :
      - si l'argument impose un format, celui-ci est utilisé
      - sinon (AUTO) le format est deviné à partir de l'instance
    """
    function Instance(infile::AbstractString; format = "AUTO")
        this = new()
        this.name = "NO_NAME"

        if !isfile(infile)
            println("\nERREUR : fichier \"$infile\" inexistant ou illisible !\n")
            exit(1)
        end
        read_alp(this, infile)

        # Mise à jour des attributs dérivés des avions
        # (seulement costs pour l'instant)
        # init_costs!.(this.planes)

        return this
    end
end
# Quelques méthodes auxiliaires simplifiant le codage
sep_max(inst::Instance) = maximum(inst.sep_mat)
ub_max(inst::Instance) = maximum(p->p.ub, inst.planes)
lb_min(inst::Instance) = minimum(p->p.lb, inst.planes)

# Plus grande valeur de pénalité (d'avance ou de retard)
etp_max(inst::Instance) = maximum(p->max(p.ep, p.tp), inst.planes)

"""
Version de haut niveau pour accéder au temps de séparation entre avion.

*Les méthodes de bas niveau peuvent accéder directement au tableau ``inst.sep_mat``.*
"""
function get_sep(inst::Instance, p1::Plane, p2::Plane)
    inst.sep_mat[p1.kind, p2.kind]
end

"""Enregistre l'instance dans un fichier au format demandé."""
function write(inst::Instance, filename::AbstractString; format="alpx")
    fh = open(filename, "w")
    println(fh, to_s_long(inst))
    close(fh)
end

# Quelques fonctions raccourcis
to_s_alp(inst::Instance) = to_s_long(inst, format="alp")
to_s_alpx(inst::Instance) = to_s_long(inst, format="alpx")

"""Génère une chaine au format demandé (ampl, alp, alpx)."""
function to_s_long(inst::Instance)
    io = IOBuffer()
    println(io, "# ALP instance version 1.0")
    println(io)
    println(io, "name ", inst.name)
    println(io, "nb_planes ", inst.nb_planes)
    println(io, "nb_kinds ", inst.nb_kinds)
    println(io, "freeze_time ", inst.freeze_time)
    println(io)
    # println(io, "#    name  kind   at     E     T     L    ep    tp ")
    println(io, to_s_alp_plane_header())
    for p in inst.planes
        println(io, to_s_long(p))
    end
    println(io)

    println(io, "# Separation time between aircraft kinds")
    for k1 = 1:inst.nb_kinds, k2 in 1:inst.nb_kinds
        println(io, "sep ", k1, " ", k2, " ", inst.sep_mat[k1,k2])
    end
    String(take!(io))
end

"""Génère une chaine de statistiques sur l'instance."""
function to_s_stats(inst::Instance)
    io = IOBuffer()
    print(io, "Statistiques sur l'instance\n")
    println(io)
    println(io, "  name: ", inst.name)
    println(io, "  nb_planes: ", inst.nb_planes, " (",
                                 inst.planes[1].name, "..",
                                 inst.planes[end].name, ")")
    println(io, "  nb_kinds: ", inst.nb_kinds)
    println(io, "  freeze_time: ", inst.freeze_time)

    println(io, "="^70)
    println(io, "Caractéristiques des avions\n")
    println(io, "  at:     ", minimum(p->p.at, inst.planes), "..",
                              maximum(p->p.at, inst.planes),)
    println(io, "  lb:     ", minimum(p->p.lb, inst.planes), "..",
                              maximum(p->p.lb, inst.planes),)
    println(io, "  target: ", minimum(p->p.target, inst.planes), "..",
                              maximum(p->p.target, inst.planes),)
    println(io, "  ub:     ", minimum(p->p.ub, inst.planes), "..",
                              maximum(p->p.ub, inst.planes),)
    println(io, "  ub-lb:  ", minimum(p->(p.ub-p.lb), inst.planes), "..",
                              maximum(p->(p.ub-p.lb), inst.planes),)
    println(io, "  ep:     ", minimum(p->p.ep, inst.planes), "..",
                              maximum(p->p.ep, inst.planes),)
    println(io, "  tp:     ", minimum(p->p.tp, inst.planes), "..",
                              maximum(p->p.tp, inst.planes),)

    viols = get_inequality_viols(inst)
    println(io, "-"^70)
    println(io, "Test de get_inequality_viols : nb_viols=$(length(viols))")
    for viol in viols
        (k1, k2, k3, sep12, sep23, sep13) = viol
        print(io,"  $k1->$k3=$sep13 >= ")
        println(io,"$k1->$k2=$(sep12) + $k2->$k3=$(sep23))")
    end
    String(take!(io))
end

"""Retourne l'avion à partir de son nom."""
function get_plane_from_name(inst::Instance, name::AbstractString)
    # TODO : utiliser plutot une fonction de recherche de Julia
    for plane in inst.planes
        if  plane.name == name
            return  plane
        end
    end
    # Cas particulier où les noms d'avion sont de la forme "p1" alors que l'on
    # recherche le nom "1"
    for plane in inst.planes
        if  plane.name == "p$name"
            return  plane
        end
    end
    error("\nAucun avion de nom : $(name)\n")
end

"""
Retourne la liste des viols de l'hypothese de l'inégalité triangulaire

```
sep(k1,k3) > sep(k1,k2) + sep(k2,k3)  => viol existe
```

En effet, l'inégalité triangulaire impose que pour toute paire de types
d'avions (k1,k3), on ne puisse pas insérer un autre avion de type k2 qui
permette de raccourcir l'écart des temps d'atterrissage entre un avion de type
k1 et un avion de type k3.

Un viol est représenté par un tuple de six entiers
`(k1, k2, k3, sep12, sep23, sep13)`
"""
function get_inequality_viols(inst::Instance)
    viols = Vector{Tuple{Int, Int, Int, Int, Int, Int}}()
    if inst.nb_kinds < 2
        # pas de problème d'inégalité triangulaire car moins de 3 types d'avions !
        return viols
    end
    for k1 in 1:inst.nb_kinds, k3 in 1:inst.nb_kinds, k2 in 1:inst.nb_kinds
        sep13 = inst.sep_mat[k1, k3]
        sep12 = inst.sep_mat[k1, k2]
        sep23 = inst.sep_mat[k2, k3]
        if sep13 > sep12+sep23
            viol = (k1, k2, k3, sep12, sep23, sep13)
            push!(viols, viol)
        end
    end
    viols
end


text:
using Printf

include("polynm.jl")

function my_sin(x)
    return sin(10*x)
end

function integ_my_sin(x)
    return -cos(10*x)/10
end

function d_my_sin(x)
    return 10*cos(10*x)
end

function test_main()
    Npoints = 10
    xa = range(1.0, 1.5, length=Npoints)
    ya = my_sin.(xa)

    m = 0
    for i in 1:Npoints
        x = xa[i]
        y = polynm(m, Npoints, xa, ya, x)
        @printf("%18.10f %18.10f %18.10f\n", x, y, ya[i])
    end

    x_plt = range(1.0, 1.5, length=5*Npoints)
    for x in x_plt
        y = polynm(m, Npoints, xa, ya, x)
        #@printf("%18.10f %18.10f %18.10f\n", x, y, my_sin(x))
        #@printf("%18.10f %18.10f %18.10f\n", x, y, integ_my_sin(x))
    end
end

test_main()


text:
# TODO: Build IPython 1.0 dependency? (wait for release?)

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

# print to stderr, since that is where Pkg prints its messages
eprintln(x...) = println(STDERR, x...)

ipyvers = try
    convert(VersionNumber, chomp(readall(`ipython --version`)))
catch e
    error("IPython is required for IJulia, got error $e")
end

if ipyvers < v"1.0.0-dev"
    error("IPython 1.0 or later is required for IJulia, got $ipyvers instead")
else
    eprintln("Found IPython version $ipyvers ... ok.")
end

#######################################################################
# Create Julia profile for IPython and fix the config options.

# create julia profile (no-op if we already have one)
eprintln("Creating julia profile in IPython...")
run(`ipython profile create julia`)

juliaprof = chomp(readall(`ipython locate profile julia`))

# set c.$s in prof file to val, or nothing if it is already set
function add_config(prof::String, s::String, val)
    p = joinpath(juliaprof, prof)
    r = Regex(string("^\\s*c\\.", replace(s, r"\.", "\\."), "\\s*="), "m")
    if isfile(p)
        if ismatch(r, readall(p))
            eprintln("(Existing $s setting in $prof is untouched.)")
        else
            eprintln("Adding $s = $val to $prof...")
            open(p, "a") do f
                print(f, """

c.$s = $val
""")
            end
        end
    else
        eprintln("Creating $prof with $s = $val...")
        open(p, "w") do f
            print(f, """
c = get_config()
c.$s = $val
""")
        end
    end
end

# add Julia kernel manager if we don't have one yet
add_config("ipython_config.py", "KernelManager.kernel_cmd",
           """["$(joinpath(JULIA_HOME,"julia-basic"))", "$(joinpath(Pkg2.dir("IJulia"),"src","kernel.jl"))", "{connection_file}"]""")

# make qtconsole require shift-enter to complete input
add_config("ipython_qtconsole_config.py",
           "IPythonWidget.execute_on_complete_input", "False")

# set Julia notebook to use a different port than IPython's 8888 by default
add_config("ipython_notebook_config.py", "NotebookApp.port", 8998)

#######################################################################
# Copying files into the correct paths in the profile lets us override
# the files of the same name in IPython.

# copy IJulia icon to profile so that IPython will use it
mkpath(joinpath(juliaprof, "static", "base", "images"))
for T in ("png", "svg")
    ipynblogo = joinpath(juliaprof, "static", "base", "images", "ipynblogo.$T")
    if !isfile(ipynblogo)
        eprintln("Copying IJulia $T logo to Julia IPython profile.")
        open(ipynblogo, "w") do f
            write(f, open(readbytes, joinpath(Pkg2.dir("IJulia"), "deps",
                                              "ijulialogo.$T")))
        end
    else
        eprintln("(Existing Julia IPython $T logo file untouched.)")
    end
end

# Use our own version of tooltip to handle identifiers ending with !
# (except for line 211, tooltip.js is identical to the IPython version)
# IPython might make his configurable later, at which point the logic
# should be moved to custom.js or a config file.
mkpath(joinpath(juliaprof, "static", "notebook", "js"))
tooltipjs = joinpath(juliaprof, "static", "notebook", "js", "tooltip.js")
if !isfile(tooltipjs)
    eprintln("Copying tooltip.js to Julia IPython profile.")
    open(tooltipjs, "w") do f
        write(f, open(readbytes, joinpath(Pkg2.dir("IJulia"), "deps",
                                          "tooltip.js")))
    end
else
    eprintln("(Existing tooltip.js file untouched.)")
end

# custom.js can contain custom js login that will be loaded
# with the notebook to add info and/or monkey-patch some javascript
# -- e.g. we use it to add .ipynb metadata that this is a Julia notebook
mkpath(joinpath(juliaprof, "static", "custom"))
customjs = joinpath(juliaprof, "static", "custom", "custom.js")
if !isfile(customjs)
    eprintln("Copying custom.js to Julia IPython profile.")
    open(customjs, "w") do f
        write(f, open(readbytes, joinpath(Pkg2.dir("IJulia"), "deps",
                                          "custom.js")))
    end
else
    eprintln("(Existing custom.js file untouched.)")
end

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


text:
v1 = Point(0,0,0)
v2 = Point(0,0,1)
v3 = Point(1,0,1)
v4 = Point(0,0,2)

@testset "Point distances" begin
    p1 = Point{3, Float32}(0,0,0)
    p2 = Point{3, Float32}(0,0,1)
    @test distance(p1, p2) == 1
end

@testset "Pair selection" begin
    vertices = [v1, v2, v3, v4]
    edges = Set([(1,2), (1,3), (2,3)])
    mesh = MeshAlgorithms.Mesh(vertices, edges, [])
    pairs = MeshAlgorithms.select_valid_pairs(mesh, 1)
    @test pairs == Set([(1,2), (1,3), (2,3), (2, 4)])
end

@testset "Mesh conversion" begin
    f1 = Face(1,2,3)
    normal = Normal(0,1,0)
    test_hm = HomogenousMesh(
            vertices = [v1, v2, v3],
            faces = [f1],
            normals = [normal, normal, normal])
    expected_mesh = Mesh([v1,v2,v3], Set([(1,2), (1,3), (2,3)]), [])
    @test convert_from_homogenous(test_hm).vertices == expected_mesh.vertices
    @test convert_from_homogenous(test_hm).edges == expected_mesh.edges

end

@testset "Error Quadrics" begin
    p = SVector{4}([1,2,3,0])
    k_p = fundamental_error_quadric(p)
    @test k_p[1,1] == p[1]^2
    @test k_p[1,2] == p[1]*p[2]
    @test k_p[1,3] == p[1]*p[3]
    @test k_p[1,4] == p[1]*p[4]
    @test k_p[2,1] == p[1]*p[2]
    @test k_p[2,2] == p[2]^2
    @test k_p[2,3] == p[2]*p[3]
    @test k_p[2,4] == p[2]*p[4]
    @test k_p[3,1] == p[1]*p[3]
    @test k_p[3,2] == p[2]*p[3]
    @test k_p[3,3] == p[3]^2
    @test k_p[3,4] == p[3]*p[4]
    @test k_p[4,1] == p[1]*p[4]
    @test k_p[4,2] == p[2]*p[4]
    @test k_p[4,3] == p[3]*p[4]
    @test k_p[4,4] == p[4]^2
end

sq_sum(v::Vector) = reduce(+, map(x->x^2, v))
@testset "Normalize" begin
    n = [Vector([1,2,3]), Vector([0,3,4]), Vector([1,0,0]), Vector([0,0,0])]
    @test MeshAlgorithms.normalize(n[1]) == [1/sqrt(14), 2/sqrt(14), 3/sqrt(14)]
    @test MeshAlgorithms.normalize(n[2]) == [0,3/5, 4/5]
    @test MeshAlgorithms.normalize(n[3]) == [1,0,0]
    @test_throws ErrorException MeshAlgorithms.normalize(n[4])
    @testset "All length 1 $a $b $c" for a in -5.0:0.5:5, b in -5:0.5:5, c in 0.1:0.2:5
        @test sq_sum(MeshAlgorithms.normalize(Vector([a,b,c]))) ≈ 1.0
    end
end


text:
##############
using Revise
using FDMIC_Geodynamics
using LoopVectorization, Printf, Base.Threads, Plots, Revise, LinearAlgebra, Statistics, SparseArrays
##############
@views function main()
# Domain
xmin      = -1.0
xmax      =  1.0
ymin      = -1.0
ymax      =  1.0
# Scales
Lc        = 1.0
Tc        = 1.0
# Numerics
ncx       = 11
ncy       = 11
dx, dy    = (xmax-xmin)/ncx, (ymax-ymin)/ncy
dt        = 1
# Initialize coordinates
xc        = LinRange( xmin+dx/2, xmax-dx/2, ncx  )
yc        = LinRange( ymin+dy/2, ymax-dy/2, ncy  )
xce       = LinRange( xmin-dx/2, xmax+dx/2, ncx+2)
yce       = LinRange( ymin-dy/2, ymax+dy/2, ncy+2)
xv        = LinRange( xmin     , xmax     , ncx+1)
yv        = LinRange( ymin     , ymax     , ncy+1)
# Parameters
Tamp      = 1.0
sig       = 0.3
k1        = 1.0
k2        = 10.0
steady    = 1
# Allocate tables
T         = zeros( ncx  , ncy  )
T0        = zeros( ncx  , ncy  )
kc        = zeros( ncx  , ncy  )
rho       =  ones( ncx  , ncy  )
Cp        =  ones( ncx  , ncy  )
H         = zeros( ncx  , ncy  )
Vx        = zeros( ncx+1, ncy  )
Vy        = zeros( ncx  , ncy+1)
# Initial configuration
xc2  = repeat(xc, 1, length(yc))
yc2  = repeat(yc, 1, length(xc))'
kc   = k1.*ones( ncx  , ncy  )
kc[ xc2.^2 .+ yc2.^2 .< sig.^2 ] .= k2
xvx2  = repeat(xv, 1, length(yc))
yvx2  = repeat(yc, 1, length(xv))'
xvy2  = repeat(xc, 1, length(yv))
yvy2  = repeat(yv, 1, length(xc))'
Lx, Ly = xmax - xmin, ymax - ymin
ax = 1/2.0; ay = 1/2
Vx    .=  -sin.(ax*2.0 .*pi.*xvx2./Lx) .* cos.(ay*2.0.*pi .* yvx2./Ly)
Vy    .=   cos.(ax*2.0 .*pi.*xvy2./Lx) .* sin.(ay*2.0.*pi .* yvy2./Ly)
Vxc    = 0.5.*(Vx[1:end-1,:] .+ Vx[2:end-0,:])
Vyc    = 0.5.*(Vy[:,1:end-1] .+ Vy[:,2:end-0])
# BC definition
BC = BoundaryConditions()
BC.T.type_W = 0;     BC.T.Dir_W = 0.0*ones(ncy)
BC.T.type_E = 0;     BC.T.Dir_E = 0.0*ones(ncy)
BC.T.type_S = 1;     BC.T.Dir_S = 1.0*ones(ncx)
BC.T.type_N = 1;     BC.T.Dir_N = 0.0*ones(ncx)
# Initialize
kx, ky = AverageConductivity( kc, ncx, ncy, BC.T )
# Residual (-> RHS)
F =  ThermalResidual( T, T0, rho, Cp, H, kx, ky, dx, dy, dt, steady, BC.T, ncx, ncy )
println("Initial residual = ", norm(F))
# Matrix assembly
K = ThermalAssembly( Cp, rho, kx, ky, dx, dy, dt, steady, BC.T, ncx, ncy )
# Solve
ThermalSolve!( T, K, F, ncx, ncy )
# Residual (Check)
F =  ThermalResidual( T, T0, rho, Cp, H, kx, ky, dx, dy, dt, steady, BC.T, ncx, ncy )
println("Solve   residual = ", norm(F))
# Visualize
p1 = Plots.heatmap(xc*Lc, yc*Lc, Array((T*Tc.-0*273.15))', aspect_ratio=1, xlims=(minimum(xv*Lc), maximum(xv)*Lc), ylims=(minimum(yv)*Lc, maximum(yv)*Lc), c=Plots.cgrad(:roma, rev = true), title="T")
# p1 = Plots.heatmap(xc*Lc, yc*Lc, Array((kc*Tc.-0*273.15))', aspect_ratio=1, xlims=(minimum(xv*Lc), maximum(xv)*Lc), ylims=(minimum(yv)*Lc, maximum(yv)*Lc), c=Plots.cgrad(:roma, rev = true), title="T")
display(Plots.plot( p1, dpi=200 ) )
end

for it=1:3
    @time main()
end


text:
using Hungarian

function compute_moment_of_inertia(b,h)
    # input: b,h in meter
    # return: I in cm^4
    return (b*h^3/12)*1e8
end

function compute_strength_capacity(b,h,Fy)
    # input: b,h in meter;
    #        Fy in kN/cm^2: allowable stress
    # return: Fmax in kN
    return b*h*1e4*Fy
end

function compute_buckling_capacity(b,h,L,E)
    # input: b,h,L in meter, E in kN/cm^2
    # return: buckling capacity R in kN
    I = compute_moment_of_inertia(b,h)
    return (1.0/3)*π^2*E*I/L^2*1e-4
end

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

function hungarian_upcycling_assignment(demand_lengths::AbstractVector, supply_lengths::AbstractVector,
    demand_loads::AbstractVector, supply_cross_sec_xs::AbstractVector, supply_cross_sec_ys::AbstractVector,
    allowable_stress::Real, elastic_modulus::Real, penalty::Real)
    demand_num = length(demand_lengths)
    supply_num = length(supply_lengths)

    if demand_num > supply_num ||
       length(demand_loads) != demand_num ||
       length(supply_cross_sec_xs) != supply_num ||
       length(supply_cross_sec_ys) != supply_num
       return nothing, nothing, nothing
    end

    cost_matrix = zeros(Float64, demand_num, supply_num)
    for i in 1:demand_num, j in 1:supply_num
        Li = demand_lengths[i]
        Lj = supply_lengths[j]
        if Li > Lj
            # ! length constraint violated
            cost_matrix[i,j] = penalty
        else
            demand_load = demand_loads[i]
            b = max(supply_cross_sec_xs[j], supply_cross_sec_ys[j])
            h = min(supply_cross_sec_xs[j], supply_cross_sec_ys[j])
            tensile_capacity = compute_strength_capacity(b,h,allowable_stress)
            capacity = 0.0
            if demand_load >= 0
                # tensile demand
                capacity = tensile_capacity
            else
                # compressive demand
                L = min(Li,Lj) # in meter
                buckling_capacity = compute_buckling_capacity(b,h,L,elastic_modulus)
                compressive_capacity = min(buckling_capacity, tensile_capacity)
                capacity = min(tensile_capacity, compressive_capacity)
            end
            if abs(demand_load) > capacity
                # capacity constraint violated
                cost_matrix[i,j] = penalty
            else
                # euclidean distance
                cost_matrix[i,j] = sqrt((Li-Lj)^2+(abs(capacity)-abs(demand_load))^2)
            end
        end
    end
    assignment, cost_star = hungarian(cost_matrix)
    material_coverage = sum(demand_lengths .* supply_cross_sec_xs[assignment] .* supply_cross_sec_ys[assignment]) / sum(supply_lengths .* supply_cross_sec_xs .* supply_cross_sec_ys)
    return assignment, cost_star, material_coverage
end


text:
# Model wrapper for hydrological response model (e.g. Gr4j or HBV)

function run_model(mdata_hydro::Hydro, prec, epot; return_all = false)

  # Number of time steps

  ntimes = size(prec, 2)

  # Allocate output arrays

  q_sim = zeros(Float64, ntimes)

  if return_all == true
    hydro_out = [mdata_hydro for i in 1:ntimes]
  end

  # Run time loop

  for itime in 1:ntimes

    get_input(mdata_hydro, prec, epot, itime)

    run_timestep(mdata_hydro)

    q_sim[itime] = mdata_hydro.q_sim

    if return_all == true
      hydro_out[itime] = deepcopy(mdata_hydro)
    end

  end

  # Return outputs

  if return_all == false
    return q_sim
  else
    return q_sim, hydro_out
  end

end

# Model wrapper for snow model and hydrological response model

function run_model(mdata_snow, mdata_hydro, tair, prec, epot; return_all = false)

  # Number of time steps

  ntimes = size(prec, 2)

  # Allocate output arrays

  q_sim = zeros(Float64, ntimes)

  if return_all == true
    snow_out  = [mdata_snow for i in 1:ntimes]
    hydro_out = [mdata_hydro for i in 1:ntimes]
  end

  # Run time loop

  for itime in 1:ntimes

    get_input(mdata_snow, prec, tair, itime)

    run_timestep(mdata_snow)

    get_input(mdata_snow, mdata_hydro, epot, itime)

    run_timestep(mdata_hydro)

    q_sim[itime] = mdata_hydro.q_sim

    if return_all == true
      snow_out[itime]  = deepcopy(mdata_snow)
      hydro_out[itime] = deepcopy(mdata_hydro)
    end

  end

  # Return outputs

  if return_all == false
    return q_sim
  else
    return q_sim, snow_out, hydro_out
  end

end

# Model wrapper for snow model

function run_model(mdata_snow, tair, prec; return_all = false)

  # Number of time steps

  ntimes = size(prec, 2)

  # Allocate output arrays

  q_sim = zeros(Float64, ntimes)

  if return_all == true
    snow_out  = [mdata_snow for i in 1:ntimes]
  end

  # Run time loop

  for itime in 1:ntimes

    get_input(mdata_snow, prec, tair, itime)

    run_timestep(mdata_snow)

    q_sim[itime] = mdata_snow.q_sim

    if return_all == true
      snow_out[itime]  = deepcopy(mdata_snow)
    end

  end

  # Return outputs

  if return_all == false
    return q_sim
  else
    return q_sim, snow_out
  end

end

max_stars_repo_name: jmgnve/Vann

max_forks_repo_name: jmgnve/Vann


text: include("quadtest.jl")