text:
using Revise
using AdFem
using PyCall
using LinearAlgebra
using ADCME
using MAT
using JLD2
using PyPlot
np = pyimport("numpy")

reset_default_graph()
function nnlaw(σ, ε)
    x = [σ ε]
    ae(x, [20,20,20,20,3], "law")
end

# end
mode = "data"
# Domain information
NT = 20
Δt = 1/NT
n = 10
m = 2*n
h = 1.0/n
bdnode = bcnode("lower", m, n, h)

b = 1.0
H_ = Variable(diagm(0=>ones(3)))
H = H_'*H_
H = H .* [1.0 1.0 0.0
                    1.0 1.0 0.0
                0.0 0.0 1.0]

bd = bcedge("upper", m, n, h)
Q, Prhs = compute_fvm_tpfa_matrix(ones(4*m*n), bd, zeros(size(bd,1)),m, n, h)
Q = SparseTensor(Q)
K = compute_fem_stiffness_matrix(H, m, n, h)
L = SparseTensor(compute_interaction_matrix(m, n, h))
M = SparseTensor(compute_fvm_mass_matrix(m, n, h))
A = [K -b*L'
b*L/Δt 1/Δt*M-Q]
A, Abd = fem_impose_coupled_Dirichlet_boundary_condition(A, bdnode, m, n, h)

U = zeros(m*n+2(m+1)*(n+1), NT+1)
x = Float64[]; y = Float64[]
for j = 1:n+1
    for i = 1:m+1
        push!(x, (i-1)*h)
        push!(y, (j-1)*h)
    end
end
injection = (div(n,2)-1)*m + 3
production = (div(n,2)-1)*m + m-3

function get_disp(ipval)

    function condition(i, tas...)
        i<=NT
    end

    function body(i, tas...)
        ta_u, ta_ε, ta_σ = tas
        u = read(ta_u, i)
        σ0 = read(ta_σ, i)
        ε0 = read(ta_ε, i)
        G = nnlaw(σ0, ε0)
        rhs1 = compute_strain_energy_term(G, m, n, h)

        rhs1 = scatter_update(rhs1, [bdnode; bdnode .+ (m+1)*(n+1)], zeros(2length(bdnode)))
        rhs2 = zeros(m*n)
        rhs2[injection] += 1.0
        rhs2[production] -= 1.0
        rhs2 = rhs2 * ipval + b*L*u[1:2(m+1)*(n+1)]/Δt +
                M * u[2(m+1)*(n+1)+1:end]/Δt + Prhs

        rhs = [rhs1;rhs2]
        o = A\rhs

        ε = eval_strain_on_gauss_pts(o, m, n, h)
        σ = G + ε*H

        ta_u = write(ta_u, i+1, o)
        ta_ε = write(ta_ε, i+1, ε)
        ta_σ = write(ta_σ, i+1, σ)
        i+1, ta_u, ta_ε, ta_σ
    end

    i = constant(1, dtype=Int32)
    ta_u = TensorArray(NT+1); ta_u = write(ta_u, 1, constant(zeros(2(m+1)*(n+1)+m*n)))
    ta_ε = TensorArray(NT+1); ta_ε = write(ta_ε, 1, constant(zeros(4*m*n, 3)))
    ta_σ = TensorArray(NT+1); ta_σ = write(ta_σ, 1, constant(zeros(4*m*n, 3)))
    _, u_out, ε_out, σ_out = while_loop(condition, body, [i, ta_u, ta_ε, ta_σ])
    u_out = stack(u_out)
    u_out.set_shape((NT+1, size(u_out,2)))
    σ_out = stack(σ_out)
    ε_out = stack(ε_out)

    upper_idx = Int64[]
    for i = 1:m+1
        push!(upper_idx, (div(n,3)-1)*(m+1)+i)
        push!(upper_idx, (div(n,3)-1)*(m+1)+i + (m+1)*(n+1))
    end
    for i = 1:m
        push!(upper_idx, (div(n,3)-1)*m+i+2(m+1)*(n+1))
    end

    u_out, σ_out
end

Us = Array{PyObject}(undef, 6)
Ss = Array{PyObject}(undef, 6)

for i = 1:5
    Us[i], Ss[i] = get_disp(0.2*i)
end
Us[6], Ss[6] = get_disp(0.5)

if isfile("invdata.mat")
    d = matread("invdata.mat")
    Us_ = d["U"]
    Ss_ = d["S"]
    global Sigma0 = Ss_[6]
    global U0 = Us_[6]
    global loss = sum([sum((Us[i][:,1:m+1] - Us_[i][:,1:m+1])^2) for i = 1:5])
end
sess = Session(); init(sess)

function visualize(i)
    visualize_von_mises_stress(Sigma_, m, n, h, name="_nn$i")
    visualize_displacement(Array(U_'), m, n, h, name="_nn0$i",
                    xlim_=[-3h, m*h+2h], ylim_=[-2h, n*h+2h])
    close("all")
    figure(figsize=(13,4))
    subplot(121)
    plot(LinRange(0, 20, NT+1), U0[:,1], "r--")
    plot(LinRange(0, 20, NT+1), U_[:,1], "r")
    plot(LinRange(0, 20, NT+1), U0[:,1+(n+1)*(m+1)], "g--")
    plot(LinRange(0, 20, NT+1), U_[:,1+(n+1)*(m+1)], "g")
    xlabel("Time")
    ylabel("Displacement")
    subplot(122)
    plot(LinRange(0, 20, NT+1), mean(Sigma0[:,1:4,1], dims=2)[:],"r--", label="\$\\sigma_{xx}\$")
    plot(LinRange(0, 20, NT+1), mean(Sigma0[:,1:4,2], dims=2)[:],"b--", label="\$\\sigma_{yy}\$")
    plot(LinRange(0, 20, NT+1), mean(Sigma0[:,1:4,3], dims=2)[:],"g--", label="\$\\sigma_{xy}\$")
    plot(LinRange(0, 20, NT+1), mean(Sigma_[:,1:4,1], dims=2)[:],"r-")
    plot(LinRange(0, 20, NT+1), mean(Sigma_[:,1:4,2], dims=2)[:],"b-")
    plot(LinRange(0, 20, NT+1), mean(Sigma_[:,1:4,3], dims=2)[:],"g-")
    legend()
    legend()
    xlabel("Time")
    ylabel("Stress")
    savefig("disp$i.png")
    savefig("disp$i.pdf")
    matwrite("nn$i.mat", Dict("U"=>U_, "S"=>Sigma_))
end

# run(sess, loss)
# gradview(sess, H_, loss, diagm(0=>ones(3))[:])
# savefig("line.png")

loss = loss*1e10
U_, Sigma_ = run(sess, [Us[6], Ss[6]])
visualize(0)
for iter = 1:5000
    BFGS!(sess, loss, 500)
    global U_, Sigma_ = run(sess, [Us[6], Ss[6]])
    visualize(iter)
end

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
using Test, OrderedCollections
include(joinpath(dirname(@__FILE__), "../../src/includeall.jl"))

## Check direct implementation of functional iteration loop matches output from solve
# Initialize
m = Li2020()
m <= Setting(:tol, 1e-2)
m <= Setting(:learning_rate, 0.6) # Set to 0.6 to speed it up relative to the Li 2020 solution
stategrid, funcvar, derivs, endo = initialize!(m)
funcvar_dict = get_functional_variables(m)
func_iter = eqcond(m)
θ = parameters_to_named_tuple(get_parameters(m))

# Get settings for the loop
max_iter     = get_setting(m, :max_iter)        # maximpum number of functional loops
verbose      = :high
func_tol     = get_setting(m, :tol)             # convergence tolerance
error_method = get_setting(m, :error_method)    # method for calculating error at end of each loop
func_errors  = Vector{Float64}(undef, max_iter) # track error for each iteration
proposal_funcvar = deepcopy(funcvar)            # so we can calculate difference between proposal and current values

# Set update method
um = (new, old) -> average_update(new, old, get_setting(m, :learning_rate))

# Start loop
total_time = 0.
error_vars = Vector{Symbol}(undef, 2)
error_vars[1] = :p
error_vars[2] = :Q̂

for iter in 1:max_iter
    begin_time = time_ns()

    # Get a new guess
    new_funcvar = func_iter(stategrid, funcvar, derivs, endo, θ; verbose = verbose)

    # Update the guess
    for (k, v) in proposal_funcvar
        proposal_funcvar[k] = um(new_funcvar[funcvar_dict[k]], v) # um = update_method
    end

    # Calculate errors
    func_errors[iter] = calculate_func_error(proposal_funcvar, funcvar, error_method;
                                             vars = error_vars)

    spaces1, spaces2, spaces3 = if iter < 10
        repeat(" ", 16), repeat(" ", 14), repeat(" ", 5)
    elseif iter < 100
        repeat(" ", 17), repeat(" ", 15), repeat(" ", 6)
    else
        repeat(" ", 18), repeat(" ", 16), repeat(" ", 7)
    end

    println("Iteration $(iter), current error:            $(func_errors[iter])")
    if verbose == :high
        for k in error_vars
            indiv_funcvar_err = calculate_func_error(proposal_funcvar[k], funcvar[k], error_method)
            indiv_space = " " ^ (28 - length(string(k)) + 1)
            println("Error for $(k):" * indiv_space * string(indiv_funcvar_err))
        end
    else
        println("")
    end

    loop_time = (time_ns() - begin_time) / 6e10 # 60 * 1e9 = 6e10
    global total_time += loop_time
    expected_time_remaining = (max_iter - iter) * loop_time
    println(verbose, :high, "Duration of loop (min):" * spaces1 * "$(round(loop_time, digits = 4))")
    println("Total elapsed time (min):" * spaces2 * "$(round(total_time, digits = 4))")
    println(verbose, :high, "Expected max remaining time (min):" * spaces3 * "$(round(expected_time_remaining, digits = 4))")
    println("\n")

    # Convergence?
    for (k, v) in proposal_funcvar
        funcvar[k] .= v
    end

    if func_errors[iter] < func_tol
        if verbose != :none
            println("Convergence achieved! Final round error: $(func_errors[iter])")
        end
        break
    end
end

println("Calculating remaining variables . . .")
aug_time = time_ns()
augment_variables!(m, stategrid, funcvar, derivs, endo)
total_time += (time_ns() - aug_time) / 6e10
println("Total elapsed time (min): $(round(total_time, digits = 4))\n")

stategrid_solve, funcvar_solve, derivs_solve, endo_solve = solve(m; verbose = verbose, vars_for_error = error_vars)

@test stategrid[:w] == stategrid_solve[:w]
for k in keys(funcvar_solve)
    @test maximum(abs.(funcvar[k][1:end - 1] - funcvar_solve[k][1:end - 1])) < 1e-11
    @test abs(funcvar[k][end] - funcvar_solve[k][end]) < 5e-6
end
for k in keys(derivs_solve)
     @test @test_matrix_approx_eq derivs[k] derivs_solve[k]
end
for k in keys(endo_solve)
    @test @test_matrix_approx_eq endo[k][1:end - 1] endo_solve[k][1:end - 1]
end

ext: jl


text:
# Create example data
x = LinRange(-2, 2, 40)
f(t) = t^3 + t^2 + t
y = f.(x)
# Plot x and y
plot(x, y)
# Plot x using the callable
plot(x, f)
# Plot y, using its indices for the x values
plot(y)
# Plot two columns
plot(x, [y 3x.^2 .- 3])
# Plot various series them with custom line specs
y2 = 3x.^2 .- 3
plot(x, y, "-r", x, y2, ":*b")


text:
struct NDMSymmBatchedCache{BC,BD} <: NNBatchedCache{NDMSymm}
    bare_cache::BC
    bare_∇logψ::BD
end

cache(net::NDMSymm, batch_sz) =
    NDMSymmBatchedCache(cache(net.bare_net, batch_sz),
                        grad_cache(net.bare_net, batch_sz))

batch_size(c::NDMSymmBatchedCache) = batch_size(c.bare_cache)

logψ!(out::AbstractMatrix, net::NDMSymm, c::NDMSymmBatchedCache, σr_r::AbstractMatrix, σc_r::AbstractMatrix) =
    logψ!(out, net.bare_net, c.bare_cache, σr_r, σc_r)

function logψ_and_∇logψ!(∇logψ, out::AbstractMatrix, net::NDMSymm, c::NDMSymmBatchedCache, σr_r::AbstractMatrix, σc_r::AbstractMatrix)
    out = logψ_and_∇logψ!(c.bare_∇logψ, out, net.bare_net, c.bare_cache, σr_r, σc_r)

    symmetrize_∇logψ_NDM_batched!(∇logψ, c.bare_∇logψ, net)
    return out
end

function symmetrize_∇logψ_NDM_batched!(∇lnψ_symm, ∇lnψ, net)
    mul!(∇lnψ_symm.b_μ, net.∇b_mat, ∇lnψ.b_μ)
    mul!(∇lnψ_symm.b_λ, net.∇b_mat, ∇lnψ.b_λ)
    mul!(∇lnψ_symm.h_μ, net.∇h_mat, ∇lnψ.h_μ)
    mul!(∇lnψ_symm.h_λ, net.∇h_mat, ∇lnψ.h_λ)
    mul!(∇lnψ_symm.d_λ, net.∇d_mat, ∇lnψ.d_λ)
    for i=1:size(∇lnψ_symm.w_μ, 3)
        mul!(_uview_vec_destride(∇lnψ_symm.w_μ, i), net.∇w_mat, _uview_vec_destride(∇lnψ.w_μ, i))
        mul!(_uview_vec_destride(∇lnψ_symm.w_λ, i), net.∇w_mat, _uview_vec_destride(∇lnψ.w_λ, i))
        mul!(_uview_vec_destride(∇lnψ_symm.u_μ, i), net.∇u_mat, _uview_vec_destride(∇lnψ.u_μ, i))
        mul!(_uview_vec_destride(∇lnψ_symm.u_λ, i), net.∇u_mat, _uview_vec_destride(∇lnψ.u_λ, i))
    end

    return ∇lnψ_symm
end


text:
using ExplainableAI
using Flux
using JLD2

const GRADIENT_ANALYZERS = Dict(
    "Gradient" => Gradient,
    "InputTimesGradient" => InputTimesGradient,
    "SmoothGrad" => model -> SmoothGrad(model, 5, 0.1, MersenneTwister(123)),
)
const LRP_ANALYZERS = Dict(
    "LRPZero" => LRPZero, "LRPEpsilon" => LRPEpsilon, "LRPGamma" => LRPGamma
)

using Random
pseudorand(T, dims...) = rand(MersenneTwister(123), T, dims...)
input_size = (224, 224, 3, 1)
img = pseudorand(Float32, input_size)

# Load VGG model:
# We run the reference test on the randomly intialized weights
# so we don't have to download ~550 MB on every CI run.
include("./vgg11.jl")
vgg11 = VGG11(; pretrain=false)
model = flatten_model(strip_softmax(vgg11.layers))

function LRPCustom(model::Chain)
    return LRP(model, [ZBoxRule(), repeat([GammaRule()], length(model.layers) - 1)...])
end

function test_vgg11(name, method; kwargs...)
    @testset "$name" begin
        # Reference test attribution
        analyzer = method(model)
        print("Timing $name cold...\t")
        @time expl = analyze(img, analyzer; kwargs...)
        attr = expl.attribution
        @test size(attr) == size(img)
        @test_reference "references/vgg11/$(name).jld2" Dict("expl" => attr) by =
            (r, a) -> isapprox(r["expl"], a["expl"]; rtol=0.05)

        # Test direct call of analyzer
        analyzer = method(model)
        print("Timing $name warm...\t")
        @time expl2 = analyzer(img; kwargs...)
        @test expl.attribution ≈ expl2.attribution

        # Test direct call of heatmap
        h1 = heatmap(expl)
        analyzer = method(model)
        h2 = heatmap(img, analyzer; kwargs...)
        @test h1 ≈ h2
        if expl.analyzer != :Gradient # TODO: remove
            @test_reference "references/heatmaps/vgg11_$(name).txt" h1
        end
    end
    @testset "$name neuron selection" begin
        analyzer = method(model)
        neuron_selection = 1
        expl = analyze(img, analyzer, neuron_selection; kwargs...)
        attr = expl.attribution

        @test size(attr) == size(img)
        @test_reference "references/vgg11/$(name)_neuron_$neuron_selection.jld2" Dict(
            "expl" => attr
        ) by = (r, a) -> isapprox(r["expl"], a["expl"]; rtol=0.05)

        analyzer = method(model)
        expl2 = analyzer(img, neuron_selection; kwargs...)
        @test expl.attribution ≈ expl2.attribution
    end
    return nothing
end

# Run analyzers
@testset "LRP analyzers" begin
    for (name, method) in LRP_ANALYZERS
        test_vgg11(name, method)
    end
end
@testset "Custom LRP composite" begin
    test_vgg11("LRPCustom", LRPCustom)
end

@testset "Gradient analyzers" begin
    for (name, method) in GRADIENT_ANALYZERS
        test_vgg11(name, method)
    end
end
# Layerwise relevances in LRP methods
@testset "Layerwise relevances" begin
    test_vgg11("LRPZero", LRPZero; layerwise_relevances=true)
end

# Test LRP constructor with no rules
a1 = LRP(model)
a2 = LRPZero(model)
@test a1.model == a2.model
@test a1.rules == a2.rules


text:
module ParticleFilters

using Distributions
using StatsBase

# import PyPlot.plot
# using PyPlot

export ParticleFilterEstimator, update!

"""
X is an MxN array representing N particles, each with M state variables
"""
type ParticleFilterEstimator
    X::Array{Float64, 2}
    ranges::Array{Float64, 2}
    anchors::Array{Float64, 2}

    function ParticleFilterEstimator(n::Integer, xrange, yrange, zrange, anchors)
        xdiff = maximum(xrange) - minimum(xrange)
        ydiff = maximum(yrange) - minimum(yrange)
        zdiff = maximum(zrange) - minimum(zrange)

        xoff = minimum(xrange)
        yoff = minimum(yrange)
        zoff = minimum(zrange)
        particlesInit = rand(3, n) .* [xdiff, ydiff, zdiff] .- [xoff, yoff, zoff]
        ranges = [xrange[1] yrange[1] zrange[1];
                  xrange[2] yrange[2] zrange[2]]
        new(particlesInit, ranges, anchors)
    end

    function ParticleFilterEstimator(n::Integer=1000)
        # anchor locations for Fluid optitrack test
        anchors = [ -1.1082  2.0964  1.9340 -0.7987;
                  2.3606  2.3635  2.3635  2.3606;
                  4.0778  3.8079 -1.8402 -1.7995]

        ParticleFilterEstimator(n, (-3, 3), (0, 3), (-3, 3), anchors)
    end
end

# function plot(pf::ParticleFilterEstimator, z=nothing)
#     if z == nothing
#         c = nothing
#     else
#         l(x) = likelihood(x, z)
#         c = vec(mapslices(l, pf.X, 1))
#     end
#
#     if size(pf.X, 1) == 2
#         # assume the first 2 state variables are x, y
#         if c != nothing
#             scatter(pf.X[1, :], pf.X[2, :], c=c, s=5)
#             colorbar()
#         else
#             scatter(pf.X[1, :], pf.X[2, :], s=5)
#         end
#     elseif size(pf.X, 1) > 2
#         # assume the first 3 state variables are x, y, z
#         if c != nothing
#             scatter3D(pf.X[1, :]', pf.X[2, :]', pf.X[3, :]', c=c, s=5)
#         else
#             scatter3D(pf.X[1, :]', pf.X[2, :]', pf.X[3, :]', s=5)
#         end
#     else
#         println("Must have at least 2 dimensions to plot")
#     end
# end
#
# """
# Plots the likelihood with the given observation
# """
# function plotlikelihood(pf::ParticleFilterEstimator, z::Vector)
#     v = linspace(0, 1, 100)
#     l = Array(Float64, 100, 100)
#     for (i, x) in enumerate(v)
#         for (j, y) in enumerate(v)
#             l[i, j] = likelihood([x, y], z)
#         end
#     end
#
#     imshow(l; origin="lower")
#     colorbar()
# end

function update!(pf::ParticleFilterEstimator, u, z, maxiters=20, thresh=[0.01, 0.02, 0.01])
    errbias = [0.1153, 0.0571, 0.0237, 0.0538]

    # copy so we don't mess it up when we subtract off the error bias
    z = copy(z)

    for i in 1:length(z)
        if z[i] != -1.0
            # correct for the bias that we measured per-anchor
            z[i] = z[i] - errbias[i]
        end
    end

    fitnorm = fit(DiagNormal, pf.X)
    iters = 0
    while iters < maxiters
        X2 = predict(pf.X, u)
        w = weights(X2, z, pf)
        pf.X, _ = resample(X2, w)
        fitnorm = fit(DiagNormal, pf.X)
        iters += 1
        # check at the end so we make sure we go at least one time through even
        # if the system started out with small covariance
        if all(diag(cov(fitnorm)) .< thresh)
            break
        end
    end

    fitnorm, iters
end

"""
Given the previous state X and control input vector u, generate a prediction for
the current state X'
"""
function predict(X, u)
    # TODO: calculate the covariance matrix from experimental data
    d = MvNormal(zeros(size(X, 1)), ones(size(X, 1)) * 0.1)

    rand(d, size(X, 2)) + X
end

"""
Computes the likelihood of the obvservation z given the state x
"""
function likelihood(x::Vector{Float64}, z::Vector{Float64}, pf::ParticleFilterEstimator)
    # approximate std dev (rounded up a bit) from measured data
    d = Normal(0.0, 0.035)

    l = 1.0
    for i in 1:length(z)
        if z[i] != -1.0
            l *= (0.2 + 0.8 * pdf(d, norm(x - pf.anchors[:, i]) - z[i]))
        else
            l *= 0.2
        end
    end

    for i in 1:3
        if !(pf.ranges[1, i] < x[i] < pf.ranges[2, i])
            # keep our solutions inside the box. This is particularly an issue
            # for reflections across the plane of the sensors, but is also to keep
            # the particle population from drifting out too much if there's no data
            l = 0.0
        end
    end

    l
end

"""
Given the current state X and measurement vector z, generate a weight vector for
the state estimates
"""
function weights(X::Array{Float64, 2}, z::Vector{Float64}, pf::ParticleFilterEstimator)
    # for i in 1:length(z)
    #     # get the difference vector from each location to anchor i
    #     diff = X .- anchors[:, i]
    #     # take the norm of the difference vector (the distance to the anchor)
    #     dist = vec(mapslices(norm, diff, 1))
    #     # find the error (difference between the expected distance and measured)
    #     err = dist - z[i]
    #     # compute the product of the distance likelihoods from each anchor
    #     w = w .* pdf(d, err)
    # end

    # create a temporary function to compute the likelihood with the given z
    l(x::Vector{Float64}) = likelihood(x, z, pf)

    w = vec(mapslices(l, X, 1))
    # now normalize the weights
    w / sum(w)
end

"""
Given the current state estimate X and weight vector w, resample a new set of
states. We use the low-variance resampling algorithm from Thrun, Burgard, and Fox's
"Probabilistic Robotics". By default it will keep the number of samples constant
"""
function resample(X, w, M=size(X, 2))
    X2 = similar(X, (size(X, 1), M))
    w2 = zeros(M)
    r = rand(Uniform(0, 1/M))
    c = w[1]
    i = 1
    for m in 1:M
        # calculate the next sample point
        U = r + (m - 1) * (1/M)
        while c < U
        # find the first weight that puts us past the sample point
            i += 1
            c = c + w[i]
        end
        X2[:, m] = X[:, i]
        w2[m] = w[i]
    end
    X2, w2
end

"""
Resample using the sample function from StatsBase
"""
function resample2(X, w, M=size(X, 2))
    X2 = similar(X, (size(X, 1), M))
    w2 = zeros(M)
    wv = WeightVec(w)

    for m in 1:M
        i = sample(wv)
        X2[:, m] = X[:, i]
        w2[m] = w[i]
    end
    X2, w2
end

end # module


text:
using DWARF
using Base.Test

# Line Table Tests

x = IOBuffer()
# stub (0x2a = 42 bytes after header field, adjust that number if changing the header)
write(x,[0x47,0x00,0x00,0x00,0x02,0x00,0x2a,0x00,0x00,0x00,0x01,0x01,0xfb,0x0e,0x0d]) # 8 bytes after header length
# opcode lengths
standard_opcode_lengths = [0x00,0x01,0x01,0x01,0x01,0x00,0x00,0x00,0x01,0x00,0x00,0x01]
write(x,standard_opcode_lengths) # 12 bytes
# include directories (11 bytes)
write(x,"/julia/test")
write(x,"\0")       # null terminate
write(x,'\0')       # end of list
# Write file entries (12 bytes)
write(x,"test.jl")
write(x,'\0')       # null terminate
write(x,ULEB128(1)) # directory index
write(x,ULEB128(0)) # timestamp
write(x,ULEB128(0)) # file size
write(x,'\0')       # end of list

# The actual line number program
write(x,[0x05,0x01,0x00,0x09,0x02,0x0e,0x30,0xec,0x2c,0x01,0x00,0x00,0x00,0x03,0x14,0x01,0x02,0xac,0x01,0x00,0x01,0x01])

seekstart(x)
header = DWARF.LineTableSupport.read_header(x)
@test header.stub == DWARF.LineTableSupport.HeaderStub{UInt32}(
    0x47, # Total Length
    0x02, # Version
    0x2a, # Header Length (bytes after this point)
    0x01, # minimum_instruction_length
    0x01, # maximum_operations_per_instruction (implicit it DWARF version 2)
    0x01, # default_is_stmt
    reinterpret(Int8, 0xfb), # line_base
    0x0e, # line_range
    0x0d) # opcode_base

@test header.standard_opcode_lengths == standard_opcode_lengths
@test header.include_directories == ["/julia/test"]
@test length(header.file_names) == 1
@test header.file_names[1] == DWARF.LineTableSupport.FileEntry("test.jl",1,0,0)
@test position(x) == header.stub.header_length+10

seekstart(x)
t = LineTable(x)
s = start(t)
r,s = next(t,s)
@test r == DWARF.LineTableSupport.RegisterState(5048643598,0,1,21,1,true,false,false,false,false,0,0)
r,s = next(t,s)
@test r == DWARF.LineTableSupport.RegisterState(5048643770,0,1,21,1,true,false,true,false,false,0,0)
@test done(t,s)


text:
include(joinpath("..","..","src","MaxwellSolver","PoissonFourier1D.jl"))
import .PoissonFourier1D

@testset "Test 1D Poisson solvers" begin
    Δx = 0.01
    ρ = [sin(2π*i*Δx) for i = 0:1000]
    e = similar(ρ)
    solve! = PoissonFourier1D.generate_solver(ρ, Δx)
    solve!(e, ρ)
    @test sum(@. (e - [-cos(2π*i*Δx)/2π for i = 0:1000])^2) ≈ 0 atol=1e-3
end


text:
@testset "nearest_rotation" begin
    @testset "$(N)-dim" for N in [2,3,4]
        for _ in 1:100
            M = randn(N,N)
            R = nearest_rotation(M)
            V = M/R
            U = R\M

            # R is rotation matrix
            @test isrotation(R)
            # U, V are (approximately) symmetric matrices
            # See [polar decomposition](https://en.wikipedia.org/wiki/Polar_decomposition)
            @test U ≈ U'
            @test V ≈ V'
            # Eigen values are non-negative, except for the first value.
            # We just need sort for Julia 1.0
            @test all(sort(eigvals(U))[2:end] .≥ 0)
            @test all(sort(eigvals(V))[2:end] .≥ 0)
        end

        for _ in 1:100
            M = @SMatrix randn(N,N)
            R = nearest_rotation(M)
            V = M/R
            U = R\M

            @test isrotation(R)
            @test U ≈ U'
            @test V ≈ V'
            # Can't calculate eigen values of general SMatrix.
            @test all(sort(eigvals(Symmetric(U)))[2:end] .≥ 0)
            @test all(sort(eigvals(Symmetric(V)))[2:end] .≥ 0)
        end
    end
end


text:
using JuMP, Gurobi

function SMORE(env, edges, capacity, flows, demand, T, Tf)

    nedges = size(edges,1)
    nflows = size(flows,1)
    ntunnels = size(T,1)

    #CREATE TUNNEL EDGE MATRIX
    L = zeros(ntunnels, nedges)
    for t in 1:ntunnels
        for e in 1:nedges
            if in(e, T[t])
                L[t,e] = 1
            end
        end
    end

    model = Model(solver=GurobiSolver(env, OutputFlag=0))
    @variable(model, Z >= 0, basename="Z")
    @variable(model, a[1:nflows,1:size(Tf[1],1)] >= 0, basename="a")

    for f in 1:nflows
        for t in 1:size(Tf[f],1)
        @constraint(model, sum(a[f,t] for t in 1:size(Tf[f],1)) == 1)   #the sum of all allocated bandwidths on every flow must be >= the total bandwidth for that flow
        end

    end

    @expression(model, U[e=1:nedges], sum(a[f,t] * L[Tf[f][t],e] for f in 1:nflows, t in 1:size(Tf[f],1)) / capacity[e])
    for e in 1:nedges
        @constraint(model, U[e] <= Z)   #overlapping flows cannot add up to the capacity of that link
    end

    @objective(model, Min, Z)
    solve(model)
    return getvalue(a)
end


text:
@reexport module ValDateFilters

using TSML.TSMLTypes
import TSML.TSMLTypes.fit! # to overload
import TSML.TSMLTypes.transform! # to overload
using TSML.Utils

using Impute
using Impute: interp, locf, nocb
using Dates
using DataFrames
using Statistics
using CSV
using CodecBzip2

using MLDataUtils: slidingwindow

export fit!,transform!

export Matrifier,Dateifier
export DateValizer,DateValgator,DateValNNer,DateValMultiNNer
export CSVDateValReader, CSVDateValWriter, DateValLinearImputer
export BzCSVDateValReader

const gAggDict = Dict(
    :median => Statistics.median,
    :mean =>   Statistics.mean,
    :maximum => Statistics.maximum,
    :minimum => Statistics.minimum,
    :sum => sum
)

"""
    Matrifier(Dict(
       Dict(
        :ahead => 1,
        :size => 7,
        :stride => 1,
      )
    )

Converts a 1-D timeseries into sliding window matrix for ML training:
- `:ahead` => steps ahead to predict
- `:size` => size of sliding window
- `:stride` => amount of overlap in sliding window

Example:

    mtr = Matrifier(Dict(:ahead=>24,:size=>24,:stride=>5))
    lower = DateTime(2017,1,1)
    upper = DateTime(2017,1,5)
    dat=lower:Dates.Hour(1):upper |> collect
    vals = 1:length(dat)
    x = DataFrame(Date=dat,Value=vals)
    fit!(mtr,x)
    res = transform!(mtr,x)

Implements: `fit!`, `transform`
"""
mutable struct Matrifier <: Transformer
  model
  args

  function Matrifier(args=Dict())
    default_args = Dict{Symbol,Any}(
        :ahead => 1,
        :size => 7,
        :stride => 1,
    )
    new(nothing,mergedict(default_args,args))
  end
end

"""
    fit!(mtr::Matrifier,xx::T,y::Vector=Vector()) where {T<:Union{Matrix,Vector,DataFrame}}

Checks and validate inputs are in correct structure
"""
function fit!(mtr::Matrifier,xx::T,y::Vector=Vector()) where {T<:Union{Matrix,Vector,DataFrame}}
  typeof(xx) <: DataFrame || error("input is not a dataframe")
  x = deepcopy(xx.Value)
  x isa Vector || error("data should be a vector")
  mtr.model = mtr.args
end

"""
    transform!(mtr::Matrifier,xx::T) where {T<:Union{Matrix,Vector,DataFrame}}

Applies the parameters of sliding windows to create the corresponding matrix
"""
function transform!(mtr::Matrifier,xx::T) where {T<:Union{Matrix,Vector,DataFrame}}
  typeof(xx) <: DataFrame || error("input is not a dataframe")
  x = deepcopy(xx.Value)
  x isa Vector || error("data should be a vector")
  mtype = eltype(x)
  res=toMatrix(mtr,x)
  resarray=convert(Array{mtype},res) |> DataFrame
  rename!(resarray,names(resarray)[end] => :output)
end

function toMatrix(mtr::Transformer, x::Vector)
  stride=mtr.args[:stride];sz=mtr.args[:size];ahead=mtr.args[:ahead]
  @assert stride>0 && sz>0 && ahead > 0
  xlength = length(x)
  xlength > sz || error("data too short for the given size of sliding window")
  ndx=collect(xlength:-1:1)
  mtuples = slidingwindow(i->(i-ahead),ndx,sz,stride)
  height=size(mtuples)[1]
  mmatrix = Array{Union{DateTime,<:Real},2}(zeros(height,sz+1))
  ctr=1
  gap = xlength - mtuples[1][2][1]
  for (s,k) in mtuples
    v = [reverse(s);k] .+ gap
    mmatrix[ctr,:].=x[v]
    ctr+=1
  end
  mmatrix
end

### ====

"""
    Dateifier(args=Dict())
       Dict(
        :ahead => 1,
        :size => 7,
        :stride => 1
       )
    )

Converts a 1-D date series into sliding window matrix for ML training

Example:

    dtr = Dateifier(Dict())
    lower = DateTime(2017,1,1)
    upper = DateTime(2018,1,31)
    dat=lower:Dates.Day(1):upper |> collect
    vals = rand(length(dat))
    x=DataFrame(Date=dat,Value=vals)
    fit!(dtr,x)
    res = transform!(dtr,x)

Implements: `'fit!`, `transform!`
"""
mutable struct Dateifier <: Transformer
  model
  args

  function Dateifier(args=Dict())
    default_args = Dict{Symbol,Any}(
        :ahead => 1,
        :size => 7,
        :stride => 1
    )
    new(nothing,mergedict(default_args,args))
  end
end

"""
    fit!(dtr::Dateifier,xx::T,y::Vector=[]) where {T<:Union{Matrix,Vector,DataFrame}}

Computes range of dates to be used during transform.
"""
function fit!(dtr::Dateifier,xx::T,y::Vector=[]) where {T<:Union{Matrix,Vector,DataFrame}}
  typeof(xx) <: DataFrame || error("input not a dataframe")
  x = deepcopy(xx.Date)
  (eltype(x) <: DateTime || eltype(x) <: Date) || error("array element types are not dates")
  dtr.args[:lower] = minimum(x)
  dtr.args[:upper] = maximum(x)
  dtr.model = dtr.args
end

"""
    transform!(dtr::Dateifier,xx::T) where {T<:Union{Matrix,Vector,DataFrame}}

Transforms to day of the month, day of the week, etc
"""
function transform!(dtr::Dateifier,xx::T) where {T<:Union{Matrix,Vector,DataFrame}}
  typeof(xx) <: DataFrame || error("input not a dataframe")
  x = deepcopy(xx.Date)
  x isa Vector || error("data should be a vector")
  @assert eltype(x) <: DateTime || eltype(x) <: Date
  res=toMatrix(dtr,x)
  endpoints = convert(Array{DateTime},res)[:,end-1]
  dt = DataFrame()
  dt.year=Dates.year.(endpoints)
  dt.month=Dates.month.(endpoints)
  dt.day=Dates.day.(endpoints)
  dt.hour=Dates.hour.(endpoints)
  dt.week=Dates.week.(endpoints)
  dt.dow=Dates.dayofweek.(endpoints)
  dt.doq=Dates.dayofquarter.(endpoints)
  dt.qoy=Dates.quarterofyear.(endpoints)
  dtr.args[:header] = names(dt)
  #convert(Matrix{Int64},dt)
  return dt
end

"""
    DateValgator(args=Dict())
       Dict(
        :dateinterval => Dates.Hour(1),
        :aggregator => :median
      )
    )

Aggregates values based on date period specified.

Example:

    # generate random values with missing data
    Random.seed!(123)
    gdate = DateTime(2014,1,1):Dates.Minute(15):DateTime(2016,1,1)
    gval = Array{Union{Missing,Float64}}(rand(length(gdate)))
    gmissing = 50000
    gndxmissing = Random.shuffle(1:length(gdate))[1:gmissing]
    X = DataFrame(Date=gdate,Value=gval)
    X.Value[gndxmissing] .= missing

    dtvlmean = DateValgator(Dict(
          :dateinterval=>Dates.Hour(1),
          :aggregator => :mean))
    fit!(dtvlmean,X)
    res = transform!(dtvlmean,X)

Implements: `fit!`, `transform!`
"""
mutable struct DateValgator <: Transformer
  model
  args
  function DateValgator(args=Dict())
    default_args = Dict{Symbol,Any}(
        :dateinterval => Dates.Hour(1),
        :aggregator => :median
    )
    new(nothing,mergedict(default_args,args))
  end
end

function validdateval!(x::T) where {T<:DataFrame}
  size(x)[2] == 2 || error("Date Val timeseries need two columns")
  (eltype(x[:,1]) <: DateTime || eltype(x[:,1]) <: Date) || error("array element types are not dates")
  eltype(x[:,2]) <: Union{Missing,Real} || error("array element types are not values")
  cnames = names(x)
  rename!(x,Dict(cnames[1]=>:Date,cnames[2]=>:Value))
end

"""
    fit!(dvmr::DateValgator,xx::T,y::Vector=[]) where {T<:Union{Matrix,DataFrame}}

Checks and validates arguments.
"""
function fit!(dvmr::DateValgator,xx::T,y::Vector=[]) where {T<:Union{Matrix,DataFrame}}
  x = deepcopy(xx)
  validdateval!(x)
  aggr = dvmr.args[:aggregator]
  aggr in keys(gAggDict) || error("aggregator function passed is not recognized: ",aggr)
  dvmr.model=dvmr.args
end

"""
    transform!(dvmr::DateValgator,xx::T) where {T<:DataFrame}

Aggregates values grouped by date-time period using aggregate
function such as mean, median, maximum, minimum. Default is mean.
"""
function transform!(dvmr::DateValgator,xx::T) where {T<:DataFrame}
  x = deepcopy(xx)
  validdateval!(x)
  # make sure aggregator function exists
  aggr = dvmr.args[:aggregator]
  aggr in keys(gAggDict) || error("aggregator function passed is not recognized: ",aggr)
  # get the Statistics function
  aggfn = gAggDict[aggr]
  # pass the aggregator function to the generic aggregator function
  fn = aggregatorclskipmissing(aggfn)
  grpby = typeof(dvmr.args[:dateinterval])
  sym = Symbol(grpby)
  x[!,sym] = round.(x.Date,grpby)
  aggr=by(x,sym,MeanValue = :Value=>fn)
  rename!(aggr,Dict(names(aggr)[1]=>:Date,names(aggr)[2]=>:Value))
  lower = round(minimum(x.Date),grpby)
  upper = round(maximum(x.Date),grpby)
  #create list of complete dates and join with aggregated data
  cdate = DataFrame(Date = collect(lower:dvmr.args[:dateinterval]:upper))
  joined = join(cdate,aggr,on=:Date,kind=:left)
  joined
end

"""
    DateValizer(
       Dict(
        :medians => DataFrame(),
        :dateinterval => Dates.Hour(1)
      )
    )

Normalizes and cleans time series by replacing `missings` with global medians
computed based on time period groupings.

Example:

    # generate random values with missing data
    Random.seed!(123)
    gdate = DateTime(2014,1,1):Dates.Minute(15):DateTime(2016,1,1)
    gval = Array{Union{Missing,Float64}}(rand(length(gdate)))
    gmissing = 50000
    gndxmissing = Random.shuffle(1:length(gdate))[1:gmissing]
    X = DataFrame(Date=gdate,Value=gval)
    X.Value[gndxmissing] .= missing

    dvzr = DateValizer(Dict(:dateinterval=>Dates.Hour(1)))
    fit!(dvzr,X)
    transform!(dvzr,X)

Implements: `fit!`, `transform!`
"""
mutable struct DateValizer <: Transformer
  model
  args

  function DateValizer(args=Dict())
    default_args = Dict{Symbol,Any}(
        :medians => DataFrame(),
        :dateinterval => Dates.Hour(1)
    )
    new(nothing,mergedict(default_args,args))
  end
end

function getMedian(t::Type{T},xx::DataFrame) where {T<:Union{TimePeriod,DatePeriod}}
  x = deepcopy(xx)
  sgp = Symbol(t)
  fn = Dict(Dates.Second=>Dates.second,
            Dates.Minute=>Dates.minute,
            Dates.Hour=>Dates.hour,
            Dates.Day=>Dates.day,
            Dates.Month=>Dates.month)
  try
    x[!,sgp]=fn[t].(x.Date)
  catch
    error("unknown dateinterval")
  end
  gpmeans = by(x,sgp,Value = :Value => skipmedian)
  gpmeans
end

function fullaggregate!(dvzr::DateValizer,xx::T) where {T<:DataFrame}
  x = deepcopy(xx)
  grpby = typeof(dvzr.args[:dateinterval])
  sym = Symbol(grpby)
  x[!,sym] = round.(x.Date,grpby)
  aggr = by(x,sym,MeanValue = :Value=>skipmedian)
  rename!(aggr,Dict(names(aggr)[1]=>:Date,names(aggr)[2]=>:Value))
  lower = minimum(x.Date)
  upper = maximum(x.Date)
  #create list of complete dates and join with aggregated data
  cdate = DataFrame(Date = collect(lower:dvzr.args[:dateinterval]:upper))
  joined = join(cdate,aggr,on=:Date,kind=:left)
  joined
end

"""
    fit!(dvzr::DateValizer,xx::T,y::Vector=[]) where {T<:DataFrame}

Validates input and computes global medians grouped by time period.
"""
function fit!(dvzr::DateValizer,xx::T,y::Vector=[]) where {T<:DataFrame}
  x = deepcopy(xx)
  validdateval!(x)
  # get complete dates and aggregate
  joined = fullaggregate!(dvzr,x)
  grpby = typeof(dvzr.args[:dateinterval])
  sym = Symbol(grpby)
  medians = getMedian(grpby,joined)
  dvzr.args[:medians] = medians
  dvzr.model=dvzr.args
end

"""
    transform!(dvzr::DateValizer,xx::T) where {T<:DataFrame}

Replaces `missing` with the corresponding global medians with respect to time period.
"""
function transform!(dvzr::DateValizer,xx::T) where {T<:DataFrame}
  x = deepcopy(xx)
  validdateval!(x)
  # get complete dates, aggregate, and get medians
  joined = fullaggregate!(dvzr,x)
  # copy medians
  medians = dvzr.args[:medians]
  grpby = typeof(dvzr.args[:dateinterval])
  sym = Symbol(grpby)
  fn = Dict(Dates.Hour=>Dates.hour,
            Dates.Minute=>Dates.minute,
            Dates.Second=>Dates.second,
            Dates.Day => Dates.day,
            Dates.Month=>Dates.month)
  try
    joined[!,sym]=fn[grpby].(joined.Date)
  catch
    error("unknown dateinterval")
  end
  # find indices of missing
  missingndx = findall(ismissing.(joined.Value))
  jmndx=joined[missingndx,sym] .+ 1 # get time period index of missing, convert 0 index time to 1 index
  missingvals::SubArray = @view joined[missingndx,:Value]
  missingvals .= medians[jmndx,:Value] # replace missing with median value
  sum(ismissing.(joined.Value)) == 0 || error("Aggregation by time period failed to replace missings")
  joined[:,[:Date,:Value]]
end

"""
    DateValNNer(
       Dict(
          :missdirection => :symmetric, #:reverse, # or :forward or :symmetric
          :dateinterval => Dates.Hour(1),
          :nnsize => 1,
          :strict => true,
          :aggregator => :median
      )
    )

Fills `missings` with their nearest-neighbors.
- `:missdirection` => direction to fill missing data (:symmetric, :reverse, :forward)
- `:dateinterval` => time period to use for grouping,
- `:nnsize` => neighborhood size,
- `:strict` => boolean value to indicate whether to be strict about replacement or not,
- `:aggregator => function to aggregate based on date interval

Example:

    Random.seed!(123)
    gdate = DateTime(2014,1,1):Dates.Minute(15):DateTime(2016,1,1)
    gval = Array{Union{Missing,Float64}}(rand(length(gdate)))
    gmissing = 50000
    gndxmissing = Random.shuffle(1:length(gdate))[1:gmissing]
    X = DataFrame(Date=gdate,Value=gval)
    X.Value[gndxmissing] .= missing

    dnnr = DateValNNer(Dict(
          :dateinterval=>Dates.Hour(1),
          :nnsize=>10,
          :missdirection => :symmetric,
          :strict=>true,
          :aggregator => :mean))
    fit!(dnnr,X)
    transform!(dnnr,X)

Implements: `fit!`, transform!`
"""
mutable struct DateValNNer <: Transformer
  model
  args

  function DateValNNer(args=Dict())
    default_args = Dict{Symbol,Any}(
        :missdirection => :symmetric, #:reverse, # or :forward or :symmetric
        :dateinterval => Dates.Hour(1),
        :nnsize => 1,
        :strict => true,
        :aggregator => :median
    )
    new(nothing,mergedict(default_args,args))
  end
end

"""
    fit!(dnnr::DateValNNer,xx::T,y::Vector=[]) where {T<:DataFrame}

Validates and checks arguments for errors.
"""
function fit!(dnnr::DateValNNer,xx::T,y::Vector=[]) where {T<:DataFrame}
  x = deepcopy(xx)
  validdateval!(x)
  aggr = dnnr.args[:aggregator]
  aggr in keys(gAggDict) || error("aggregator function passed is not recognized: ",aggr)
  dnnr.model=dnnr.args
end

"""
    transform!(dnnr::DateValNNer,xx::T) where {T<:DataFrame}

Replaces `missings` by nearest neighbor looping over the dataset until all missing values are gone.
"""
function transform!(dnnr::DateValNNer,xx::T) where {T<:DataFrame}
  x = deepcopy(xx)
  validdateval!(x)
  # make sure aggregator function exists
  aggr = dnnr.args[:aggregator]
  aggr in keys(gAggDict) || error("aggregator function pass is not recognized: ",aggr)
  # get the Statistics function
  aggfn = gAggDict[aggr]
  # pass the aggregator function to the generic aggregator function
  fn = aggregatorclskipmissing(aggfn)
  grpby = typeof(dnnr.args[:dateinterval])
  sym = Symbol(grpby)
  # aggregate by time period
  x[!,sym] = round.(x.Date,grpby)
  aggr = by(x,sym,MeanValue = :Value=>fn)
  rename!(aggr,Dict(names(aggr)[1]=>:Date,names(aggr)[2]=>:Value))
  lower = round(minimum(x.Date),grpby)
  upper = round(maximum(x.Date),grpby)
  #create list of complete dates and join with aggregated data
  cdate = DataFrame(Date = collect(lower:dnnr.args[:dateinterval]:upper))
  joined = join(cdate,aggr,on=:Date,kind=:left)
  missingcount = sum(ismissing.(joined.Value))
  dnnr.args[:missingcount] = missingcount
  res = transform_worker!(dnnr,joined)
  count=1
  if dnnr.args[:missdirection] == :symmetric
    while sum(ismissing.(res.Value)) > 0
      res = transform_worker!(dnnr,res)
      count += 1
    end
  end
  dnnr.args[:loopcount] = count
  res
end

function transform_worker!(dnnr::DateValNNer,joinc::T) where {T<:DataFrame}
  joined = deepcopy(joinc)
  maxrow = size(joined)[1]

  # to fill-in with nearest neighbors
  nnsize::Int64 = dnnr.args[:nnsize]
  themissing = findall(ismissing.(joined.Value))
  # ==== symmetric nearest neighbor
  missingndx = DataFrame()
  if dnnr.args[:missdirection] == :symmetric
    missed = themissing |> reverse
    missingndx.Missed = missed
    # get lower:upper range
    missingndx.neighbors = map(missingndx.Missed) do m
      lower = (m-nnsize >= 1) ? (m-nnsize) : 1
      upper = (m+nnsize <= maxrow) ? m+nnsize : maxrow
      lower:upper
    end
  else
    # ===== reverse and forward
    missed = (dnnr.args[:missdirection] == :reverse) ? (themissing |> reverse) : themissing
    missingndx.Missed = missed
    # dealing with boundary exceptions, default to range until the maxrow
    missingndx.neighbors = (m->((m+1>=maxrow) || (m+nnsize>=maxrow)) ? (m+1:maxrow) : (m+1:m+nnsize)).(missingndx.Missed) # NN ranges
  end
  #joined[missingndx[:Missed],:Value] = (r -> skipmedian(joined[r,:Value])).(missingndx[:neighbors]) # iterate to each range
  missingvals::SubArray = @view joined[missingndx.Missed,:Value] # get view of only missings
  missingvals .=  (r -> skipmedian(joined[r,:Value])).(missingndx.neighbors) # replace with nn medians
  dnnr.args[:strict] && (sum(ismissing.(joined.Value)) == 0 || error("Nearest Neigbour algo failed to replace missings"))
  joined
end

"""
    CSVDateValReader(
       Dict(
          :filename => "",
          :dateformat => ""
       )
    )

Reads csv file and parse date using the given format.
- `:filename` => complete path including filename of csv file
- `:dateformat` => date format to parse

Example:

    inputfile =joinpath(dirname(pathof(TSML)),"../data/testdata.csv")
    csvreader = CSVDateValReader(Dict(:filename=>inputfile,:dateformat=>"d/m/y H:M"))
    fit!(csvreader)
    df = transform!(csvreader)

    # using pipeline workflow
    filter1 = DateValgator()
    filter2 = DateValNNer(Dict(:nnsize=>1))
    mypipeline = Pipeline(Dict(
          :transformers => [csvreader,filter1,filter2]
      )
    )
    fit!(mypipeline)
    res=transform!(mypipeline)

Implements: `fit!`, `transform!`
"""
mutable struct CSVDateValReader <: Transformer
    model
    args
    function CSVDateValReader(args=Dict())
        default_args = Dict(
            :filename => "",
            :dateformat => ""
        )
        new(nothing,mergedict(default_args,args))
    end
end

"""
    fit!(csvrdr::CSVDateValReader,x::T=[],y::Vector=[]) where {T<:Union{DataFrame,Vector,Matrix}}

Makes sure filename and dateformat are not empty strings.
"""
function fit!(csvrdr::CSVDateValReader,x::T=[],y::Vector=[]) where {T<:Union{DataFrame,Vector,Matrix}}
    fname = csvrdr.args[:filename]
    fmt = csvrdr.args[:dateformat]
    (fname != "" && fmt != "") || error("missing filename or date format")
    csvrdr.model = csvrdr.args
end

"""
    transform!(csvrdr::CSVDateValReader,x::T=[]) where {T<:Union{DataFrame,Vector,Matrix}}

Uses CSV package to read the csv file and converts it to dataframe.
"""
function transform!(csvrdr::CSVDateValReader,x::T=[]) where {T<:Union{DataFrame,Vector,Matrix}}
    fname = csvrdr.args[:filename]
    fmt = csvrdr.args[:dateformat]
    df = CSV.read(fname) |> DataFrame
    ncol(df) == 2 || error("dataframe should have only two columns: Date,Value")
    rename!(df,names(df)[1]=>:Date,names(df)[2]=>:Value)
    if !(eltype(df.Date) <: DateTime )
      df.Date = DateTime.(df.Date,fmt)
    end
    df
end

"""
    CSVDateValWriter(
       Dict(
          :filename => "",
          :dateformat => ""
       )
    )

Writes the time series dataframe into a file with the given date format.

Example:

    inputfile =joinpath(dirname(pathof(TSML)),"../data/testdata.csv")
    outputfile = joinpath("/tmp/test.csv")
    csvreader = CSVDateValReader(Dict(:filename=>inputfile,:dateformat=>"d/m/y H:M"))
    csvwtr = CSVDateValWriter(Dict(:filename=>outputfile,:dateformat=>"d/m/y H:M"))
    filter1 = DateValgator()
    filter2 = DateValNNer(Dict(:nnsize=>1))
    mypipeline = Pipeline(Dict(
          :transformers => [csvreader,filter1,filter2,csvwtr]
      )
    )
    fit!(mypipeline)
    res=transform!(mypipeline)

    # read back what was written to validate
    csvreader = CSVDateValReader(Dict(:filename=>outputfile,:dateformat=>"y-m-d HH:MM:SS"))
    fit!(csvreader)
    transform!(csvreader)

Implements: `fit!`, `transform!`
"""
mutable struct CSVDateValWriter <: Transformer
    model
    args
    function CSVDateValWriter(args=Dict())
        default_args = Dict(
            :filename => "",
            :dateformat => ""
        )
        new(nothing,mergedict(default_args,args))
    end
end

"""
    fit!(csvwtr::CSVDateValWriter,x::T=[],y::Vector=[]) where {T<:Union{DataFrame,Vector,Matrix}}

Makes sure filename and dateformat are not empty strings.
"""
function fit!(csvwtr::CSVDateValWriter,x::T=[],y::Vector=[]) where {T<:Union{DataFrame,Vector,Matrix}}
    fname = csvwtr.args[:filename]
    fmt = csvwtr.args[:dateformat]
    fname != ""  || error("missing filename")
    csvwtr.model = csvwtr.args
end

"""
    transform!(csvwtr::CSVDateValWriter,x::T) where {T<:Union{DataFrame,Vector,Matrix}}

Uses CSV package to write the dataframe into a csv file.
"""
function transform!(csvwtr::CSVDateValWriter,x::T) where {T<:Union{DataFrame,Vector,Matrix}}
    fname = csvwtr.args[:filename]
    fmt = csvwtr.args[:dateformat]
    df = deepcopy(x) |> DataFrame
    if ncol(df) == 2
      rename!(df,names(df)[1]=>:Date,names(df)[2]=>:Value)
      eltype(df.Date) <: DateTime || error("Date format error")
    end
    df |> CSV.write(fname)
    return df
end

"""
    BzCSVDateValReader(
       Dict(
          :filename => "",
          :dateformat => ""
       )
    )

Reads Bzipped csv file and parse date using the given format.
- `:filename` => complete path including filename of csv file
- `:dateformat` => date format to parse

Example:

    inputfile =joinpath(dirname(pathof(TSML)),"../data/testdata.csv.bz2")
    csvreader = BzCSVDateValReader(Dict(:filename=>inputfile,:dateformat=>"d/m/y H:M"))
    filter1 = DateValgator()
    filter2 = DateValNNer(Dict(:nnsize=>1))
    mypipeline = Pipeline(Dict(
          :transformers => [csvreader,filter1,filter2]
      )
    )
    fit!(mypipeline)
    res=transform!(mypipeline)

Implements: `fit!`, `transform!`
"""
mutable struct BzCSVDateValReader <: Transformer
    model
    args
    function BzCSVDateValReader(args=Dict())
        default_args = Dict(
            :filename => "",
            :dateformat => ""
        )
        new(nothing,mergedict(default_args,args))
    end
end

"""
    fit!(bzcsvrdr::BzCSVDateValReader,x::T=[],y::Vector=[]) where {T<:Union{DataFrame,Vector,Matrix}}

Makes sure filename and dateformat are not empty strings.
"""
function fit!(bzcsvrdr::BzCSVDateValReader,x::T=[],y::Vector=[]) where {T<:Union{DataFrame,Vector,Matrix}}
    fname = bzcsvrdr.args[:filename]
    fmt = bzcsvrdr.args[:dateformat]
    (fname != "" && fmt != "") || error("missing filename or date format")
    bzcsvrdr.model = bzcsvrdr.args
end

"""
    transform!(bzcsvrdr::BzCSVDateValReader,x::T=[]) where {T<:Union{DataFrame,Vector,Matrix}}

Uses CodecBzip2 package to read the csv file and converts it to dataframe.
"""
function transform!(bzcsvrdr::BzCSVDateValReader,x::T=[]) where {T<:Union{DataFrame,Vector,Matrix}}
    fname = bzcsvrdr.args[:filename]
    fmt = bzcsvrdr.args[:dateformat]
    stream = Bzip2DecompressorStream(open(fname))
    df = CSV.read(stream) |> DataFrame
    ncol(df) == 2 || error("dataframe should have only two columns: Date,Value")
    rename!(df,names(df)[1]=>:Date,names(df)[2]=>:Value)
    df.Date = DateTime.(df.Date,fmt)
    df
end

"""
    DateValMultiNNer(
       Dict(
          :type => :knn # :linear
          :missdirection => :symmetric, #:reverse, # or :forward or :symmetric
          :dateinterval => Dates.Hour(1),
          :nnsize => 1,
          :strict => true,
          :aggregator => :median
      )
    )

Fills `missings` with their nearest-neighbors. It assumes that first column is a Date class
and the other columns are Union{Missings,Real}. It uses DateValNNer and DateValizer+Impute to
process each numeric column concatendate with the Date column.
- `:type` => type of imputation which can be a linear interpolation or nearest neighbor
- `:missdirection` => direction to fill missing data (:symmetric, :reverse, :forward)
- `:dateinterval` => time period to use for grouping,
- `:nnsize` => neighborhood size,
- `:strict` => boolean value to indicate whether to be strict about replacement or not,
- `:aggregator => function to aggregate based on date interval

Example:

    Random.seed!(123)
    gdate = DateTime(2014,1,1):Dates.Minute(15):DateTime(2016,1,1)
    gval1 = Array{Union{Missing,Float64}}(rand(length(gdate)))
    gval2 = Array{Union{Missing,Float64}}(rand(length(gdate)))
    gval3 = Array{Union{Missing,Float64}}(rand(length(gdate)))
    gmissing = 50000
    gndxmissing1 = Random.shuffle(1:length(gdate))[1:gmissing]
    gndxmissing2 = Random.shuffle(1:length(gdate))[1:gmissing]
    gndxmissing3 = Random.shuffle(1:length(gdate))[1:gmissing]
    X = DataFrame(Date=gdate,Temperature=gval1,Humidity=gval2,Ozone=gval3)
    X.Temperature[gndxmissing1] .= missing
    X.Humidity[gndxmissing2] .= missing
    X.Ozone[gndxmissing3] .= missing

    dnnr = DateValMultiNNer(Dict(
          :type=>:linear,
          :dateinterval=>Dates.Hour(1),
          :nnsize=>10,
          :missdirection => :symmetric,
          :strict=>true,
          :aggregator => :mean))
    fit!(dnnr,X)
    transform!(dnnr,X)

Implements: `fit!`, transform!`
"""
mutable struct DateValMultiNNer <: Transformer
  model
  args

  function DateValMultiNNer(args=Dict())
    default_args = Dict{Symbol,Any}(
        :type => :knn,
        :missdirection => :symmetric, #:reverse, # or :forward or :symmetric
        :dateinterval => Dates.Hour(1),
        :nnsize => 1,
        :strict => true,
        :aggregator => :median
    )
    new(nothing,mergedict(default_args,args))
  end
end

function multivalidateval(x::DataFrame)
  size(x)[2] > 2 || error("Multi Date Val timeseries need more than two columns")
  (eltype(x[:,1]) <: DateTime || eltype(x[:,1]) <: Date) || error("array element types are not dates")
  sum(broadcast(y->eltype(y)<:Union{Missing,Real},eachcol(x[:,2:end]))) == ncol(x)-1 || error("columns 2:end should be real numbers")
end

"""
    fit!(dnnr::DateValMultiNNer,xx::T,y::Vector=[]) where {T<:DataFrame}

Validates and checks arguments for errors.
"""
function fit!(dnnr::DateValMultiNNer,xx::T,y::Vector=[]) where {T<:DataFrame}
  x = deepcopy(xx)
  # validate it's multi-dimensional and first column is date
  multivalidateval(x)
  cnames = names(x)
  rename!(x,Dict(cnames[1]=>:Date))
  aggr = dnnr.args[:aggregator]
  aggr in keys(gAggDict) || error("aggregator function passed is not recognized: ",aggr)
  dnnr.model=dnnr.args
end

"""
    transform!(dnnr::DateValMultiNNer,xx::T) where {T<:DataFrame}

Replaces `missings` by nearest neighbor or linear interpolation by looping over the dataset
for each column until all missing values are gone.
"""
function transform!(dnnr::DateValMultiNNer,xx::T) where {T<:DataFrame}
  x = deepcopy(xx)
  # make sure data is valid
  multivalidateval(x)
  # loop each column and call ValdateNNer to impute
  df = DataFrame()
  if dnnr.args[:type] == :knn
    df = knnimpute(dnnr,x)
  elseif dnnr.args[:type] == :linear
    df = linearimpute(dnnr,x)
  end
  return df
end

function knnimpute(dnnr::DateValMultiNNer,x::DataFrame)
  valnner = DateValNNer(dnnr.args)
  df = DataFrame(Date=x.Date)
  cnames = names(x)
  for y in eachcol(x[:,2:end])
    input = DataFrame(Date=x.Date,Value=y)
    fit!(valnner,input)
    res=transform!(valnner,input)
    df = join(df,res,on=:Date,makeunique=true)
  end
  names!(df,cnames)
  return df
end

function linearimpute(dnnr::DateValMultiNNer,x::DataFrame)
  valgator = DateValgator(dnnr.args)
  linearputer = DateValLinearImputer(dnnr.args)
  df = DataFrame(Date=x.Date)
  cnames = names(x)
  for y in eachcol(x[:,2:end])
    input = DataFrame(Date=x.Date,Value=y)
    fit!(valgator,input)
    agg=transform!(valgator,input)
    fit!(linearputer,agg)
    res=transform!(linearputer,agg)
    df = join(df,res,on=:Date,makeunique=true)
  end
  names!(df,cnames)
  return df
end

"""
    DateValLinearImputer(
       Dict(
          :dateinterval => Dates.Hour(1),
      )
    )

Fills `missings` by linear interpolation.
- `:dateinterval` => time period to use for grouping,

Example:

    Random.seed!(123)
    gdate = DateTime(2014,1,1):Dates.Minute(15):DateTime(2016,1,1)
    gval = Array{Union{Missing,Float64}}(rand(length(gdate)))
    gmissing = 50000
    gndxmissing = Random.shuffle(1:length(gdate))[1:gmissing]
    X = DataFrame(Date=gdate,Value=gval)
    X.Value[gndxmissing] .= missing

    dnnr = DateValLinearImputer()
    fit!(dnnr,X)
    transform!(dnnr,X)

Implements: `fit!`, transform!`
"""
mutable struct DateValLinearImputer <: Transformer
  model
  args

  function DateValLinearImputer(args=Dict())
    default_args = Dict{Symbol,Any}(
        :dateinterval => Dates.Hour(1),
    )
    new(nothing,mergedict(default_args,args))
  end
end

"""
    fit!(dnnr::DateValLinearImputer,xx::T,y::Vector=[]) where {T<:DataFrame}

Validates and checks arguments for errors.
"""
function fit!(dnnr::DateValLinearImputer,xx::T,y::Vector=[]) where {T<:DataFrame}
  x = deepcopy(xx)
  validdateval!(x)
  dnnr.model=dnnr.args
end

"""
    transform!(dnnr::DateValLinearImputer,xx::T) where {T<:DataFrame}

Replaces `missings` by linear interpolation.
"""
function transform!(dnnr::DateValLinearImputer,xx::T) where {T<:DataFrame}
  x = deepcopy(xx)
  validdateval!(x)
  valgator = DateValgator(dnnr.args)
  fit!(valgator,x)
  df=transform!(valgator,x)
  df.Value = Impute.interp(df.Value) |> Impute.locf() |> Impute.nocb()
  return df
end

end


text:
macro def(name, definition) # thanks to http://www.stochasticlifestyle.com/type-dispatch-design-post-object-oriented-programming-julia/
    return esc(quote
        macro $name()
            esc($(Expr(:quote, definition)))
        end
    end)
end

macro cthreads(condition::Symbol,loop::Expr) #does not work well because of #15276, but seems to work on Julia v0.7
    return esc(quote
        if $condition
            Threads.@threads $loop
        else
            $loop
        end
    end)
end

macro mthreads(ex)
    if THR
       return esc(:(Threads.@threads $ex))
    else
        return esc(ex)
    end
end

macro msimd(ex)
    return esc(ex)
    #return esc(:(@simd $ex))
end

macro gen(ex)
    ex2 = macroexpand(ex)
    return esc(:(@generated $ex2))
end

const sim_par = (:EquationType,:VelocityTimeStepType,:PassiveScalarType,:DensityStratificationType,:LESModelType,:ForcingType,:HyperViscosityType)
# Nx,Ny,Nz size of the grid in Fourier Space
# Lcs = Nx*Ny*Nz; Lcv = Lcs*3
# Nrx size of x direction in Real Space skipping padding. The same as 1:(:nx) in global file
# Nrrx size of x direction in Real Space NOT skipping padding.
# Lrs = (Nrx+2)*Ny*Nz; Lrv = 3*Nrx
# kx,ky,kz wavenumber vectors.
# Gdirec can be :x, :y, :z and gives the gravity direction. Odly, only if :z seems to work
# Thr is true or false for threads.
# Nt number of threads.

const sim_fields = (:passivescalar,:densitystratification,:lesmodel,:forcing,:hyperviscosity)

"""
    @par(t::Symbol)
Add the parameters defined in `sim_par` to the Symbol.

# Examples
```julia-repl
julia> sim_par = (:Nx,:Ny,:Nz,:Lcs,:Lcv,:Nrx,:Lrs,:Lrv,:Tr,:Integrator)
(:Nx, :Ny, :Nz, :Lcs, :Lcv, :Nrx, :Lrs, :Lrv, :Tr, :Integrator)

julia> @macroexpand @par(Parameters)
:(Parameters{Nx, Ny, Nz, Lcs, Lcv, Nrx, Lrs, Lrv, Tr, Integrator})

julia> @macroexpand @par(BoussinesqParameters)
:(BoussinesqParameters{Nx, Ny, Nz, Lcs, Lcv, Nrx, Lrs, Lrv, Tr, Integrator})
```
"""
macro par(t::Symbol)
    return esc(Expr(:curly,t,sim_par...))
end

"""
---
    @par(ex:Expr)
Add the parameters defined in `sim_par` after a `where` on a function definition.

# Examples
```julia-repl
julia> sim_par = (:Nx,:Ny,:Nz,:Lcs,:Lcv,:Nrx,:Lrs,:Lrv,:Tr,:Integrator)
(:Nx, :Ny, :Nz, :Lcs, :Lcv, :Nrx, :Lrs, :Lrv, :Tr, :Integrator)

julia> @macroexpand @par function test(s::@par(Parameters),b::Array{T,N}) where {T,N}; return 1; end
:(function test(s::Parameters{Nx, Ny, Nz, Lcs, Lcv, Nrx, Lrs, Lrv, Tr, Integrator}, b::Array{T, N}) where {Nx, Ny, Nz, Lcs, Lcv, Nrx, Lrs, Lrv, Tr, Integrator, T, N} # REPL[1], line 1:
return 1
end)
```
"""
macro par(ex::Expr)
    if ex.head == :function || ex.head == :(=)
        if ex.args[1].head == :where
            # append!(ex.args[1].args,sim_par)
            ex.args[1].args = vcat(ex.args[1].args[1], sim_par..., ex.args[1].args[2:end]) # This way the parameters are available to the existing parameters
        elseif ex.args[1].head == :call
            ex.args[1] = Expr(:where, ex.args[1], sim_par...)
        end
    end
    return esc(ex)
end


text:
struct ResourceManagement
    "External resources used by one or more shaders for a given widget."
    resources::Vector{Any}
    "Whether resources needs an update."
    needs_update::Function
    "Function that updates resources."
    update_resources::Function
end

struct WidgetDependencies
    shaders::Vector{ShaderSpecification}
    resource_management::Optional{ResourceManagement}
end

WidgetDependencies(shaders) = WidgetDependencies(shaders, nothing)

"""
Make the link between a description of widgets and rendering.
"""
struct GUIRenderer
    device::Created{Device,DeviceCreateInfo}
    disp::QueueDispatch
    swapchain_info::SwapchainCreateInfoKHR
    descriptor_allocator::DescriptorAllocator
    command_pools::ThreadedCommandPool
    shader_cache::ShaderCache
    descriptor_set_layout_cache::DescriptorSetLayoutCache
    pipeline_layout_cache::PipelineLayoutCache
    render_set::RenderSet
end

device(gr::GUIRenderer) = handle(gr.device)

function GUIRenderer(device, disp::QueueDispatch)
    GUIRenderer(
        device,
        disp,
        DescriptorAllocator(device),
        ThreadedCommandPool(device, disp, dictionary([1 => QUEUE_COMPUTE_BIT | QUEUE_GRAPHICS_BIT])),
        ShaderCache(device),
        DescriptorSetLayoutCache(device),
        PipelineLayoutCache(device),
        RenderSet(device),
    )
end

function add_widget!(gr::GUIRenderer, w::Widget, widget_dependencies::WidgetDependencies)
    dset_layouts = get!(gr.descriptor_set_layout_cache, widget_dependencies.shaders)
    dsets = allocate_descriptor_sets!(gr.descriptor_allocator, dset_layouts)
    append!(shader_dependencies.descriptor_sets, dsets)
    insert!(gr.render_set, w, RenderInfo(shader_dependencies, (nvertices(w),1,0,0); push_data = push_data(w)))
    resources = widget_dependencies.resource_management.resources
    !isempty(resources) && update_descriptor_sets(device(gr), resources)
    nothing
end

function needs_resource_update(gr::GUIRenderer)
    any(gr.widget_dependencies) do deps
        isnothing(deps.resource_management) && return false
        rm = deps.resource_management
        rm.needs_update(rm.resources...)
    end
end

function update_resources(gr::GUIRenderer)
    foreach(zip(gr.shader_dependencies, gr.widget_dependencies)) do (sdeps, wdeps)
        !isnothing(wdeps.resource_management) || return
        rm = wdeps.resource_management
        if rm.needs_update(resources...)
            rm.resources .= rm.update(resources...)
        end
        update_descriptor_sets(device(gr), sdeps)
    end
end


text:
#=
BibText
@TECHREPORT{Benner95acollection,
author = {Peter Benner and Alan J. Laub and Volker Mehrmann},
title = {A Collection of Benchmark Examples for the Numerical Solution of
Algebraic Riccati Equations II: Discrete-Time Case},
institution = {FAK. F. MATHEMATIK, TU CHEMNITZ--ZWICKAU},
year = {1995}
}

=#

# ------------------------------------------------------------------
# ejm1 collection of benchmark examples
A1 = [4. 3.;-4.5 -3.5]
B1 = [1. -1]'
R1 = 1.
Q1 = [9. 6.;6. 4.]
X1 = (1+sqrt(5))/2*[9 6;6 4] # theoretical solucion

X11 = dare(A1,B1,Q1,R1) # computed solution

# --------------------------------------------------------------------
# ejm computed with mathematica
Am = [1. -1 1;0 1 1;0 0 1]
Bm = [1. 0;1 0;0 1]
Rm = [10 0;0 0.1]
Qm = [10 0 0;0 1 0;0 0 0.1]
Xm = [42.2835 -68.5247 -3.94783;-68.5247 154.043 16.0017;-3.94783 16.0017 8.33197] # mathematica solution

Xmm = dare(Am,Bm,Qm,Rm) # computed solution

#####

print("Testing dare: ")
if (sum((X1-X11).^2) < 1e-6) && (sum((Xm-Xmm).^2) < 1e-6)
	println("OK")
end


text:
using LazyArtifacts
using PkgUtility
using PlotPlants
using PyPlot
using Test

println("");
@testset "PlotPlants --- canvas related" begin
    _fig,_axes = create_canvas("New canvas", ncol=2, nrow=2);
    set_titles!(_axes);
    set_titles!(_axes; usetex=false);
    set_xylabels!(_axes, ["X" for i in 1:4], ["Y" for i in 1:4]);
    set_xylabels!(_axes, "X", "Y");
    set_xylims!(_axes, [[0,1] for i in 1:4], [[0,1] for i in 1:4]);
    set_xylims!(_axes, [0,1], [0,1]);
    set_xyticks!(_axes, [[0,1] for i in 1:4], [[0,1] for i in 1:4]);
    set_xyticks!(_axes, [0,1], [0,1]);
    set_xyticklabels!(_axes, [["",""] for i in 1:4], [["",""] for i in 1:4]);
    set_xyticklabels!(_axes, ["",""], ["",""]);
    save_canvas!(_fig, "_test.pdf", true);
    rm("_test.pdf");
    @test true;
    close("all");
    _figs = Figure[];
    push!(_figs, deepcopy(_fig));
    push!(_figs, deepcopy(_fig));
    save_gif!(_figs, 1, "test.gif");
    rm("test.gif");
    @test true;
end

println("");
@testset "PlotPlants --- plots related" begin
    _fig,_axes = create_canvas("New canvas"; ncol=3, nrow=2);
    _ax1,_ax2,_ax3,_ax4,_ax5,_ax6 = _axes;

    # density plot
    plot_density!(_ax1, rand(100), rand(100));
    @test true;

    # linear regression
    plot_line_regress!(_ax2, rand(100), rand(100), interval=true);
    @test true;

    # preview data
    preview_data(rand(5), rand(5));
    predownload_artifact!.(["CH_20X_1Y_V1","LAI_4X_1M_V1"], "../Artifacts.toml");
    preview_dataset!(_ax3, artifact"LAI_4X_1M_V1" * "/LAI_4X_1M_V1.nc", "LAI");
    preview_dataset!(_ax3, artifact"CH_20X_1Y_V1" * "/CH_20X_1Y_V1.nc", "Band1");
    @info "Expecting a warning here!";
    preview_dataset!(_ax3, artifact"CH_20X_1Y_V1" * "/CH_20X_1Y_V1.nc", 1);
    @test true;

    # artists
    plot_ellipse!(_ax4, (0,0));
    plot_stoma!(_ax4, (0,0));
    @test true;

    # hexbin
    plot_hexbin!(_ax5, rand(100), rand(100));
    plot_hexbin!(_ax5, rand(100), rand(100); logbins=true);
    plot_hexbin!(_ax6, rand(100), rand(100), [-1,1], [-1,1]);
    close("all");
end

println("");
@testset "PlotPlants --- regression related" begin
    line_regress(rand(100), rand(100); intercept=false);
    line_regress(rand(100), rand(100); intercept=true);
    @test true;
end

println("");
@testset "PlotPlants --- statistics related" begin
    # calculate density
    calculate_density(rand(100), rand(100));
    calculate_density(rand(1001), rand(1001));
    calculate_density(rand(100), rand(100), 2);
    calculate_density(rand(1001), rand(1001), 2);
    @test true;

    # test slope
    test_slope(rand(100), rand(100); intercept=false);
    test_slope(rand(100), rand(100); intercept=true);
    @test true;
end

println("");
@testset "PlotPlants --- string related" begin
    # symbol
    latex_symbol("H", sub="sub", sup="sup");
    latex_symbol("HaHa", presub="presub", presup="presup");
    @test true;

    # unit
    latex_unit("A");
    latex_unit("E");
    latex_unit("E_MMOL");
    latex_unit("G");
    latex_unit("PAR");
    latex_unit("T");
    latex_unit("WUE"; squared=true, space=true);
    latex_unit("WUE"; squared=false, space=false);
    latex_unit("LAI");
    @test true;

    # parse symbol
    syms = ["alpha", "Alpha", "beta", "Beta", "gamma", "Gamma", "delta",
            "Delta", "epsilon", "Epsilon", "zeta", "Zeta", "eta", "Eta",
            "theta", "Theta", "iota", "Iota", "kappa", "Kappa", "lambda",
            "Lambda", "mu", "Mu", "nu", "Nu", "xi", "Xi", "pi", "Pi", "rho",
            "Rho", "sigma", "Sigma", "tau", "Tau", "upsilon", "Upsilon", "phi",
            "Phi", "chi", "Chi", "psi", "Psi", "omega", "Omega", "deg", "degc",
            "nothing"]
    parse_symbol.(syms, false);
    parse_symbol.(syms, true);
    @test true;

    # sub and superscripts
    subscript(-1);
    subscript(-1; standalone=false);
    superscript(-1);
    superscript(-1; standalone=false);
    @test true;
end

println("");
@testset "PlotPlants --- string related" begin
    # fonts
    use_sans();
    use_serif();
    @test true;

    # fonts and renders
    use_sans_tex();
    use_sans_text();
    use_serif_tex();
    use_serif_text();
    @test true;
end


text:
using Turing
using Test
using Random

setadsafe(true)

include("../utility.jl")
Random.seed!(100)

D = [1.0 1.0 4.0 4.0]

@model MoGtest(D) = begin
  mu1 ~ Normal(1, 1)
  mu2 ~ Normal(4, 1)
  z1 ~ Categorical(2)
  if z1 == 1
    D[1] ~ Normal(mu1, 1)
  else
    D[1] ~ Normal(mu2, 1)
  end
  z2 ~ Categorical(2)
  if z2 == 1
    D[2] ~ Normal(mu1, 1)
  else
    D[2] ~ Normal(mu2, 1)
  end
  z3 ~ Categorical(2)
  if z3 == 1
    D[3] ~ Normal(mu1, 1)
  else
    D[3] ~ Normal(mu2, 1)
  end
  z4 ~ Categorical(2)
  if z4 == 1
    D[4] ~ Normal(mu1, 1)
  else
    D[4] ~ Normal(mu2, 1)
  end
  z1, z2, z3, z4, mu1, mu2
end

gibbs = Gibbs(500, PG(10, 1, :z1, :z2, :z3, :z4), HMC(3, 0.15, 3, :mu1, :mu2))
chain = sample(MoGtest(D), gibbs)

@test mean(chain[:z1]) ≈ 1.0 atol=0.1
@test mean(chain[:z2]) ≈ 1.0 atol=0.1
@test mean(chain[:z3]) ≈ 2.0 atol=0.1
@test mean(chain[:z4]) ≈ 2.0 atol=0.1
@test mean(chain[:mu1]) ≈ 1.0 atol=0.1
@test mean(chain[:mu2]) ≈ 4.0 atol=0.1

setadsafe(false)


text:
using StaticTaylorSeries
using Test

@testset "StaticTaylorSeries.jl" begin
    t = StaticTaylor(1, 2, 3)

    @test t + 1 == StaticTaylor(2, 2, 3)
    @test t - 1 == StaticTaylor(0, 2, 3)
    @test t * 2 == StaticTaylor(2, 4, 6)

    @test t^2 == StaticTaylor(1, 4, 10)

    t = StaticTaylor(1.0, 2.0, 3.0)
    @test t / 2 == StaticTaylor(0.5, 1.0, 1.5)
end

max_stars_repo_name: dpsanders/StaticTaylorSeries

max_issues_repo_name: dpsanders/StaticTaylorSeries

max_forks_repo_name: dpsanders/StaticTaylorSeries


text:
"""

$(TYPEDEF)

Structure to contain the density values obtained from the calculation.

$(TYPEDFIELDS)

"""
@with_kw mutable struct Density
    solute::Float64 = 0.0
    solvent::Float64 = 0.0
    solvent_bulk::Float64 = 0.0
end

function reset!(d::Density)
    d.solute = 0.0
    d.solvent = 0.0
    d.solvent_bulk = 0.0
    return nothing
end

#function Base.show(io::IO, d::Density )
#  println(" Mean solute density: $(d.solute) ")
#  println(" Mean solvent density: $(d.solvent) ")
#  println(" Mean solvent bulk density: $(d.solvent_bulk) ")
#end


text:
include("../../LibSPP/librarySPP.jl")
include("./spp_dm3.jl")
include("../../../setSPP.jl")

function main()
    println("Etudiant(e)s : LE et OUSMAN")

    # Collecting the names of instances to solve located in the folder Data ----
    target = "../dat"
    fnames = getfname(target)
    println()

    fres = splitdir(splitdir(pwd())[end-1])[end]
    io = open("../res/"*fres*".res", "w")
    for instance = 1:length(fnames)

        C, A = loadSPP(string(target,"/",fnames[instance]))

        iterMax,iterPenaliser,length_tm = 100000,rand(35:45),rand(5:9)

        println("Instance : ",fnames[instance])
        println("iterMax = ", iterMax, " | iterPenaliser = ", iterPenaliser, " | |TM| init = ", length_tm)
        sleep(5)        #pause 5s pour lire les informations

        t1 = @elapsed xInit, zInit, Einit = GreedyConstruction(C,A)
        println("z(x0) = ",zInit)
        t2 = @elapsed xBest, zBest, length_tm2 = Tabou(C,A,iterMax,iterPenaliser,length_tm,xInit,zInit,Einit)

        println(io, "Instance : ",fnames[instance])
        println(io, "iterMax = ", iterMax, " | iterPenaliser = ", iterPenaliser, " | |TM| init = ", length_tm, " | |TM| final = ", length_tm2)
        if isAdmissible(C,A,xBest)
            println(io, "Feasible : yes | Σ(x_j) = ", length(findall(isequal(1),xBest)), " | z(x0) = ", zInit, " | z(x) = ", zBest)
        else
            println(io, "Feasible : no")
        end
        println(io, "CPU time = ",t1 + t2)
        println(io)
        println()
    end
    close(io)

end

main()

max_stars_repo_name: ducanh-le/metaSPPstu

max_issues_repo_name: ducanh-le/metaSPPstu

max_forks_repo_name: ducanh-le/metaSPPstu


text:
############################################
################ Structures ################
############################################

function glb(types...)
    # If a lower bound is in the types then it is greatest
    # else we just return Unknown for now
    for type in types
        all(type <: t_ for t_ in types) && return type
    end
    return Unknown
end

function input_target_scitypes(models, metalearner)
    # The target scitype is defined as the greatest lower bound of the
    # metalearner and the base models in the library
    all_tg_scitypes = [target_scitype(m) for m in models]
    tg_scitype = glb(target_scitype(metalearner), all_tg_scitypes...)
    # The input scitype is defined as the greatest lower bound of the
    # base models in the library
    inp_scitype = glb([input_scitype(m) for m in models]...)

    return inp_scitype, tg_scitype
end

mutable struct DeterministicStack{modelnames, inp_scitype, tg_scitype} <: DeterministicComposite
   models::Vector{Supervised}
   metalearner::Deterministic
   resampling
   function DeterministicStack(modelnames, models, metalearner, resampling)
        inp_scitype, tg_scitype = input_target_scitypes(models, metalearner)
        return new{modelnames, inp_scitype, tg_scitype}(models, metalearner, resampling)
   end
end

mutable struct ProbabilisticStack{modelnames, inp_scitype, tg_scitype} <: ProbabilisticComposite
    models::Vector{Supervised}
    metalearner::Probabilistic
    resampling
    function ProbabilisticStack(modelnames, models, metalearner, resampling)
        inp_scitype, tg_scitype = input_target_scitypes(models, metalearner)
        return new{modelnames, inp_scitype, tg_scitype}(models, metalearner, resampling)
    end
 end

const Stack{modelnames, inp_scitype, tg_scitype} =
    Union{DeterministicStack{modelnames, inp_scitype, tg_scitype},
            ProbabilisticStack{modelnames, inp_scitype, tg_scitype}}

"""
    Stack(;metalearner=nothing, resampling=CV(), named_models...)

Implements the generalized Stack algorithm introduced by Wolpert
in https://www.sciencedirect.com/science/article/abs/pii/S0893608005800231 and
generalized by Van der Laan et al in https://biostats.bepress.com/ucbbiostat/paper222/.

Instead of using your favorite model, use them all! The Stack is a metalearning algorithm
covered by theoretical guarantees.

The stack in a nutshell:

- The data is split in training/validation sets
- Each model in the library is trained on each training set and outputs predictions on the validation sets
- The metalearner is subsequently trained on those predictions and finds the best combination of the models
in the library.
- Each model is then retrained on the full data
- "Stacking" those models and the metalearner results in an end to end fully trained model

You are not exempt from evaluating the stack.

We currently provide two different stack types the `DeterministicStack` and the `ProbabilisticStack`.
The type of which is automatically chosen by the constructor based on the provided metalearner.

# Arguments
- `metalearner::Model`: The model that will optimize the desired criterion based on its internals.
                        For instance, a LinearRegression model will optimize the squared error.
- `resampling::Union{CV, StratifiedCV}`: The resampling strategy used to train the metalearner.
- `named_models`: The models that will be part of the library

# Examples

Let's build a simple DeterministicStack for a continuous target, we show that:

- Some members of the library may very well be `Probabilistic` models even though the stack
is `Deterministic`, the expected value will be taken.
- Rather than running a hyperparameter search, you can integrate each model in the stack

```julia
using MLJ

DecisionTreeRegressor = @load DecisionTreeRegressor pkg=DecisionTree
EvoTreeRegressor = @load EvoTreeRegressor
XGBoostRegressor = @load XGBoostRegressor
KNNRegressor = @load KNNRegressor pkg=NearestNeighborModels
LinearRegressor = @load LinearRegressor pkg=MLJLinearModels

X, y = make_regression(500, 5)

stack = Stack(;metalearner=LinearRegressor(),
                resampling=CV(),
                constant=ConstantRegressor(),
                tree_2=DecisionTreeRegressor(max_depth=2),
                tree_3=DecisionTreeRegressor(max_depth=3),
                evo=EvoTreeRegressor(),
                knn=KNNRegressor(),
                xgb=XGBoostRegressor())

mach = machine(stack, X, y)
evaluate!(mach; resampling=Holdout(), measure=rmse)

```

"""
function Stack(;metalearner=nothing, resampling=CV(), named_models...)
    metalearner === nothing &&
        throw(ArgumentError("No metalearner specified. Use Stack(metalearner=...)"))

    nt = NamedTuple(named_models)
    modelnames = keys(nt)
    models = collect(nt)

    if metalearner isa Deterministic
        stack =  DeterministicStack(modelnames, models, metalearner, resampling)
    elseif metalearner isa Probabilistic
        stack = ProbabilisticStack(modelnames, models, metalearner, resampling)
    else
        throw(ArgumentError("The metalearner should be a subtype
                    of $(Union{Deterministic, Probabilistic})"))
    end
    # Issuing clean! statement
    message = MMI.clean!(stack)
    isempty(message) || @warn message

    # Warning if either input_scitype/target_scitype is
    # Unknown at construction time
    params = typeof(stack).parameters
    params[end-1] == Unknown && @warn "Could not infer input_scitype of the stack"
    params[end] == Unknown && @warn "Could not infer target_scitype of the stack"

    return stack
end

function MMI.clean!(stack::Stack{modelnames, inp_scitype, tg_scitype}) where {modelnames,inp_scitype,tg_scitype}
    # We only carry checks and don't try to correct the arguments here
    message = ""
    # Checking target_scitype and input_scitype have not been changed from the original stack
    glb_inp_scitype, glb_tg_scitype = input_target_scitypes(getfield(stack, :models), stack.metalearner)
    glb_inp_scitype == inp_scitype ||
            throw(DomainError(inp_scitype, "The newly inferred input_scitype of the stack doesn't
            match its original one. You have probably changed one of the base models or the metalearner
            to a non compatible type."))
    glb_tg_scitype == tg_scitype ||
            throw(DomainError(tg_scitype, "The newly inferred target_scitype of the stack doesn't
            match its original one. You have probably changed one of the base model or the metalearner
            to a non compatible type."))
    # Checking the target scitype is consistent with either Probabilistic/Deterministic Stack
    target_scitype(stack.metalearner) <: Union{AbstractArray{<:Continuous}, AbstractArray{<:Finite}} ||
        throw(ArgumentError("The metalearner should have target_scitype:
                $(Union{AbstractArray{<:Continuous}, AbstractArray{<:Finite}})"))

    return message
end

Base.propertynames(::Stack{modelnames}) where modelnames = tuple(:resampling, :metalearner, modelnames...)

function Base.getproperty(stack::Stack{modelnames}, name::Symbol) where modelnames
    name === :metalearner && return getfield(stack, :metalearner)
    name === :resampling && return getfield(stack, :resampling)
    models = getfield(stack, :models)
    for j in eachindex(modelnames)
        name === modelnames[j] && return models[j]
    end
    error("type Stack has no property $name")
end

function Base.setproperty!(stack::Stack{modelnames}, _name::Symbol, val) where modelnames
    _name === :metalearner && return setfield!(stack, :metalearner, val)
    _name === :resampling && return setfield!(stack, :resampling, val)
    idx = findfirst(==(_name), modelnames)
    idx isa Nothing || return getfield(stack, :models)[idx] = val
    error("type Stack has no property $name")
end

MMI.target_scitype(::Type{<:Stack{modelnames, input_scitype, target_scitype}}) where
    {modelnames, input_scitype, target_scitype} = target_scitype

MMI.input_scitype(::Type{<:Stack{modelnames, input_scitype, target_scitype}}) where
    {modelnames, input_scitype, target_scitype} = input_scitype

MLJBase.load_path(::Type{<:ProbabilisticStack}) = "MLJBase.ProbabilisticStack"
MLJBase.load_path(::Type{<:DeterministicStack}) = "MLJBase.DeterministicStack"
MLJBase.package_name(::Type{<:Stack}) = "MLJBase"
MLJBase.package_uuid(::Type{<:Stack}) = "a7f614a8-145f-11e9-1d2a-a57a1082229d"
MLJBase.package_url(::Type{<:Stack}) = "https://github.com/alan-turing-institute/MLJBase.jl"
MLJBase.package_license(::Type{<:Stack}) = "MIT"

###########################################################
################# Node operations Methods #################
###########################################################

getfolds(y::AbstractNode, cv::CV, n::Int) =
    source(train_test_pairs(cv, 1:n))

getfolds(y::AbstractNode, cv::StratifiedCV, n::Int) =
    node(YY->train_test_pairs(cv, 1:n, YY), y)

trainrows(X::AbstractNode, folds::AbstractNode, nfold) =
    node((XX, ff) -> selectrows(XX, ff[nfold][1]), X, folds)

testrows(X::AbstractNode, folds::AbstractNode, nfold) =
    node((XX, ff) -> selectrows(XX, ff[nfold][2]), X, folds)

pre_judge_transform(ŷ::Node, ::Type{<:Probabilistic}, ::Type{<:AbstractArray{<:Finite}}) =
    node(ŷ -> pdf(ŷ, levels(first(ŷ))), ŷ)

pre_judge_transform(ŷ::Node, ::Type{<:Probabilistic}, ::Type{<:AbstractArray{<:Continuous}}) =
    node(ŷ->mean.(ŷ), ŷ)

pre_judge_transform(ŷ::Node, ::Type{<:Deterministic}, ::Type{<:AbstractArray{<:Continuous}}) =
    ŷ

function oos_set(m::Stack, folds::AbstractNode, Xs::Source, ys::Source)
    Zval = []
    yval = []
    # Loop over the cross validation folds to build a training set for the metalearner.
    for nfold in 1:m.resampling.nfolds
        Xtrain = trainrows(Xs, folds, nfold)
        ytrain = trainrows(ys, folds, nfold)
        Xtest = testrows(Xs, folds, nfold)
        ytest = testrows(ys, folds, nfold)

        # Train each model on the train fold and predict on the validation fold
        # predictions are subsequently used as an input to the metalearner
        Zfold = []
        for model in getfield(m, :models)
            mach = machine(model, Xtrain, ytrain)
            ypred = predict(mach, Xtest)
            # Dispatch the computation of the expected mean based on
            # the model type and target_scytype
            ypred = pre_judge_transform(ypred, typeof(model), target_scitype(model))
            push!(Zfold, ypred)
        end

        Zfold = hcat(Zfold...)

        push!(Zval, Zfold)
        push!(yval, ytest)
    end

    Zval = MLJBase.table(vcat(Zval...))
    yval = vcat(yval...)

    Zval, yval
end

#######################################
################# Fit #################
#######################################
"""
    fit(m::Stack, verbosity::Int, X, y)
"""
function fit(m::Stack, verbosity::Int, X, y)
    n = nrows(y)

    Xs = source(X)
    ys = source(y)

    folds = getfolds(ys, m.resampling, n)

    Zval, yval = oos_set(m, folds, Xs, ys)

    metamach = machine(m.metalearner, Zval, yval)

    # Each model is retrained on the original full training set
    Zpred = []
    for model in getfield(m, :models)
        mach = machine(model, Xs, ys)
        ypred = predict(mach, Xs)
        ypred = pre_judge_transform(ypred, typeof(model), target_scitype(model))
        push!(Zpred, ypred)
    end

    Zpred = MLJBase.table(hcat(Zpred...))
    ŷ = predict(metamach, Zpred)

    # We can infer the Surrogate by two calls to supertype
    mach = machine(supertype(supertype(typeof(m)))(), Xs, ys; predict=ŷ)

    return!(mach, m, verbosity)

end