text:
using Printf

abstract type SubspaceCache{T} end
abstract type HermitianSubspaceCache{T} <: SubspaceCache{T} end

include("stegr_cache.jl")

get_subspace_cache(Ks::KrylovSubspace{B,T,U}) where {B,T,U<:Complex} =
    error("Subspace exponential caches not yet available for non-Hermitian matrices.")
get_subspace_cache(Ks::KrylovSubspace{B,T,U}) where {B,T,U<:Real} =
    StegrCache(Ks)

"""
    expv!(w, t, A, b, Ks, cache)

Alternative interface for calculating the action of `exp(t*A)` on the
vector `b`, storing the result in `w`. The Krylov iteration is
terminated when an error estimate for the matrix exponential in the
generated subspace is below the requested tolerance. `Ks` is a
`KrylovSubspace` and `typeof(cache)<:HermitianSubspaceCache`, the
exact type decides which algorithm is used to compute the subspace
exponential.
"""
function expv!(w::AbstractVector{T}, t::Number, A, b::AbstractVector{T},
               Ks::KrylovSubspace{B, T, B}, cache::HSC;
               atol::B=1.0e-8, rtol::B=1.0e-4,
               m=min(Ks.maxiter, size(A,1)),
               verbose::Bool=false) where {B, T <: Number, HSC <: HermitianSubspaceCache}
    if m > Ks.maxiter
        resize!(Ks, m)
    else
        Ks.m = m # might change if error estimate is below requested tolerance
    end

    V, H = getV(Ks), getH(Ks)
    Ks.beta = norm(b)
    @. V[:, 1] = b / Ks.beta

    ε = atol + rtol * Ks.beta
    verbose && @printf("Initial norm: β₀ %e, stopping threshold: %e\n", Ks.beta, ε)

    α = @diagview(H)
    β = @diagview(H,-1)
    n = size(V, 1)

    for j ∈ 1:m
        lanczos_step!(j, A, V, α, β)
        expT!(@view(α[1:j]), @view(β[1:j]), t, cache)

        # This is practical error estimate Er₂ from
        #
        #   Saad, Y. (1992). Analysis of some Krylov subspace
        #   approximations. SIAM Journal on Numerical Analysis.
        σ = β[j]*Ks.beta*abs(cache.v[j])
        verbose && @printf("iter %d, α[%d] %e, β[%d] %e, σ %e\n",j, j, α[j], j, β[j], σ)
        if σ < ε
            Ks.m = j
            break
        end
    end
    verbose && println("Krylov subspace size: ", Ks.m)

    lmul!(Ks.beta, mul!(w, @view(Ks.V[:,1:Ks.m]), @view(cache.v[1:Ks.m])))
end


text:
using Documenter
using TimeSeries

makedocs(
    format = :html,
    sitename = "TimeSeries.jl",
    modules = [TimeSeries],
    pages = [
        "index.md",
        "getting_started.md",
        "timearray.md",
        "indexing.md",
        "split.md",
        "modify.md",
        "operators.md",
        "apply.md",
        "combine.md",
        "readwrite.md",
        "dotfile.md",
        "plotting.md",
    ]
)

deploydocs(
    repo = "github.com/JuliaStats/TimeSeries.jl.git"
)

ext: jl


text:
include("Delaunay.jl")
using .Delaunay

#TODO:: Needs revision. Doesn't work.
function TriangulationNonSimplicialFaces(VerticesNSFaces,SimplexIndexNS,SimplexFaceIndexNS,ConvexExpIntVert,n)
    #println("\n\n@TriangulationNonSimplicialFaces...\n\n")
    TriangNSFaces = [0]
    numofNSFaces = size(VerticesNSFaces, 2)

    for a = 1:numofNSFaces
        # Contains the indices of the elements in IntVert generating the corresponding
        # polytope face.
        Indices = findall(x->x!=0, VerticesNSFaces[:, a])

        FaceVert = complementary(SimplexFaceIndexNS[a], n + 1) .+ (SimplexIndexNS[a] - 1) * (n + 1)
        BarCoordinates = ConvexExpIntVert[Indices, FaceVert]

        # The rows of this array contain the (strictly positive) concvex exp. coefficients of
        # the IntPoints generating the polytope face, in terms of the vertices
        # generating the simplex sharing the boundary with the polytope face
        i = size(BarCoordinates, 1)
        BarCoordinates = [ones(1, i); transpose(BarCoordinates)]

        # The columns of Null form a basis of the Null space of BarCoordinates
        # PermC is a permutation of the columns of BarCoordinates (or the
        # rows of Null) such that Null(PermC,:)=[M;eye(N)] with M some
        # matrix and N the number of columns in Null.
        Null, PermC = NullSpace(BarCoordinates, n)

        CopBarCoord = zeros(n, i)
        CopBarCoord[:, PermC[1:n]] = Matrix(1.0I, n, n)
        CopBarCoord[:, PermC[(n + 1):i]] = - Null[PermC[1:n], :]

        CoplanarCoord = CopBarCoord[2:n, :]
        CoplanarCoord[:, PermC[1]] .= 0.0

        # Triangulating
        simplices = delaunayn(copy(transpose(CoplanarCoord)))
        # Notice that the order in the argument of delaunayn
        # corresponds to the order given by Indices
        # The rows of simplices contain the indices of the vertices generating
        # the corresponding simplex in the order given by Indices, on the other hand
        # the entries of Indices correspond to the indices of the original intersecting
        # points in IntVert
        d = size(simplices, 1)

        # Triangulation of the non simplicial face a. The indices of the vertices generating
        # the simplices now correspond to the order given by IntVert
        Aux = reshape(transpose(simplices), n * d, 1)
        simplices = transpose(reshape(Indices[Aux], n, d))

        if a == 1
            TriangNSFaces = simplices;
        else
            TriangNSFaces = [TriangNSFaces; simplices]
        end
    end
    return TriangNSFaces
end


text:
module AdventOfCode

mutable struct Node
    name::String
    children::Array{Node}
    parent::Union{Nothing, Node}
end

depth(node::Node) = isnothing(node.parent) ? 0 : (depth(node.parent) + 1)
ancestors(node::Node) = isnothing(node.parent) ? [] : [node.parent; ancestors(node.parent)]

"""Return Dict(name => Node) from list of [parent, child] pairs"""
function build_tree(orbits)
    nodes = Dict()
    for (orbitee, orbiter) in orbits
        orbiter_node = get!(nodes, orbiter, Node(orbiter, [], nothing))
        orbitee_node = get!(nodes, orbitee, Node(orbitee, [], nothing))
        orbiter_node.parent = orbitee_node
        push!(orbitee_node.children, orbiter_node)
    end
    return nodes
end

# Part 1 - build orbit tree and compute depth of each node
orbits = open(f -> split.(readlines(f), ")"), "day06_input.txt")
nodes = build_tree(orbits)
println("Total direct/indirect orbits: $(sum(depth.(values(nodes))))")

# Part 2 - minimum orbital transfers from YOU to SAN
you_nodes = ancestors(nodes["YOU"])
santa_nodes = ancestors(nodes["SAN"])
non_common_nodes = symdiff(you_nodes, santa_nodes)
print("Transfers: $(length(non_common_nodes))")

end

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
import Base.show

function show(x::T where T<:CountModel)
    println("")
    println("Distribution: "*"\t\t"*reduce(*, x.distr.*" "))
    println("Link: "*"\t\t\t"*reduce(*, x.link.*" "))
    if typeof(x) != IIDModel
        if length(x.pastObs) > 0
            println("Past Obs: "*"\t\t"*reduce(*, string.(x.pastObs).*" "))
        end
    end

    if typeof(x) == INGARCHModel
        if length(x.pastMean) > 0
            println("Past Mean: "*"\t\t"*reduce(*, string.(x.pastMean).*" "))
        end
    end

    if length(x.X) > 0
        println("X: "*"\t\t\t"*string(size(x.X)))
    else
        println("X: "*"\t\t\t"*"None")
    end
    if length(x.external) > 0
        println("External: "*"\t\t"*reduce(*, ["N", "Y"][x.external .+ 1].*" "))
    end
    println("Zero Inflation: "*"\t"*ifelse(true, "Y", "N"))
end

function show(x::T where T<:INGARCHModel)
    println("")
    println("Distribution: "*"\t\t"*x.distr)
    println("Link: "*"\t\t\t"*x.link)

    if typeof(x) != IIDModel
        if length(x.pastObs) > 0
            println("Past Obs: "*"\t\t"*reduce(*, string.(x.pastObs).*" "))
        end
    end
    if typeof(x) == INGARCHModel
        if length(x.pastMean) > 0
            println("Past Mean: "*"\t\t"*reduce(*, string.(x.pastMean).*" "))
        end
    end
    if length(x.X) > 0
        println("X: "*"\t\t\t"*string(size(x.X)))
    else
        println("X: "*"\t\t\t"*"None")
    end
    if length(x.external) > 0
        println("External: "*"\t\t"*reduce(*, ["N", "Y"][x.external .+ 1].*" "))
    end
    println("Zero Inflation: "*"\t"*ifelse(true, "Y", "N"))
end

function show(x::T where T<:INARMA)
    println("")
    println("Distribution: "*"\t\t"*x.distr[1]*" and "*x.distr[2])
    println("Link: "*"\t\t\t"*x.link[1]*" and "*x.link[2])
    if typeof(x) != INMAModel
        if length(x.pastObs) > 0
            println("Past Obs: "*"\t\t"*reduce(*, string.(x.pastObs).*" "))
        end
    end
    if typeof(x) != INARModel
        if length(x.pastMean) > 0
            println("Past Mean: "*"\t\t"*reduce(*, string.(x.pastMean).*" "))
        end
    end
    if length(x.X) > 0
        println("X: "*"\t\t\t"*string(size(x.X)))
    else
        println("X: "*"\t\t\t"*"None")
    end
    if length(x.external) > 0
        println("External: "*"\t\t"*reduce(*, ["N", "Y"][x.external .+ 1].*" "))
    end
    println("Zero Inflation: "*"\t"*ifelse(true, "Y", "N"))
end

function show(x::parameter)
    println("")
    println("β0: "*"\t"*string(round(x.β0, digits = 4)))
    if length(x.α) > 0
        println("α: "*"\t"*reduce(*, string.(round.(x.α, digits = 4)).*" "))
    end

    if length(x.β) > 0
        println("β: "*"\t"*reduce(*, string.(round.(x.β, digits = 4)).*" "))
    end

    if length(x.η) > 0
        println("η: "*"\t"*reduce(*, string.(round.(x.η, digits = 4)).*" "))
    end

    if length(x.ϕ) > 0
        println("ϕ: "*"\t"*reduce(*, string.(round.(x.ϕ, digits = 4)).*" "))
    end

    if x.ω > 0
        println("ω: "*"\t"*reduce(*, string.(round.(x.ω, digits = 4)).*" "))
    end
end

function show(x::INGARCHresults)
    tStat = x.θ ./ x.se
    pVals = cdf.(Normal(), -abs.(tStat)).*2
    p1 = ifelse.(pVals .< 0.05, "*", "")
    p2 = ifelse.(pVals .< 0.01, "*", "")
    p3 = ifelse.(pVals .< 0.001, "*", "")
    stars = string.(p1, p2, p3)

    r = length(x.model.external)
    nb = x.model.distr == "NegativeBinomial"
    if typeof(x.model) == INGARCHModel
        name = ["β0"; string.("α", x.model.pastObs); string.("β", x.model.pastMean); string.("η", 1:r); ifelse(nb, "ϕ", []); ifelse(x.model.zi, "ω", [])]
    elseif typeof(x.model) == INARCHModel
        name = ["β0"; string.("α", x.model.pastObs); string.("η", 1:r); ifelse(nb, "ϕ", []); ifelse(x.model.zi, "ω", [])]
    else
        name = ["β0"; string.("η", 1:r); ifelse(nb, "ϕ", []); ifelse(x.model.zi, "ω", [])]
    end

    pr1 = round.(x.θ, digits = 4)
    pr2 = round.(x.se, digits = 4)
    pr3 = round.(pVals, digits = 4)
    pr4 = round.(x.CI[1, :], digits = 4)
    pr5 = round.(x.CI[2, :], digits = 4)
    println("")
    println("Results: Estimates, Standard Errors, p-values, Conf. Intervals")
    for i = 1:x.nPar
        println(name[i], "\t", pr1[i], "\t", pr2[i], "\t", pr3[i], "\t(", pr4[i], ", ", pr5[i], ")", stars[i])
    end
end

function show(x::INARMAresults)
    tStat = x.θ ./ x.se
    pVals = cdf.(Normal(), -abs.(tStat)).*2
    p1 = ifelse.(pVals .< 0.05, "*", "")
    p2 = ifelse.(pVals .< 0.01, "*", "")
    p3 = ifelse.(pVals .< 0.001, "*", "")
    stars = string.(p1, p2, p3)

    r = length(x.model.external)
    nb1 = x.model.distr[1] == "NegativeBinomial"
    nb2 = x.model.distr[2] == "NegativeBinomial"

    if typeof(x.model) == INARMAModel
        name = ["β0"; string.("α", x.model.pastObs); string.("β", x.model.pastMean); string.("η", 1:r); ifelse(nb1, "ϕ1", []); ifelse(nb2, "ϕ2", []); ifelse(x.model.zi, "ω", [])]
    elseif typeof(x.model) == INARModel
        name = ["β0"; string.("α", x.model.pastObs); string.("η", 1:r); ifelse(nb1, "ϕ1", []); ifelse(nb2, "ϕ2", []); ifelse(x.model.zi, "ω", [])]
    elseif typeof(x.model) == INMAModel
        name = ["β0"; string.("β", x.model.pastMean); string.("η", 1:r); ifelse(nb1, "ϕ1", []); ifelse(nb2, "ϕ2", []); ifelse(x.model.zi, "ω", [])]
    else
        name = ["β0"; string.("η", 1:r); ifelse(nb1, "ϕ1", []); ifelse(nb2, "ϕ2", []); ifelse(x.model.zi, "ω", [])]
    end

    pr1 = round.(x.θ, digits = 4)
    pr2 = round.(x.se, digits = 4)
    pr3 = round.(pVals, digits = 4)
    pr4 = round.(x.CI[1, :], digits = 4)
    pr5 = round.(x.CI[2, :], digits = 4)
    println("")
    println("Results: Estimates, Standard Errors, p-values, Conf. Intervals")
    for i = 1:x.nPar
        println(name[i], "\t", pr1[i], "\t", pr2[i], "\t", pr3[i], "\t(", pr4[i], ", ", pr5[i], ")", stars[i])
    end
end


text:
# # Workspace
#
#md # [`notebook`](@__NBVIEWER_ROOT_URL__notebooks/two_problems.ipynb)
#
using WaterWaves1D, Plots, FFTW, Statistics;
gr();
#include("../src/dependencies.jl")

#include("../src/models/WaterWaves.jl")
#include("../src/models/DeepWaterWaves.jl")
#include("../src/models/WWn.jl")

#include("../src/models/SerreGreenNaghdi.jl")
#include("../src/models/Nonhydrostatic.jl")
#include("../src/models/SquareRootDepth.jl")

#include("../src/Figures.jl")
#
#
# if there is an error at this step, try commenting the lines above and commenting out the line below
#include("../src/dependencies.jl")

#---- parameters
param = (
    μ = 0.1,
    ϵ = 0.01,
    N = 2^10, # number of collocation points
    L = 10,# size of the mesh (-L,L)
    T = 5,# final time of computation
    dt = 0.001,  # timestep
);
param = (
    μ = Inf,
    ϵ = 0.1,
    N = 2^12, # number of collocation points
    L = 5,# size of the mesh (-L,L)
    T = 5,# final time of computation
    dt = 0.001,  # timestep
);

#---- initial data
ζ(x) = exp.(-abs.(x).^2 );
v2(x) = zero(x);
init = Init(ζ, v2);
# ζ et v sont des fonctions.
# Init les mets sous une forme utilisable par le programme

#---- models to compare
models = AbstractModel[]
push!(models, WaterWaves(param; dealias = 1))  # dealias = 1 pour dealiasing
#push!(models, Boussinesq(param; dealias = 1))  # dealias = 1 pour dealiasing
#push!(models, WhithamBoussinesq(param; dealias = 1))  # dealias = 1 pour dealiasing
#push!(models, WWn(param; n=2,δ=0.1, dealias = 1, label = "WW2, δ=0.1"))  # dealias = 1 pour dealiasing
push!(models, WWn(param; n=2,δ=0.01, dealias = 1, label = "WW2, δ=0.01"))  # dealias = 1 pour dealiasing
push!(models, WWn(param; n=2,δ=0.001, dealias = 1, label = "WW2, δ=0.001"))  # dealias = 1 pour dealiasing
#push!(models, WWn(param; n=3,δ=0.1, dealias = 1, label = "WW3, δ=0.1"))  # dealias = 1 pour dealiasing
push!(models, WWn(param; n=3,δ=0.01, dealias = 1, label = "WW3, δ=0.01"))  # dealias = 1 pour dealiasing
push!(models, WWn(param; n=3,δ=0.001, dealias = 1, label = "WW3, δ=0.001"))  # dealias = 1 pour dealiasing

problems = Problem[]
for model in models
    push!(problems, Problem(model, init, param))
end

#---- computation
for problem in problems
    @time solve!(problem)
end

#---- visualization

plt = plot_solution(problems[2:end]; fourier = true)
savefig(plt,"plotWWn.pdf") # ou "plot.png" par exemple
pairs=Tuple{Problem, Problem}[]
for i in 2:5
    push!(pairs,(problems[1],problems[i]))
end
plt2 = plot_difference(pairs; fourier = false)
savefig(plt2,"diffWWn.pdf") # ou "plot.png" par exemple

# pour visualiser autre chose, les différences par exemple
#(η1,v1,x1,t)=solution(problems[1])
#(η2,v2,x2,t)=solution(problems[2],x=x1)
#(η3,v3,x3,t)=solution(problems[3],x=x1)

#plt = plot(x1,η1-η2,label="diff 1")
#plot!(plt,x1,η1-η3,label="diff 2")
#plot!(xlims=(0,10))
#plot!(ylims=(-0.1,0.1))

#plot_solution(problems;fourier=false)
# pour faire un film
anim = create_animation(problems;fourier=false,ylims=(-0.5,1.1))
gif(anim, "animWWn.gif", fps=15)
#mov(anim, "SGNvsWGN.mov", fps=15)
#mp4(anim, "SGNvsWGN.mp4", fps=15)


text:
mutable struct MajoranaProduct{T<:Integer}
    #=
    Product of Majorana operators.
    =#

    sites::Vector{T} # Majorana operators on these sites

end

function MajoranaProduct()
    #=
    Construct an empty product (the identity).
    =#

    sites = Int64[]

    MajoranaProduct(sites)

end

function show(io::IO, mp::MajoranaProduct)
    #=
    Formatted printing.
    =#

    if isempty(mp.sites)
        print(io, "I")
    else
        for s in mp.sites
            print(io, "X{", s, "}")
        end
    end

end

length(op::MajoranaProduct) = length(op.sites)

# Is the operator a constant (propto identity)?
isconst(op::MajoranaProduct) = isempty(op.sites)

function iscanonical(op::MajoranaProduct)
    #=
    Check if term is in canonical form.
    =#

    ic = true # in canonical form

    for i in 1:length(op)-1

        if op.sites[i+1] <= op.sites[i] # out of canonical form
            ic = false
            break
        end

    end

    ic

end

function canonicalize!(op::MajoranaProduct)
    #=
    Put product of Majorana operators into canonical form.
    =#

    num = 1 # will produce a factor of +1 or -1 from anticommuting Majorana oparators

    # normal order if not already true
    if !iscanonical(op)

        # bubble sort and keep track of minus signs from commuting majoranas
        l = length(op)

        for i in 1:l
            for j in 1:l-i

                if op.sites[j] > op.sites[j+1]
                    op.sites[j:1:j+1] = op.sites[j+1:-1:j] # swap sites
                    num *= -1 # majoranas anticommute
                end

            end
        end

        # pair annihilate repeated sites (majoranas square to 1)
        keep = trues(l)

        i = 1
        while i < l

            if op.sites[i] == op.sites[i+1]
                keep[i:i+1] .= false
                i += 2
            else
                i += 1
            end

        end

        op.sites = op.sites[keep]

    end

    num

end

function *(op1::MajoranaProduct, op2::MajoranaProduct)
    #=
    Multiply two Majorana products.
    =#

    newsites = [op1.sites; op2.sites]
    newop = MajoranaProduct(newsites)

    # put the new result into canonical form
    num = canonicalize!(newop)

    newop, num

end

function commute(op1::MajoranaProduct, op2::MajoranaProduct)
    #=
    Do these products of Majorana operators commute?
    =#

    # note op12 will equal op21
    op12, num12 = op1 * op2
    op21, num21 = op2 * op1

    # they commute if there's no relative minus sign
    num12 == num21

end

function (==)(op1::MajoranaProduct, op2::MajoranaProduct)
    #=
    Are the operators equal?
    =#

    _equalarrays(op1.sites, op2.sites)

end

max_stars_repo_name: aormorningstar/SymbolicMajoranas

max_issues_repo_name: aormorningstar/SymbolicMajoranas

max_forks_repo_name: aormorningstar/SymbolicMajoranas


text:
loadifpresent(filename, default=Float32[]) = isfile(filename) ? deserialize(filename) : default

"""
    PlotFitness(plotfun, rootdir, subdir="PlotFitness")

Plots best and average fitness for each generation.

Also serializes data so that plotting can be resumed if evolution is aborted.

# Examples
```julia-repl
julia> using NaiveGAflux, Plots

julia> gr();

julia> cb=PlotFitness(plot, "models/test");
```
"""
struct PlotFitness
    best::Vector{Float32}
    avg::Vector{Float32}
    plt
    basedir
end

function PlotFitness(plotfun, rootdir, subdir="PlotFitness")
    basedir = joinpath(rootdir, subdir)
    best = loadifpresent(joinpath(basedir, "best.jls"))
    avg = loadifpresent(joinpath(basedir, "avg.jls"))
    plt = plotfun(collect(Int, 1:length(best)), hcat(best,avg), label=["Best" "Avg"], xlabel="Generation", ylabel="Fitness", m=:circle, legend=:bottomright)
    return PlotFitness(best, avg, plt, basedir)
end

function plotfitness(p::PlotFitness, population)
    fits = fitness.(population)
    best = maximum(fits)
    avg = mean(fits)
    push!(p.best, best)
    push!(p.avg, avg)
    push!(p.plt, length(p.best), [best, avg])
end

plotgen(p::PlotFitness, gen=length(p.best)) = p.plt # Plot already loaded with data...

function (p::PlotFitness)(population)
    plotfitness(p, population)
    mkpath(p.basedir)
    serialize(joinpath(p.basedir, "best.jls"), p.best)
    serialize(joinpath(p.basedir, "avg.jls"), p.avg)
    return p.plt
end

"""
    ScatterPop(plotfun, rootdir, subdir="ScatterPop")

Scatter plot of number of vertices in model vs fitness vs number of parameters for each candidate.

Also serializes data so that plotting can be resumed if evolution is aborted.

# Examples
```julia-repl
julia> using NaiveGAflux, Plots

julia> cb=ScatterPop(scatter, "models/test");
```
"""
struct ScatterPop
    plotfun
    data::Vector{Array{Float32, 2}}
    basedir
end

function ScatterPop(plotfun, rootdir::AbstractString, subdir="ScatterPop")
    basedir = joinpath(rootdir, subdir)
    data = loadifpresent(joinpath(basedir, "nvfitnp.jls"), Vector{Array{Float32, 2}}(undef, 0))
    return ScatterPop(plotfun, data, basedir)
end

function plotfitness(p::ScatterPop, population)
    fits = fitness.(population)
    nverts = nvertices.(model.(population))
    npars = nparams.(population)

    push!(p.data, hcat(nverts, fits, npars))
    plotgen(p)
end

function plotgen(p::ScatterPop, gen=length(p.data))
    gen == 0 && return p.plotfun()
    data = p.data[gen]
    nverts = data[:,1]
    fits = data[:,2]
    npars = data[:,3]
    return p.plotfun(nverts, fits, zcolor=log10.(npars), m=(:heat, 0.8), xlabel="Number of vertices", ylabel="Fitness", colorbar_title="Number of parameters (log10)", label=nothing)
end

function(p::ScatterPop)(population)
    plt = plotfitness(p, population)
    mkpath(p.basedir)
    serialize(joinpath(p.basedir, "nvfitnp.jls"), p.data)
    return plt
end

"""
    ScatterOpt(plotfun, rootdir, subdir="ScatterOpt")

Scatter plot of learning rate vs fitness vs optimizer type for each candidate.

Also serializes data so that plotting can be resumed if evolution is aborted.

# Examples
```julia-repl
julia> using NaiveGAflux, Plots

julia> cb=ScatterOpt(scatter, "models/test");
```
"""
struct ScatterOpt
    plotfun
    data::Vector{Array{Any, 2}}
    basedir
end

function ScatterOpt(plotfun, rootdir::AbstractString, subdir="ScatterOpt")
    basedir = joinpath(rootdir, subdir)
    data = loadifpresent(joinpath(basedir, "fitlropt.jls"), Vector{Array{Any, 2}}(undef, 0))
    return ScatterOpt(plotfun, data, basedir)
end

function plotfitness(p::ScatterOpt, population)
    fits = fitness.(population)
    opts = opt.(population)
    lrs = map(o -> lr(o), opts)
    ots = map(o -> ot(o), opts)

    push!(p.data, hcat(fits, lrs, ots))
    plotgen(p)
end

lr(o) = learningrate(o)
ot(o::Flux.Optimiser) = ot(o.os[1])
ot(o::ShieldedOpt{T}) where T = T
ot(o) = typeof(o)

function plotgen(p::ScatterOpt, gen = length(p.data))
    gen == 0 && return p.plotfun()
    data = p.data[gen]
    fits = data[:,1]
    lrs = data[:,2]
    ots = map(ostr -> last(split(ostr, ".")), string.(data[:,3]))

    uots = sort(unique(ots))
    inds = map(o -> o .== ots, uots)

    fitso = map(indv -> fits[indv], inds)
    lrso = map(indv -> log10.(lrs[indv]), inds)

    return p.plotfun(lrso, fitso, xlabel="Learning rate (log10)", ylabel="Fitness", label=reshape(uots,1,:), legend=:outerright, legendfontsize=5)
end

function(p::ScatterOpt)(population)
    plt = plotfitness(p, population)
    mkpath(p.basedir)
    serialize(joinpath(p.basedir, "fitlropt.jls"), p.data)
    return plt
end

"""
    MultiPlot(plotfun, plts...)

Multiple plots in the same figure.

# Examples
```julia-repl
julia> using NaiveGAflux, Plots

julia> pdir = "models/test"

julia> cb = MultiPlot(display ∘ plot, PlotFitness(plot, pdir), ScatterPop(scatter, pdir), ScatterOpt(scatter, pdir)));
```
"""
struct MultiPlot
    plotfun
    plts
end
function MultiPlot(plotfun, plts...;init=true)
    mp = MultiPlot(plotfun, plts)
    if init
        plotgen(mp)
    end
    return mp
end

(p::MultiPlot)(population) = p.plotfun(map(pp -> pp(population), p.plts)...)

plotgen(p::MultiPlot) = p.plotfun(map(pp -> plotgen(pp), p.plts)...)
plotgen(p::MultiPlot, gen) = p.plotfun(map(pp -> plotgen(pp, gen), p.plts)...)

"""
    struct CbAll
    CbAll(cbs...)

Aggregates several callbacks into one struct.

# Examples
```julia-repl
julia> using NaiveGAflux, Plots

julia> pdir = "models/test"

julia> cb = CbAll(persist, MultiPlot(display ∘ plot, PlotFitness(plot, pdir), ScatterPop(scatter, pdir), ScatterOpt(scatter, pdir)))
```
"""
struct CbAll
    cbs
end
CbAll(cbs...) = CbAll(cbs)

(cba::CbAll)(population) = foreach(cb -> cb(population), cba.cbs)


text:
@testset "552.student-attendance-record-ii.jl" begin
    for (i, j) in zip(
        [
            2
            10958
            17780
            66339
            85599
            82009
            69923
            35625
            77608
            4945
            74685
            62753
            75243
            74415
            42810
            53515
            8347
            28621
            69387
            92734
            38902
            34193
            94858
            70774
            27668
        ],
        [
            8
            649450795
            986891059
            561907615
            364265078
            483789422
            510526871
            108511862
            961167969
            146076425
            30706201
            644547908
            960151554
            675210797
            421499587
            656217676
            473875754
            621142329
            39572738
            618690594
            815584001
            411744682
            24911794
            99736734
            418675188
        ],
    )
        @test check_record(i) == j
    end
end


text:
using Datetime
using Base.Test
y = year(1)
m = months(1)
w = weeks(1)
d = days(1)
h = hour(1)
mi = minute(1)
s = second(1)

#Period arithmetic
@test typeof(y+m) ==  Month{ISOCalendar}
@test typeof(m+y) ==  Month{ISOCalendar}
@test typeof(y+w) ==  Week{ISOCalendar}
@test typeof(y+d) ==  Day{ISOCalendar}
@test typeof(y+h) ==  Hour{ISOCalendar}
@test typeof(y+mi) == Minute{ISOCalendar}
@test typeof(y+s) ==  Second{ISOCalendar}
@test_throws MethodError typeof(m+w)
@test_throws MethodError typeof(m+d)
@test_throws MethodError typeof(m+h)
@test_throws MethodError typeof(m+mi)
@test_throws MethodError typeof(m+s)
@test typeof(w+d) ==  Day{ISOCalendar}
@test typeof(w+h) ==  Hour{ISOCalendar}
@test typeof(w+mi) ==  Minute{ISOCalendar}
@test typeof(w+s) ==  Second{ISOCalendar}
@test typeof(d+h) ==  Hour{ISOCalendar}
@test typeof(d+mi) ==  Minute{ISOCalendar}
@test typeof(d+s) ==  Second{ISOCalendar}
@test typeof(h+mi) ==  Minute{ISOCalendar}
@test typeof(h+s) ==  Second{ISOCalendar}
@test typeof(mi+s) ==  Second{ISOCalendar}
@test y * m*months(2) == months(24)
@test y * w*weeks(2) == weeks(104)
@test y * d*days(2) == days(730)
@test y > m
@test d < w
@test typemax(Year) == 2147483647
@test typemax(Year) + y == -2147483648
@test typemin(Year) == -2147483648
q = convert(Year{Calendar},1)
@test typeof(q) == Year{Calendar}
y,q = promote(y,q) #check periods w/ diff calendars
@test typeof(q) == Year{ISOCalendar}
@test typeof(y) == Year{ISOCalendar}
@test q + y == years(2)
#Period-Real arithmetic
@test y + 1 == 2
@test 1 + y == 2
@test y + 1.0 == 2.0
@test y * 4 == 4
@test y * 4f0 == 4.0f0
@test y * 3//4 == year(1)
@test y / 2 == year(1)
@test 2 / y == 2
@test y / y == 1
@test y*10 % 5 == 0
@test 5 % y*10 == 0
@test !(y > 3)
@test 4 > y
#PeriodRange
t = year(100)
r = y:t
@test size(r) == (100,)
@test length(r) == 100
@test step(r) == year(1)
@test first(r) == y
@test last(r) == t
@test first(r + y) == year(2)
@test last(r + y) == year(101)
@test last(y:year(2):t) == year(99)
@test [year(0):year(25):t][4] == year(75)
@test_throws MethodError y:m:t

#traits
@test zero(typeof(y)) == year(0)
@test one(typeof(y)) == y
@test months(y) == months(12)
@test weeks(y) == weeks(52)
@test days(y) == days(365)

# Wrapping arithemtic for Months
# This ends up being trickier than expected because
# the user might do 2014-01-01 + Month(-14)
# monthwrap figures out the resulting month
# when adding/subtracting months from a date
@test Datetime.monthwrap(1,-14) == 11
@test Datetime.monthwrap(1,-13) == 12
@test Datetime.monthwrap(1,-12) == 1
@test Datetime.monthwrap(1,-11) == 2
@test Datetime.monthwrap(1,-10) == 3
@test Datetime.monthwrap(1,-9) == 4
@test Datetime.monthwrap(1,-8) == 5
@test Datetime.monthwrap(1,-7) == 6
@test Datetime.monthwrap(1,-6) == 7
@test Datetime.monthwrap(1,-5) == 8
@test Datetime.monthwrap(1,-4) == 9
@test Datetime.monthwrap(1,-3) == 10
@test Datetime.monthwrap(1,-2) == 11
@test Datetime.monthwrap(1,-1) == 12
@test Datetime.monthwrap(1,0) == 1
@test Datetime.monthwrap(1,1) == 2
@test Datetime.monthwrap(1,2) == 3
@test Datetime.monthwrap(1,3) == 4
@test Datetime.monthwrap(1,4) == 5
@test Datetime.monthwrap(1,5) == 6
@test Datetime.monthwrap(1,6) == 7
@test Datetime.monthwrap(1,7) == 8
@test Datetime.monthwrap(1,8) == 9
@test Datetime.monthwrap(1,9) == 10
@test Datetime.monthwrap(1,10) == 11
@test Datetime.monthwrap(1,11) == 12
@test Datetime.monthwrap(1,12) == 1
@test Datetime.monthwrap(1,13) == 2
@test Datetime.monthwrap(1,24) == 1
@test Datetime.monthwrap(12,-14) == 10
@test Datetime.monthwrap(12,-13) == 11
@test Datetime.monthwrap(12,-12) == 12
@test Datetime.monthwrap(12,-11) == 1
@test Datetime.monthwrap(12,-2) == 10
@test Datetime.monthwrap(12,-1) == 11
@test Datetime.monthwrap(12,0) == 12
@test Datetime.monthwrap(12,1) == 1
@test Datetime.monthwrap(12,2) == 2
@test Datetime.monthwrap(12,11) == 11
@test Datetime.monthwrap(12,12) == 12
@test Datetime.monthwrap(12,13) == 1

# yearwrap figures out the resulting year
# when adding/subtracting months from a date
@test Datetime.yearwrap(2000,1,-3600) == 1700
@test Datetime.yearwrap(2000,1,-37) == 1996
@test Datetime.yearwrap(2000,1,-36) == 1997
@test Datetime.yearwrap(2000,1,-35) == 1997
@test Datetime.yearwrap(2000,1,-25) == 1997
@test Datetime.yearwrap(2000,1,-24) == 1998
@test Datetime.yearwrap(2000,1,-23) == 1998
@test Datetime.yearwrap(2000,1,-14) == 1998
@test Datetime.yearwrap(2000,1,-13) == 1998
@test Datetime.yearwrap(2000,1,-12) == 1999
@test Datetime.yearwrap(2000,1,-11) == 1999
@test Datetime.yearwrap(2000,1,-2) == 1999
@test Datetime.yearwrap(2000,1,-1) == 1999
@test Datetime.yearwrap(2000,1,0) == 2000
@test Datetime.yearwrap(2000,1,1) == 2000
@test Datetime.yearwrap(2000,1,11) == 2000
@test Datetime.yearwrap(2000,1,12) == 2001
@test Datetime.yearwrap(2000,1,13) == 2001
@test Datetime.yearwrap(2000,1,23) == 2001
@test Datetime.yearwrap(2000,1,24) == 2002
@test Datetime.yearwrap(2000,1,25) == 2002
@test Datetime.yearwrap(2000,1,36) == 2003
@test Datetime.yearwrap(2000,1,3600) == 2300
@test Datetime.yearwrap(2000,2,-2) == 1999
@test Datetime.yearwrap(2000,3,10) == 2001
@test Datetime.yearwrap(2000,4,-4) == 1999
@test Datetime.yearwrap(2000,5,8) == 2001
@test Datetime.yearwrap(2000,6,-18) == 1998
@test Datetime.yearwrap(2000,6,-6) == 1999
@test Datetime.yearwrap(2000,6,6) == 2000
@test Datetime.yearwrap(2000,6,7) == 2001
@test Datetime.yearwrap(2000,6,19) == 2002
@test Datetime.yearwrap(2000,12,-3600) == 1700
@test Datetime.yearwrap(2000,12,-36) == 1997
@test Datetime.yearwrap(2000,12,-35) == 1998
@test Datetime.yearwrap(2000,12,-24) == 1998
@test Datetime.yearwrap(2000,12,-23) == 1999
@test Datetime.yearwrap(2000,12,-14) == 1999
@test Datetime.yearwrap(2000,12,-13) == 1999
@test Datetime.yearwrap(2000,12,-12) == 1999
@test Datetime.yearwrap(2000,12,-11) == 2000
@test Datetime.yearwrap(2000,12,-2) == 2000
@test Datetime.yearwrap(2000,12,-1) == 2000
@test Datetime.yearwrap(2000,12,0) == 2000
@test Datetime.yearwrap(2000,12,1) == 2001
@test Datetime.yearwrap(2000,12,11) == 2001
@test Datetime.yearwrap(2000,12,12) == 2001
@test Datetime.yearwrap(2000,12,13) == 2002
@test Datetime.yearwrap(2000,12,24) == 2002
@test Datetime.yearwrap(2000,12,25) == 2003
@test Datetime.yearwrap(2000,12,36) == 2003
@test Datetime.yearwrap(2000,12,37) == 2004
@test Datetime.yearwrap(2000,12,3600) == 2300

#Date tests
dt = date(1,1,1)
@test convert(Date,1) == dt
@test convert(Date,int32(1)) == dt
@test convert(Date,float32(1)) == dt
dt2 = convert(Date{Calendar},float64(1))
@test typeof(dt2) == Date{Calendar}
a,b = promote(dt,dt2) #checks 2 diff calendars promotion
@test typeof(a) == Date{ISOCalendar}
@test typeof(b) == Date{ISOCalendar}
@test dt == dt2
@test dt - dt2 == days(0)
@test !(dt > dt2)
@test length(dt) == 1
@test convert(Date,0) == date(year(0),month(12),day(31))
@test convert(Date,-1) == date(0,12,30)
@test_throws MethodError dt2 + y #Diff calendars
dt3 = datetime(1,1,1)
@test convert(DateTime,86400000) == dt3
dt4 = convert(DateTime{Calendar,UTC},float64(86400000))
@test typeof(dt4) == DateTime{Calendar,UTC}
a,b = promote(dt3,dt4) #checks 2 diff calendars promotion
@test typeof(a) == DateTime{ISOCalendar,UTC}
@test typeof(b) == DateTime{ISOCalendar,UTC}
@test dt3 == dt4
@test dt3 - dt4 == seconds(0)
@test !(dt3 > dt4)
@test length(dt3) == 1
a,b = promote(dt,dt3) #promote Date => DateTime
@test typeof(a) == DateTime{ISOCalendar,UTC}
@test typeof(b) == DateTime{ISOCalendar,UTC}
@test date(dt3) == dt
@test datetime(dt) == dt3
@test dt == dt3
@test !(dt > dt3)
@test dt - dt3 == seconds(0)

check = (-252522163911150,-1,0,1,252522163911149)
for (i,x) in enumerate(check)
	yd = Datetime.yeardays(x-1)
	md = Datetime.monthdays[1]
	d = 1
	td = Datetime.totaldays(x,1,1)
	@test d + md + yd - 306 == td
	@test td == (-92231826452318568,-730,-365,1,92231826452317474)[i]
	@test Datetime._year(td) == x
	@test Datetime._month(td) == 1
	@test Datetime._day(td) == 1
	@test Datetime._day2date(td) == (x,1,1)
	@test Datetime._week(td) == (1,1,1,1,53)[i]
end
@test isleap(-4)
@test isleap(0)
@test isleap(4)
@test !isleap(100)
@test isleap(400)
@test isleap(2000)
@test !isleap(2100)

dt = date(2013,7,1)
@test year(dt) == 2013
@test month(dt) == 7
@test day(dt) == 1
@test Datetime._days(dt) == 735050 == int64(dt)
@test string(dt) == "2013-07-01"
@test week(dt) == 27
@test isdate(dt)
@test calendar(dt) == ISOCalendar
@test typemax(Date) == date(252522163911149,12,31)
@test typemin(Date) == date(-252522163911150,1,1)
@test isleap(dt) == false
@test lastdayofmonth(dt) == 31
@test dayofweek(dt) == 1
@test dayofweekinmonth(dt) == 1
@test daysofweekinmonth(dt) == 5
@test firstdayofweek(dt) == dt
@test lastdayofweek(dt) == dt + days(6)
@test dayofyear(dt) == 182
@test typeof(today()) <: Date
dt = date(2012,2,29)
dt2 = date(2000,2,1)
@test dt > dt2
@test dt != dt2
@test year(-dt) == -2012
@test +dt == dt
@test dt - dt2 == 4411
@test dt > dt2
@test dt2 < dt
@test dt2 - year(3) == date(1997,2,1)
@test dt + year(1) == date(2013,2,28)
@test year(3) - dt2 == dt2 - year(3)
@test dt2 - months(3) == date(1999,11,1)
@test months(11) + dt == date(2013,1,29)
@test months(8) + dt == dt + months(8)
@test dt2 + days(4411) == date(2012,2,29)
@test dt2 + days(4412) == date(2012,3,1)
@test dt2 + weeks(52) == date(2001,1,30)
@test dt2 + weeks(104) == date(2002,1,29)
@test dt2 - year(4) + days(366) == date(1997,2,1)
@test dt2 - (year(4) - days(366)) == date(1997,2,2) #non-associative
dt = date(1999,12,27)
check = (52,52,52,52,52,52,52,1,1,1,1,1,1,1,2,2,2,2,2,2,2)
for i = 1:21
	@test week(dt) == check[i]
	dt = dt + day(1)
end
dt = date(2000,12,25)
for i = 1:21
	@test week(dt) == check[i]
	dt = dt + day(1)
end
dt = date(2030,12,23)
for i = 1:21
	@test week(dt) == check[i]
	dt = dt + day(1)
end
dt = date(2004,12,20)
check = (52,52,52,52,52,52,52,53,53,53,53,53,53,53,1,1,1,1,1,1,1)
for i = 1:21
	@test week(dt) == check[i]
	dt = dt + day(1)
end
#Simple Date parsing
a = "2005-01-31"
b = "2005/31/01"
c = "01.31.2005"
d = "31 01 2005"
aa = "2005,1,31"
bb = "2005-31-1"
cc = "1-31-2005"
dd = "31-1-2005"
aaa = "05-01-31"
bbb = "05/31/01"
ccc = "01.31.05"
ddd = "31 01 05"
e = "20050131"
ee = "2005131"
for i in (a,b,c,d,aa,bb,cc,dd,ccc,ddd)
	@test date(i) == date(2005,1,31)
end
@test date(aaa) == date(2031,5,1)
@test date(bbb) == date(2001,5,31)

#DateRange
dt1 = date(2000,1,1)
dt2 = date(2010,1,1)
y = year(1)
r = dt1:y:dt2
@test first(r) == date(2000,1,1)
@test last(r) == date(2010,1,1)
@test typeof(r.step) == Year{ISOCalendar}
@test length(r) == 11
dt2 = date(2001,1,1)
r = dt1:month(1):dt2
@test length(r) == 13
@test last(r) == date(2001,1,1)
@test typeof(r.step) == Month{ISOCalendar}
r = dt1:weeks(2):dt2
@test length(r) == 27
@test last(r) == date(2000,12,30)
@test typeof(r.step) == Week{ISOCalendar}
dt2 = date(2000,3,1)
r = dt2:day(-1):dt1
@test length(r) == 61
@test last(r) == date(2000,1,1)
@test typeof(r.step) == Day{ISOCalendar}
@test last(r + year(1)) == date(2000,12,31)
dt = date(2009,1,1)
dt2 = date(2013,1,1)
MemorialDay = x->month(x)==May && dayofweek(x)==Monday && dayofweekinmonth(x)==daysofweekinmonth(x)
f = dt:MemorialDay:dt2
@test typeof(f) == DateRange1{ISOCalendar}
@test [f] == [date(2009,5,25),date(2010,5,31),date(2011,5,30),date(2012,5,28)]
@test length(f) == 4
@test last(f) == date(2012,5,28)
@test typeof(f.step) == Function
Thanksgiving = x->dayofweek(x)==Thursday && month(x)==November && dayofweekinmonth(x)==4
f = dt:Thanksgiving:dt2
@test [f] == [date(2009,11,26),date(2010,11,25),date(2011,11,24),date(2012,11,22)]
@test length(f) == 4
@test last(f) == date(2012,11,22)
ff = recur(dt:dt2) do x
	dayofweek(x)==4 &&
	month(x)==11 &&
	dayofweekinmonth(x)==4
end
@test [ff] == [f]
#DateTime
dt = datetime(2013,7,1)
@test year(dt) == 2013
@test month(dt) == 7
@test day(dt) == 1
@test hour(dt) == 0
@test minute(dt) == 0
@test second(dt) == 0
@test string(dt) == "2013-07-01T00:00:00 UTC"
@test Datetime._days(dt) == 735050 == fld(int64(dt),86400000)
@test week(dt) == 27
@test isdatetime(dt)
@test calendar(dt) == ISOCalendar
@test timezone(dt) == Datetime.Zone0
@test typemax(DateTime) == datetime(292277024,12,31,23,59,59)
@test typemin(DateTime) == datetime(-292277023,1,1,0,0,0)
@test isleap(dt) == false
@test lastdayofmonth(dt) == 31
@test dayofweek(dt) == 1
@test dayofweekinmonth(dt) == 1
@test daysofweekinmonth(dt) == 5
@test firstdayofweek(dt) == dt
@test lastdayofweek(dt) == dt + days(6)
@test dayofyear(dt) == 182
@test typeof(now()) <: DateTime
@test int64(dt) == 63508320025000
@test second(datetime(1972,6,30,22,58,60)) == 0 #entering "invalid" periods just rolls the date forward
dt = datetime(2012,2,29)
dt2 = datetime(2000,2,1)
@test dt > dt2
@test dt != dt2
@test isdatetime(dt)
@test calendar(dt) == ISOCalendar
@test year(-dt) == -2012
@test +dt == dt
@test dt - dt2 == 381110402000
@test dt > dt2
@test dt2 < dt
@test dt2 - year(3) == datetime(1997,2,1)
@test dt + year(1) == datetime(2013,2,28)
@test year(3) - dt2 == dt2 - year(3)
@test dt2 - months(3) == datetime(1999,11,1)
@test months(11) + dt == datetime(2013,1,29)
@test months(8) + dt == dt + months(8)
@test dt2 + days(4411) == datetime(2012,2,29)
@test dt2 + days(4412) == datetime(2012,3,1)
@test dt2 + weeks(52) == datetime(2001,1,30)
@test dt2 + weeks(104) == datetime(2002,1,29)
@test dt2 - year(4) + days(366) == datetime(1997,2,1)
@test dt2 - (year(4) - days(366)) == datetime(1997,2,2) #non-associative
dt = datetime(1999,12,27,0,0,0)
check = (52,52,52,52,52,52,52,1,1,1,1,1,1,1,2,2,2,2,2,2,2)
for i = 1:21
	@test week(dt) == check[i]
	dt = dt + day(1)
end
dt = datetime(2000,12,25,0,0,0)
for i = 1:21
	@test week(dt) == check[i]
	dt = dt + day(1)
end
dt = datetime(2030,12,23,0,0,0)
for i = 1:21
	@test week(dt) == check[i]
	dt = dt + day(1)
end
dt = datetime(2004,12,20,0,0,0)
check = (52,52,52,52,52,52,52,53,53,53,53,53,53,53,1,1,1,1,1,1,1)
for i = 1:21
	@test week(dt) == check[i]
	dt = dt + day(1)
end

#Leap second support
# _lm =
# [1972 6 30 23 59 59;
# 1972 12 31 23 59 59;
# 1973 12 31 23 59 59;
# 1974 12 31 23 59 59;
# 1975 12 31 23 59 59;
# 1976 12 31 23 59 59;
# 1977 12 31 23 59 59;
# 1978 12 31 23 59 59;
# 1979 12 31 23 59 59;
# 1981 6 30 23 59 59;
# 1982 6 30 23 59 59;
# 1983 6 30 23 59 59;
# 1985 6 30 23 59 59;
# 1987 12 31 23 59 59;
# 1989 12 31 23 59 59;
# 1990 12 31 23 59 59;
# 1992 6 30 23 59 59;
# 1993 6 30 23 59 59;
# 1994 6 30 23 59 59;
# 1995 12 31 23 59 59;
# 1997 6 30 23 59 59;
# 1998 12 31 23 59 59;
# 2005 12 31 23 59 59;
# 2008 12 31 23 59 59;
# 2012 6 30 23 59 59]
# #Generate leap moments in microsends
# _leapmoments = ref(Int64) #1 sec before a leap second in a timeline that ignores leap seconds
# _leapmoments1 = ref(Int64) #leap second moment in leap second timeline
# for i in 1:size(_lm)[1]
# 	l = 1000*(_lm[i,6] +
# 		60*_lm[i,5] + 3600*_lm[i,4] +
# 		86400*totaldays(_lm[i,1],_lm[i,2],_lm[i,3]))
# 	push!(_leapmoments,l)
# 	push!(_leapmoments1,l+(1000*length(_leapmoments)))
# end
q = datetime(1972,6,30,23,59,58)
t = datetime(1972,6,30,23,59,59)
@test t - q == 1000
@test t == t
r = datetime(1972,7,1,0,0,0)
@test r - t == 2000 #Includes leap second
n = datetime(1972,6,30,23,59,60)
@test n - t == 1000
@test r - n == 1000
a = datetime(1972,12,31,23,59,59)
b = datetime(1972,12,31,23,59,60)
c = datetime(1973,1,1,0,0,0)
d = datetime(1973,1,1,0,0,1)
@test c - b == 1000
@test b - a == 1000
a = datetime(1982,6,30,23,59,59)
b = datetime(1982,6,30,23,59,60)
c = datetime(1982,7,1,0,0,0)
@test c - b == 1000
@test b - a == 1000

#Offsets
TZ = CST
a = datetime(1972,12,31,17,59,59,0,TZ)
b = datetime(1972,12,31,17,59,60,0,TZ)
c = datetime(1972,12,31,18,0,0,0,TZ)
d = datetime(1972,12,31,18,0,1,0,TZ)
@test c - b == 1000
@test b - a == 1000
@test c - a == 2000 #Includes leap second
bb = timezone(b,PST)
@test bb == b
bb = timezone(b,"Pacific/Honolulu")
@test bb == b
bb = offset(b,Offset{-360})
@test bb == b
@test string(bb) == "1972-12-31T17:59:60 -06:00"
aa = offset(a,Offset{-360})
@test aa == a
@test offset(a,Offset{-600}) == timezone(a,"Pacific/Honolulu")
@test offset(a,offset(hours(10))) == timezone(a,"Pacific/Honolulu")
@test timezone(datetime(2013,7,6,0,0,0,0,"America/Chicago")) == CDT

#Daylight savings time
t1 = datetime(1920,6,13,1,59,59,0,CST)
t2 = datetime(1920,6,13,2,0,0,0,CST)
t3 = datetime(1920,6,13,3,0,0,0,CST)
@test t2 - t1 == 1000
@test t3 - t1 == 1000
@test t3 == t2
@test t1 + second(1) == t2
@test t1 + second(1) == t3
t4 = datetime(1920,6,13,4,0,0,0,CST)
@test t4 - hour(1) == t3

#Conversion from Unix time
# _ = time()
# tt = unix2datetime(1000*int64(_),TZ)
# ttt = TmStruct(_)
# @test year(tt) == int("20"*string(digits(ttt.year)[2])*string(digits(ttt.year)[1]))
# @test month(tt) == ttt.month+1
# @test day(tt) == ttt.mday
# @test hour(tt) == ttt.hour
# @test minute(tt) == ttt.min
# @test second(tt) == ttt.sec #second may be off by 1, not sure why
#DateTimeRange
dt1 = datetime(2000,1,1)
dt2 = datetime(2010,1,1)
y = year(1)
r = dt1:y:dt2
@test first(r) == datetime(2000,1,1)
@test last(r) == datetime(2010,1,1)
@test typeof(r.step) == Year{ISOCalendar}
@test length(r) == 11
dt2 = datetime(2001,1,1)
r = dt1:month(1):dt2
@test length(r) == 13
@test last(r) == datetime(2001,1,1)
@test typeof(r.step) == Month{ISOCalendar}
r = dt1:weeks(2):dt2
@test length(r) == 27
@test last(r) == datetime(2000,12,30)
@test typeof(r.step) == Week{ISOCalendar}
dt2 = datetime(2000,3,1)
r = dt2:day(-1):dt1
@test length(r) == 61
@test last(r) == datetime(2000,1,1)
@test typeof(r.step) == Day{ISOCalendar}
@test last(r + year(1)) == datetime(2000,12,31)
dt2 = datetime(2000,1,31)
r = dt1:hours(24):dt2
@test length(r) == 31
@test last(r) == datetime(2000,1,31)
@test typeof(r.step) == Hour{ISOCalendar}
r = dt1:minutes(1440):dt2
@test length(r) == 31
@test last(r) == datetime(2000,1,31)
@test typeof(r.step) == Minute{ISOCalendar}
r = dt1:seconds(86400):dt2
@test length(r) == 31
@test last(r) == datetime(2000,1,31)
@test typeof(r.step) == Second{ISOCalendar}

#Test for 32 bit issues #20
datetime(int32(2013),int32(10),int32(27),int32(11),int32(10),int32(9))
datetime(int32(2013),int32(10),int32(27),int32(11),int32(10),int32(9), int32(8), timezone("UTC"))

#Test for date->DateTime conversion w/ leap seconds
@test date(datetime(date(2012,7,1))) == date(2012,7,1)

#Test for date->DateTime conversion w/ Timezone
@test datetime(date(2012,7,1),EST) == datetime(2012,7,1,4,0,0,0)
@test datetime(date(2012,7,1),"America/New_York") == datetime(2012,7,1,4,0,0,0)


text:
using BoardGameGeek
using Dates

name = "bggjulia"

user = userinfo(name)
@test user.id == 3162756
@test user.name == name
@test user.country == "Poland"
@test user.yearregistered == 2022

buds = buddies(name)
@test length(buds) == 2

reviews = userreviews(name)
@test length(reviews) == 4

r = first(sort(reviews; by=r -> -r.rating))
@test r.id == 188
@test r.name == "Go"
@test r.username == "bggjulia"
@test r.rating == 10.0
@test r.comment == ""

vuser = userinfo([name])
vbuds = buddies([name])
vreviews = userreviews([name])
@test first(vuser) == user
@test first(vbuds) == buds
@test vreviews == reviews


text:
export AbstractMeterProvider,
    AbstractInstrument, global_meter_provider, global_meter_provider!, Meter

"""
A meter provider defines how to collect and update [`Measurement`](@ref)s. Each
meter provider should have the following interfaces implemented:

  - `Base.push!(provider, m::Meter)`, register a meter.
  - `Base.push!(provider, ins::AbstractInstrument)`, register an instrument.
  - `Base.push!(provider, (ins::AbstractInstrument, m::Measurement))`, update a measurement.
"""
abstract type AbstractMeterProvider end

Base.@kwdef struct DummyMeterProvider <: AbstractMeterProvider end

Base.push!(p::DummyMeterProvider, x) = nothing

resource(::DummyMeterProvider) = Resource()

const GLOBAL_METER_PROVIDER = Ref{AbstractMeterProvider}(DummyMeterProvider())

"""
    global_meter_provider()

Get the global meter provider.
"""
global_meter_provider() = GLOBAL_METER_PROVIDER[]

"""
    global_meter_provider!(p)

Set the global meter provider to `p`.
"""
global_meter_provider!(p) = GLOBAL_METER_PROVIDER[] = p

"""
`AbstractInstrument`` is the super type of all instruments, which are used to report [`Measurement`](@ref)s.

Each concrete subtype should at least have the following fields:

  - `name`
  - `description`
  - `unit`
  - `meter`, the associated [`Meter`](@ref)

See also [the specification](https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/metrics/api.md#instrument):
"""
abstract type AbstractInstrument{T} end

provider(ins::AbstractInstrument) = provider(ins.meter)
resource(ins::AbstractInstrument) = resource(ins.meter)

"""
    Meter(name::String;kw...)

Meter is responsible for creating instruments.

## Keyword Arguments:

  - `provider::P = global_meter_provider()`
  - `instrumentation_info = InstrumentationInfo()`
"""
struct Meter{P<:AbstractMeterProvider}
    name::String
    provider::P
    instrumentation_info::InstrumentationInfo
    instruments::Vector{AbstractInstrument}
    function Meter(
        name::String;
        provider::P = global_meter_provider(),
        instrumentation_info = InstrumentationInfo(),
    ) where {P<:AbstractMeterProvider}
        m = new{P}(name, provider, instrumentation_info, AbstractInstrument[])
        push!(provider, m)
        m
    end
end

function Base.push!(m::Meter, ins::AbstractInstrument)
    push!(m.instruments, ins)
    push!(m.provider, ins)
end

Base.push!(m::Meter, ms::Pair{<:AbstractInstrument}) = push!(m.provider, ms)

provider(m::Meter) = m.provider

resource(m::Meter) = resource(provider(m))


text:
draw(PNG("square-bevel.png", 140, 100),
     Axes(border=0.2),
     Paper("lightgrey"),
     Pen(0.3, cap="square", join="bevel")) do
    move(0, 0)
    line(0.5, 0)
    line(0.5, 0.5)
    line(1, 1)
end


text:
function line_search(f, x0, g, h, opt_type="steepest", max_iter=50, tol=1e-6)
    """
    f: objective function
    x0: starting value of iterate
    grad_f: gradient of the objective
    hess_f: Hessian of the objective, not required for all methods
    opt_type: type of step direction
        "steepest"
        "newton"
        "newton_CG"
    """
    println("Using method ", opt_type)
    n = length(x0)
    xvals = Vector{Vector{Float64}}(); push!(xvals, x0)
    ncalls = 0; Bk = nothing

    for i in 1:max_iter
        i % 100 == 0 && println(i)
        xcurr = xvals[end]
        xprev = i > 1 ? xvals[end-1] : nothing

        p, alpha, nc, Bk = compute_steps(xcurr, f, g, h, opt_type, Bk, xprev)
        xnew = xcurr + alpha * p
        push!(xvals, xnew)
        ncalls += nc

        if abs(mean(xnew - xcurr)) <= tol
            println("Number of indefinite fixes ", ncalls)
            println("Number of iterations ", i)
            return(xvals)
        end
    end

    println("Finished algorithm without converging.")
    return(xvals)
end

function compute_steps(xcurr, f, g, h, opt_type, Bk, xprev)
    """
    Computes the step size and step direction depending on the type of method

    xcurr: current iterate
    f: objective function
    grad_f: gradient of the objective
    hess_f: Hessian of the objective, not required for all methods
    opt_type: type of step direction
    """

    fk = f(xcurr)
    gk = g(xcurr)
    ncalls = 0

    if opt_type == "steepest"
        p = steepest_descent(gk)
        alpha = get_step_size(3, 0.9, xcurr, f, g, p)
    elseif opt_type == "newton"
        Bk = h(xcurr)
        if !check_posdef(Bk)
            ncalls += 1
            Bk = cholesky_mod(1e-3, Bk)
        end
        p = newton(gk, Bk)
        alpha = get_step_size(1, 0.9, xcurr, f, g, p)

    elseif opt_type == "newton_CG"
        Bk = h(xcurr)
        if !check_posdef(Bk)
            ncalls += 1
            Bk = cholesky_mod(1e-3, Bk)
        end
        p = newton_CG(gk, Bk)
        alpha = get_step_size(1, 0.9, xcurr, f, g, p)
    elseif opt_type == "BFGS"
        # here Bk is actually the approx inverse Hessian
        if Bk == nothing
            Hk = h(xcurr)
            if !check_posdef(Hk)
                ncalls += 1
                Hk = cholesky_mod(1e-3, Hk)
            end
            Bk = inv(Hk)
            p = -Bk * gk
        else
        #if xprev != nothing
            sk = xcurr - xprev
            yk = gk - g(xprev)
            Bk = BFGS(sk, yk, Bk)
            p = -Bk * gk
        end

        alpha = get_step_size(1, 0.9, xcurr, f, g, p)
    else
        println(1)
    end

    return p, alpha, ncalls, Bk
end

function get_step_size(alpha_init, rho, xcurr, f, grad_f, p_k, c=1e-4)

    alpha = alpha_init
    fk = f(xcurr)
    gradfk = grad_f(xcurr)

    while f(xcurr + alpha*p_k) > fk + c * alpha * sum(gradfk'*p_k)
        alpha = rho * alpha
    end
    return(alpha)
end

function steepest_descent(gk)
    return(-gk)
end

function newton(gk, Bk)
    return(-Bk \ gk)
end

function newton_CG(gk, Bk)
    gradnorm = norm(gk, 2)[1]
    tol = min(0.5, sqrt(gradnorm)) * gradnorm

    z = zeros(gk); r = gk; d = -gk

    j = 1
    while true
        if (d' * Bk * d)[1] <= 0
            if j == 1
                return -fk
            else
                return z
            end
        end
        a = (r'*r ./ (d' * Bk * d))[1]
        z += a * d
        rold = r
        r += a * Bk * d
        if norm(r, 2) <= tol
            return z
        end
        beta = ((r' * r) ./ (rold' * rold))[1]
        d = -r + beta * d
        j += 1
    end
end

function BFGS(sk, yk, Bk)
"""
    Performs BFGS update
    returns search direction approx **inverse Hessian** Bk

    sk: xnew - xcurr
    yk: gradnew - gradcurr
    Bk: approx to inverse Hessian
"""
    rho = (1 ./ (yk' * sk))[1]
    return (eye(Bk) - rho * sk * yk') * Bk * (eye(Bk) - rho * yk * sk') + rho * sk * sk'
end

function cholesky_mod(beta, H)
    @assert beta > 0
    if minimum(diag(H)) > 0
        tau = 0
    else
        tau = -minimum(diag(H)) + beta
    end
    while true
        candidate = H + tau * eye(H)
        check_posdef(candidate) && return candidate
        tau = max(2*tau, beta)
    end
end

function check_posdef(A)
    try
        L = chol(A)
        return true
    catch
        return false
    end
end


text:
function _pairwise!(::Val{:none}, f, dest::AbstractMatrix, x, y, symmetric::Bool)
    @inbounds for (i, xi) in enumerate(x), (j, yj) in enumerate(y)
        symmetric && i > j && continue

        # For performance, diagonal is special-cased
        if f === cor && eltype(dest) !== Union{} && i == j && xi === yj
            # TODO: float() will not be needed after JuliaLang/Statistics.jl#61
            dest[i, j] = float(cor(xi))
        else
            dest[i, j] = f(xi, yj)
        end
    end
    if symmetric
        m, n = size(dest)
        @inbounds for j in 1:n, i in (j+1):m
            dest[i, j] = dest[j, i]
        end
    end
    return dest
end

function check_vectors(x, y, skipmissing::Symbol)
    m = length(x)
    n = length(y)
    if !(all(xi -> xi isa AbstractVector, x) && all(yi -> yi isa AbstractVector, y))
        throw(ArgumentError("All entries in x and y must be vectors " *
                            "when skipmissing=:$skipmissing"))
    end
    if m > 1
        indsx = keys(first(x))
        for i in 2:m
            keys(x[i]) == indsx ||
                throw(ArgumentError("All input vectors must have the same indices"))
        end
    end
    if n > 1
        indsy = keys(first(y))
        for j in 2:n
            keys(y[j]) == indsy ||
                throw(ArgumentError("All input vectors must have the same indices"))
        end
    end
    if m > 1 && n > 1
        indsx == indsy ||
            throw(ArgumentError("All input vectors must have the same indices"))
    end
end

function _pairwise!(::Val{:pairwise}, f, dest::AbstractMatrix, x, y, symmetric::Bool)
    check_vectors(x, y, :pairwise)
    @inbounds for (j, yj) in enumerate(y)
        ynminds = .!ismissing.(yj)
        @inbounds for (i, xi) in enumerate(x)
            symmetric && i > j && continue

            if xi === yj
                ynm = view(yj, ynminds)
                # For performance, diagonal is special-cased
                if f === cor && eltype(dest) !== Union{} && i == j
                    # TODO: float() will not be needed after JuliaLang/Statistics.jl#61
                    dest[i, j] = float(cor(xi))
                else
                    dest[i, j] = f(ynm, ynm)
                end
            else
                nminds = .!ismissing.(xi) .& ynminds
                xnm = view(xi, nminds)
                ynm = view(yj, nminds)
                dest[i, j] = f(xnm, ynm)
            end
        end
    end
    if symmetric
        m, n = size(dest)
        @inbounds for j in 1:n, i in (j+1):m
            dest[i, j] = dest[j, i]
        end
    end
    return dest
end

function _pairwise!(::Val{:listwise}, f, dest::AbstractMatrix, x, y, symmetric::Bool)
    check_vectors(x, y, :listwise)
    m, n = size(dest)
    nminds = .!ismissing.(first(x))
    @inbounds for xi in Iterators.drop(x, 1)
        nminds .&= .!ismissing.(xi)
    end
    if x !== y
        @inbounds for yj in y
            nminds .&= .!ismissing.(yj)
        end
    end

    # Computing integer indices once for all vectors is faster
    nminds′ = findall(nminds)
    # TODO: check whether wrapping views in a custom array type which asserts
    # that entries cannot be `missing` (similar to `skipmissing`)
    # could offer better performance
    return _pairwise!(Val(:none), f, dest,
                      [view(xi, nminds′) for xi in x],
                      [view(yi, nminds′) for yi in y],
                      symmetric)
end

function _pairwise!(f, dest::AbstractMatrix, x, y;
                    symmetric::Bool=false, skipmissing::Symbol=:none)
    if !(skipmissing in (:none, :pairwise, :listwise))
        throw(ArgumentError("skipmissing must be one of :none, :pairwise or :listwise"))
    end

    x′ = x isa Union{AbstractArray, Tuple, NamedTuple} ? x : collect(x)
    y′ = y isa Union{AbstractArray, Tuple, NamedTuple} ? y : collect(y)
    m = length(x′)
    n = length(y′)

    size(dest) != (m, n) &&
        throw(DimensionMismatch("dest has dimensions $(size(dest)) but expected ($m, $n)"))

    Base.has_offset_axes(dest) && throw("dest indices must start at 1")

    return _pairwise!(Val(skipmissing), f, dest, x′, y′, symmetric)
end

function _pairwise(::Val{skipmissing}, f, x, y, symmetric::Bool) where {skipmissing}
    x′ = x isa Union{AbstractArray, Tuple, NamedTuple} ? x : collect(x)
    y′ = y isa Union{AbstractArray, Tuple, NamedTuple} ? y : collect(y)
    m = length(x′)
    n = length(y′)

    T = Core.Compiler.return_type(f, Tuple{eltype(x′), eltype(y′)})
    Tsm = Core.Compiler.return_type((x, y) -> f(disallowmissing(x), disallowmissing(y)),
                                     Tuple{eltype(x′), eltype(y′)})

    if skipmissing === :none
        dest = Matrix{T}(undef, m, n)
    elseif skipmissing in (:pairwise, :listwise)
        dest = Matrix{Tsm}(undef, m, n)
    else
        throw(ArgumentError("skipmissing must be one of :none, :pairwise or :listwise"))
    end

    # Preserve inferred element type
    isempty(dest) && return dest

    _pairwise!(f, dest, x′, y′, symmetric=symmetric, skipmissing=skipmissing)

    if isconcretetype(eltype(dest))
        return dest
    else
        # Final eltype depends on actual contents (consistent with map and broadcast)
        U = mapreduce(typeof, promote_type, dest)
        # V is inferred (contrary to U), but it only gives an upper bound for U
        V = promote_type(T, Tsm)
        return convert(Matrix{U}, dest)::Matrix{<:V}
    end
end

"""
    pairwise!(f, dest::AbstractMatrix, x[, y];
              symmetric::Bool=false, skipmissing::Symbol=:none)

Store in matrix `dest` the result of applying `f` to all possible pairs
of entries in iterators `x` and `y`, and return it. Rows correspond to
entries in `x` and columns to entries in `y`, and `dest` must therefore
be of size `length(x) × length(y)`.
If `y` is omitted then `x` is crossed with itself.

As a special case, if `f` is `cor`, diagonal cells for which entries
from `x` and `y` are identical (according to `===`) are set to one even
in the presence `missing`, `NaN` or `Inf` entries.

# Keyword arguments
- `symmetric::Bool=false`: If `true`, `f` is only called to compute
  for the lower triangle of the matrix, and these values are copied
  to fill the upper triangle. Only allowed when `y` is omitted.
  Defaults to `true` when `f` is `cor` or `cov`.
- `skipmissing::Symbol=:none`: If `:none` (the default), missing values
  in inputs are passed to `f` without any modification.
  Use `:pairwise` to skip entries with a `missing` value in either
  of the two vectors passed to `f` for a given pair of vectors in `x` and `y`.
  Use `:listwise` to skip entries with a `missing` value in any of the
  vectors in `x` or `y`; note that this might drop a large part of entries.
  Only allowed when entries in `x` and `y` are vectors.

# Examples
```jldoctest
julia> using StatsBase, Statistics

julia> dest = zeros(3, 3);

julia> x = [1 3 7
            2 5 6
            3 8 4
            4 6 2];

julia> pairwise!(cor, dest, eachcol(x));

julia> dest
3×3 Matrix{Float64}:
  1.0        0.744208  -0.989778
  0.744208   1.0       -0.68605
 -0.989778  -0.68605    1.0

julia> y = [1 3 missing
            2 5 6
            3 missing 2
            4 6 2];

julia> pairwise!(cor, dest, eachcol(y), skipmissing=:pairwise);

julia> dest
3×3 Matrix{Float64}:
  1.0        0.928571  -0.866025
  0.928571   1.0       -1.0
 -0.866025  -1.0        1.0
```
"""
function pairwise!(f, dest::AbstractMatrix, x, y=x;
                   symmetric::Bool=false, skipmissing::Symbol=:none)
    if symmetric && x !== y
        throw(ArgumentError("symmetric=true only makes sense passing " *
                            "a single set of variables (x === y)"))
    end

    return _pairwise!(f, dest, x, y, symmetric=symmetric, skipmissing=skipmissing)
end

"""
    pairwise(f, x[, y];
             symmetric::Bool=false, skipmissing::Symbol=:none)

Return a matrix holding the result of applying `f` to all possible pairs
of entries in iterators `x` and `y`. Rows correspond to
entries in `x` and columns to entries in `y`. If `y` is omitted then a
square matrix crossing `x` with itself is returned.

As a special case, if `f` is `cor`, diagonal cells for which entries
from `x` and `y` are identical (according to `===`) are set to one even
in the presence `missing`, `NaN` or `Inf` entries.

# Keyword arguments
- `symmetric::Bool=false`: If `true`, `f` is only called to compute
  for the lower triangle of the matrix, and these values are copied
  to fill the upper triangle. Only allowed when `y` is omitted.
  Defaults to `true` when `f` is `cor` or `cov`.
- `skipmissing::Symbol=:none`: If `:none` (the default), missing values
  in inputs are passed to `f` without any modification.
  Use `:pairwise` to skip entries with a `missing` value in either
  of the two vectors passed to `f` for a given pair of vectors in `x` and `y`.
  Use `:listwise` to skip entries with a `missing` value in any of the
  vectors in `x` or `y`; note that this might drop a large part of entries.
  Only allowed when entries in `x` and `y` are vectors.

# Examples
```jldoctest
julia> using StatsBase, Statistics

julia> x = [1 3 7
            2 5 6
            3 8 4
            4 6 2];

julia> pairwise(cor, eachcol(x))
3×3 Matrix{Float64}:
  1.0        0.744208  -0.989778
  0.744208   1.0       -0.68605
 -0.989778  -0.68605    1.0

julia> y = [1 3 missing
            2 5 6
            3 missing 2
            4 6 2];

julia> pairwise(cor, eachcol(y), skipmissing=:pairwise)
3×3 Matrix{Float64}:
  1.0        0.928571  -0.866025
  0.928571   1.0       -1.0
 -0.866025  -1.0        1.0
```
"""
function pairwise(f, x, y=x; symmetric::Bool=false, skipmissing::Symbol=:none)
    if symmetric && x !== y
        throw(ArgumentError("symmetric=true only makes sense passing " *
                            "a single set of variables (x === y)"))
    end

    return _pairwise(Val(skipmissing), f, x, y, symmetric)
end

# cov(x) is faster than cov(x, x)
_cov(x, y) = x === y ? cov(x) : cov(x, y)

pairwise!(::typeof(cov), dest::AbstractMatrix, x, y;
          symmetric::Bool=false, skipmissing::Symbol=:none) =
    pairwise!(_cov, dest, x, y, symmetric=symmetric, skipmissing=skipmissing)

pairwise(::typeof(cov), x, y; symmetric::Bool=false, skipmissing::Symbol=:none) =
    pairwise(_cov, x, y, symmetric=symmetric, skipmissing=skipmissing)

pairwise!(::typeof(cov), dest::AbstractMatrix, x;
         symmetric::Bool=true, skipmissing::Symbol=:none) =
    pairwise!(_cov, dest, x, x, symmetric=symmetric, skipmissing=skipmissing)

pairwise(::typeof(cov), x; symmetric::Bool=true, skipmissing::Symbol=:none) =
    pairwise(_cov, x, x, symmetric=symmetric, skipmissing=skipmissing)

pairwise!(::typeof(cor), dest::AbstractMatrix, x;
          symmetric::Bool=true, skipmissing::Symbol=:none) =
    pairwise!(cor, dest, x, x, symmetric=symmetric, skipmissing=skipmissing)

pairwise(::typeof(cor), x; symmetric::Bool=true, skipmissing::Symbol=:none) =
    pairwise(cor, x, x, symmetric=symmetric, skipmissing=skipmissing)


text:
using ProjAssistant
@quickactivate "Chemostat_Rath2017"

@time begin

    using ProgressMeter
    using Base.Threads

    import Chemostat
    import Chemostat.MetNets
    import Chemostat.MetLP
    import Chemostat.MetEP
    const Ch = Chemostat

    import Chemostat_Rath2017
    const ChR = Chemostat_Rath2017
    const Rd = ChR.RathData
    const H1 = ChR.Human1
    const HG = H1.HumanGEM
    const AG = H1.HartGEMs

    using Plots
    using Statistics
end

## ---------------------------------------------------------------------
function _find_epout_files(modelid, tissue)
    modelid = string(modelid)
    tissue = string(tissue)
    files = Dict()
    for file in readdir(datdir(AG); join = true)
        heads, params, ext = parse_dfname(file)
        !("maxent_ep_beta0" in heads) && continue
        !("epout" in heads) && continue
        get(params, "tissue", nothing) != tissue && continue
        get(params, "modelid", nothing) != modelid && continue
        it = get(params, "it", nothing)
        isnothing(it) && continue
        files[it] = file
    end
    return files
end

## ---------------------------------------------------------------------
# stoi err
function _plot_stoi_err(modelid, tissue)

    files = _find_epout_files(modelid, tissue)

    model = AG.load_model(;modelid, tissue, uncompress = false)
    M, N = size(model)

    # get fluxs
    cid = (:AVS_CACHE, modelid, tissue)
    # delcache(AG, cid) # Reset
    its, avs = lcache(AG, cid) do
        _avs = Vector{Vector{Float64}}()
        _its = Vector{Int}()
        for (it, file) in sort!(collect(files)) # Test
            _epout = ldat(file)
            push!(_its, it)
            push!(_avs, MetEP.av(_epout))
        end
        return _its, _avs
    end

    # sub range
    itrange = 1:80
    its = its[itrange]
    avs = avs[itrange]

    cid = (:ERRS_CACHE, modelid, tissue, its)
    # delcache(AG, cid) # Reset
    abs_errs, norm_errs = lcache(AG, cid) do
        _abs_errs, _norm_errs = [], []
        @showprogress for (it, av) in zip(its, avs)
            _abs_err = MetNets.stoi_err(model, av)
            _norm_err = MetNets.norm_stoi_err(model, av; normfun = (metv) -> maximum(abs, metv))
            push!(_abs_errs, _abs_err)
            push!(_norm_errs, _norm_err)
        end
        return _abs_errs, _norm_errs
    end

    # get errs
    stoi_errav, stoi_errstd = Float64[], Float64[]
    stoi_norm_errav, stoi_norm_errstd = Float64[], Float64[]
    stoi_norm2_errav, stoi_norm2_errstd = Float64[], Float64[]
    errs_samples = Dict{Int, Vector{Float64}}()
    samples_its = its[1:end ÷ 4:end]
    @show samples_its
    @showprogress for (it, abs_err, norm_err) in zip(its, abs_errs, norm_errs)

        # trend

        push!(stoi_errav, maximum(abs, abs_err))
        push!(stoi_errstd, std(abs_err))

        push!(stoi_norm_errav, maximum(abs, norm_err))
        push!(stoi_norm_errstd, std(norm_err))

        push!(stoi_norm2_errav, mean(abs, norm_err))
        push!(stoi_norm2_errstd, 0.0)

        # hist
        if it in samples_its
            errs_samples[it] = abs.(norm_err)
        end

    end

    # Err trend
    _fun(arr) = log10.(arr)
    p1 = plot(its, _fun(stoi_errav);
        label = "", lw = 3, c = :black, xlabel = "EP it", ylabel = "abs stoi err", title = modelid,
        ylim = [0.0, Inf]
    )
    plot!(p1, its, _fun(stoi_errav .- stoi_errstd); label = "", lw = 1, ls = :dash, c = :black)
    plot!(p1, its, _fun(stoi_errav .+ stoi_errstd); label = "", lw = 1, ls = :dash, c = :black)

    p2 = plot(its, _fun(stoi_norm_errav);
        label = "", lw = 3, c = :black, xlabel = "EP it", ylabel = "maximum norm stoi err", title = modelid,
        ylim = [0.0, Inf]
    )
    plot!(p2, its, _fun(stoi_norm_errav .- stoi_norm_errstd); label = "", lw = 1, ls = :dash, c = :black)
    plot!(p2, its, _fun(stoi_norm_errav .+ stoi_norm_errstd); label = "", lw = 1, ls = :dash, c = :black)

    p3 = plot(its, _fun(stoi_norm2_errav);
        label = "", lw = 3, c = :black, xlabel = "EP it", ylabel = "mean norm stoi err", title = modelid
    )
    plot!(p3, its, _fun(stoi_norm2_errav .- stoi_norm2_errstd); label = "", lw = 1, ls = :dash, c = :black)
    plot!(p3, its, _fun(stoi_norm2_errav .+ stoi_norm2_errstd); label = "", lw = 1, ls = :dash, c = :black)

    ps = [p1, p2, p3]
    sfig(AG, ps,
        @fileid, "maxent_ep_beta0", "stoi_err", (;modelid, tissue), ".png";
        layout = (1, length(ps))
    )

    # Err hist
    ps = Plots.Plot[]
    for (it, errs) in sort!(collect(errs_samples))

        # Show big errs
        if it == last(samples_its)
            _max = maximum(abs, errs)
            _th = 0.1
            _big_erris = findall(errs .> _max * _th)

            println()
            println("-"^54)
            println("-"^54)
            @info("Doing", it, errth =  _max * _th)
            for meti in _big_erris
                println("-"^54)
                MetNets.summary(model, model.mets[meti])
                println()
                println("rel stoi err: ", errs[meti])
                println()
            end
            println()
        end

        p = histogram(_fun(errs);
            label = "", title = string(modelid, " it: ", it),
            bins = 200, xlabel = "stoi err"
        )
        push!(ps, p)
    end

    sfig(AG, ps,
        @fileid, "maxent_ep_beta0", "stoi_err_errs_samples", (;modelid, tissue), ".png";
        layout = (1, length(ps))
    )

    # Err vs met_rxns
    lit = last(samples_its)
    lerrs = errs_samples[lit]
    p = plot(; xlabel = "log rxn count", ylabel = "log stoi err", title = string(modelid, " it: ", lit))
    rxns_counts = Int[]
    _errs = Float64[]
    for (meti, err) in enumerate(lerrs)
        rxnis = MetNets.met_rxns(model, meti)
        rxns_count = length(rxnis)
        rxns_count == 0 && continue
        err == 0 && continue
        if rxns_count > 100
            MetNets.summary(model, model.mets[meti])
        end
        push!(rxns_counts, rxns_count)
        push!(_errs, err)
    end
    scatter!(p, log10.(rxns_counts), log10.(_errs); label = "", c = :black)
    sfig(AG, p,
        @fileid, "maxent_ep_beta0", "stoi_err_vs_rxns_count", (;modelid, tissue, it = lit), ".png"
    )

end

let
    tissue = "GBM"
    for modelid in ["fva_base", "fva_scaled"]
        _plot_stoi_err(modelid, tissue)
    end
end

## ---------------------------------------------------------------------
let
    tissue = "GBM"
    modelid = "base"
    # model = AG.load_model(;modelid, tissue, uncompress = false)
    model = HG.load_humangem_base_model(;uncompress = false)

    big_err_mets = [ "m01401s", "m02394c", "m02394s", "m02680c", "m02680s", "m02871m", "m02877m" ]
    for met in big_err_mets
        println("-"^54)
        MetNets.summary(model, met)
        println()
    end
    # MetNets.balance_str.([model], big_err_mets)
end

## ---------------------------------------------------------------------
# let
#     tissue = "GBM"
#     modelid = "fva_base"
#     model = AG.load_model(modelid, tissue; uncompress = false)
#     MetNets.similar_rxns(model)
# end


text:
using SimpleSolvers
using Test
using GeometricIntegrators
using GeometricIntegrators.Utils
using GeometricProblems.LotkaVolterra3d

SimpleSolvers.set_config(:nls_atol, 8eps())
SimpleSolvers.set_config(:nls_rtol, 2eps())

const Δt = 0.01
const nt = 1000

const ref = [0.39947308320241187, 1.9479527336244262, 2.570183075433086]

@testset "$(rpad("Lotka-Volterra 3D",80))" begin
    ode  = lotka_volterra_3d_ode()

    int = Integrator(ode, TableauGauss(1), Δt)
    sol = integrate(ode, int, nt)
    H, ΔH = compute_energy_error(sol.t, sol.q)
    C, ΔC = compute_casimir_error(sol.t, sol.q)
    @test relative_maximum_error(sol.q, ref) < 5E-4
    @test ΔH[end] < 4E-6
    @test ΔC[end] < 8E-6

    int = Integrator(ode, TableauGauss(2), Δt)
    sol = integrate(ode, int, nt)
    H, ΔH = compute_energy_error(sol.t, sol.q)
    C, ΔC = compute_casimir_error(sol.t, sol.q)
    @test relative_maximum_error(sol.q, ref) < 2E-9
    @test ΔH[end] < 5E-11
    @test ΔC[end] < 2E-11

end


text:
abstract type AbstractObjective end

clear!(::Function) = nothing

mutable struct UnivariateObjective{TF, TD, Tf, Td, Tx} <: AbstractObjective
    F::TF     # objective
    D::TD     # derivative of objective

    f::Tf     # cache for f output
    d::Td     # cache for d output

    x_f::Tx   # x used to evaluate F (stored in f)
    x_d::Tx   # x used to evaluate D (stored in d)

    f_calls::Int
    d_calls::Int
end

function UnivariateObjective(F::Callable, D::Callable, x::Number;
                             f::Real = alloc_f(x),
                             d::Number = alloc_d(x))
    UnivariateObjective(F, D, f, d, alloc_x(x), alloc_x(x), 0, 0)
end

function UnivariateObjective(F::Callable, x::Number; kwargs...)
    D = (x) -> ForwardDiff.derivative(F,x)
    UnivariateObjective(F, D, x; kwargs...)
end

UnivariateObjective(F::Callable, D::Nothing, x::Number; kwargs...) = UnivariateObjective(F, x; kwargs...)

"""
Evaluates the objective value at `x`.
Returns `f(x)`, but does *not* store the value in `obj.F`
"""
function value(obj::UnivariateObjective, x)
    obj.f_calls += 1
    return obj.F(x)
end

"Get the most recently evaluated objective value of `obj`."
value(obj::UnivariateObjective) = obj.f

"""
Force (re-)evaluation of the objective value at `x`.
Returns `f(x)` and stores the value in `obj.F`
"""
function value!!(obj::UnivariateObjective, x)
    obj.x_f = x
    obj.f = obj.F(x)
    obj.f_calls += 1
    value(obj)
end

"""
Evaluates the objective value at `x`.
Returns `f(x)` and stores the value in `obj.F`
"""
function value!(obj::UnivariateObjective, x)
    if x != obj.x_f
        value!!(obj, x)
    end
    value(obj)
end

"""
Evaluates the derivative of the objective at `x`.
Returns `f'(x)`, but does *not* store the derivative in `obj.D`
"""
function derivative(obj::UnivariateObjective, x)
    obj.d_calls += 1
    return obj.D(x)
end

"Get the most recently evaluated derivative of the objective of `obj`."
derivative(obj::UnivariateObjective) = obj.d

"""
Force (re-)evaluation of the derivative of the objective at `x`.
Returns `f'(x)` and stores the derivative in `obj.D`
"""
function derivative!!(obj::UnivariateObjective, x)
    obj.x_d = x
    obj.d = obj.D(x)
    obj.d_calls += 1
    derivative(obj)
end

"""
Evaluates the derivative of the objective at `x`.
Returns `f'(x)` and stores the derivative in `obj.D`
"""
function derivative!(obj::UnivariateObjective, x)
    if x != obj.x_d
        derivative!!(obj, x)
    end
    derivative(obj)
end

(obj::UnivariateObjective)(x) = value!(obj, x)

function _clear_f!(obj::UnivariateObjective)
    obj.f_calls = 0
    obj.f = typeof(obj.f)(NaN)
    obj.x_f = typeof(obj.x_f)(NaN)
    nothing
end

function _clear_d!(obj::UnivariateObjective)
    obj.d_calls = 0
    obj.d = eltype(obj.d)(NaN)
    obj.x_d = eltype(obj.x_d)(NaN)
    nothing
end

function clear!(obj::UnivariateObjective)
    _clear_f!(obj)
    _clear_d!(obj)
    nothing
end

mutable struct MultivariateObjective{TF, TG, Tf, Tg, Tx} <: AbstractObjective
    F::TF
    G::TG

    f::Tf
    g::Tg

    x_f::Tx
    x_g::Tx

    f_calls::Int
    g_calls::Int
end

function MultivariateObjective(F, G,
                               x::AbstractArray;
                               f::Real = alloc_f(x),
                               g::AbstractArray = alloc_g(x))
    MultivariateObjective(F, G, f, g, alloc_x(x), alloc_x(x), 0, 0)
end

function MultivariateObjective(F, x::AbstractArray; kwargs...)
    G = GradientParametersAD(F, x)
    MultivariateObjective(F, G, x; kwargs...)
end

MultivariateObjective(F, G::Nothing, x::AbstractArray; kwargs...) = MultivariateObjective(F, x; kwargs...)

"""
Evaluates the objective value at `x`.
Returns `f(x)`, but does *not* store the value in `obj.f`
"""
function value(obj::MultivariateObjective, x)
    obj.f_calls += 1
    return obj.F(x)
end

"Get the most recently evaluated objective value of `obj`."
value(obj::MultivariateObjective) = obj.f

"""
Force (re-)evaluation of the objective at `x`.
Returns `f(x)` and stores the value in `obj.f`
"""
function value!!(obj::MultivariateObjective, x)
    copyto!(obj.x_f, x)
    obj.f = obj.F(x)
    obj.f_calls += 1
    value(obj)
end

"""
Evaluates the objective at `x`.
Returns `f(x)` and stores the value in `obj.f`
"""
function value!(obj::MultivariateObjective, x)
    if x != obj.x_f
        value!!(obj, x)
    end
    value(obj)
end

"Get the most recently evaluated gradient of `obj`."
gradient(obj::MultivariateObjective) = obj.g

"""
Evaluates the gradient at `x`.
This does *not* update `obj.g` or `obj.x_g`.
"""
function gradient(obj::MultivariateObjective, x)
    ḡ = copy(obj.g)
    obj.G(ḡ, x)
    obj.g_calls += 1
    return ḡ
end

"""
Force (re-)evaluation of the gradient at `x`.
Returns ∇f(x) and stores the value in `obj.g`.
"""
function gradient!!(obj::MultivariateObjective, x)
    copyto!(obj.x_g, x)
    obj.G(obj.g, x)
    obj.g_calls += 1
    gradient(obj)
end

"""
Evaluates the gradient at `x`.
Returns ∇f(x) and stores the value in `obj.g`.
"""
function gradient!(obj::MultivariateObjective, x)
    if x != obj.x_g
        gradient!!(obj, x)
    end
    gradient(obj)
end

(obj::MultivariateObjective)(x) = value!(obj, x)

function _clear_f!(obj::MultivariateObjective)
    obj.f_calls = 0
    obj.f = typeof(obj.f)(NaN)
    obj.x_f .= eltype(obj.x_f)(NaN)
    nothing
end

function _clear_g!(obj::MultivariateObjective)
    obj.g_calls = 0
    obj.g .= eltype(obj.g)(NaN)
    obj.x_g .= eltype(obj.x_g)(NaN)
    nothing
end

function clear!(obj::MultivariateObjective)
    _clear_f!(obj)
    _clear_g!(obj)
    nothing
end

f_calls(o::AbstractObjective) = error("f_calls is not implemented for $(summary(o)).")
f_calls(o::UnivariateObjective) = o.f_calls
f_calls(o::MultivariateObjective) = o.f_calls

d_calls(o::AbstractObjective) = error("d_calls is not implemented for $(summary(o)).")
d_calls(o::UnivariateObjective) = o.d_calls

g_calls(o::AbstractObjective) = error("g_calls is not implemented for $(summary(o)).")
g_calls(o::MultivariateObjective) = o.g_calls


text:
struct MachOStrTab{H <: MachOHandle} <: StrTab{H}
    cmd::MachOLoadCmdRef{H,MachOSymtabCmd{H}}
end

function StrTab(cmd::MachOLoadCmdRef{H,MachOSymtabCmd{H}}) where {H <: MachOHandle}
    return MachOStrTab(cmd)
end
function StrTab(oh::MachOHandle)
    cmds = findall(MachOLoadCmds(oh), [MachOSymtabCmd])
    if isempty(cmds)
        error("Mach-O file does not contain Symtab load commands")
    end
    return StrTab(first(cmds))
end

handle(strtab::MachOStrTab) = handle(strtab.cmd)
strtab_lookup(strtab::MachOStrTab, index) = strtab_lookup(strtab.cmd, index)


text:
using Tk
using Winston

win = Toplevel("TestCanvas", 400, 200)
c = Canvas(win, 400, 200)
pack(c, expand=true, fill="both")

p = FramedPlot()
x = linspace(0,10,1001)
y = sin(x)
add(p, Curve(x, y, "color", "green"))

Winston.display(c, p)


text:
"""
Gallery of RTG models
"""

abstract type RTGmodel end

include("utils.jl")
include("actM1.jl")
include("rtgM1.jl")
include("rtgM2.jl")
include("rtgM3.jl")
include("rtgM4.jl")