text:
# # Image Classification With Cubical Persistent Homology

# In this example, we will show how to use Ripserer in an image classification
# context. Persistent homology is not a predictive algorithm, but it can be used to extract
# useful features from data.

# ## Setting up

using Ripserer
using PersistenceDiagrams
using Images # also required: ImageIO to read .png files
using Plots
using ProgressMeter
using Random
Random.seed!(1337)

data_dir = joinpath(@__DIR__, "../assets/data/malaria") # replace with the correct path.
nothing # hide

# Let's load the data. We will use a a [data
# set](https://lhncbc.nlm.nih.gov/publication/pub9932) with microscope images of healthy
# cells and cells infected with malaria. The original data set is quite large, but we can
# pretend we were only given 200 images to work with. We have chosen the 200 images
# randomly.

uninfected = shuffle!(Images.load.(readdir(joinpath(data_dir, "uninfected"); join=true)))
infected = shuffle!(Images.load.(readdir(joinpath(data_dir, "infected"); join=true)))

images = [uninfected; infected]
classes = [fill(false, length(uninfected)); fill(true, length(infected))]
nothing # hide

# Let's see what the images look like.

plot(
    plot(uninfected[1]; title="Healthy"),
    plot(uninfected[2]; title="Healthy"),
    plot(infected[1]; title="Infected"),
    plot(infected[2]; title="Infected"),
)

# To make the images work with Ripserer, we convert them to floating gray scale values. We
# do not have to resize the images. Maybe some additional preprocessing, such as
# normalization would help, but we'll skip it for this example.

inputs = [Gray.(image) for image in images]
nothing # hide

# Now we can compute persistence diagrams. Since we are working with images, we have to use
# the `Cubical` filtration type. Cubical persistent homology should detect the dark spots
# (local minima) in the images. It's pretty efficient, so this should only take a few
# seconds.

diagrams = @showprogress [ripserer(Cubical(i)) for i in inputs]

# This is what some of the diagrams look like.

plot(plot(images[1]; title="Healthy"), plot(diagrams[1]))

#

plot(plot(images[end]; title="Infected"), plot(diagrams[end]))

# Notice that there is a lot more going on in the middle of the infected diagram, especially
# in ``H_0``.

# The persistence diagrams might look nice, but are hard to use with machine learning
# algorithms. The number of points in the diagram may be different for every image, even
# when images are of the same size. We can solve this problem by using a vectorization
# method, such as converting all diagrams to persistence images.

# Persistence images work by weighting each point in the diagram with a distribution. The
# distribution defaults to a Gaussian, but any function of two arguments can be used. Each
# point is also weighted by a weighting function that should be equal to zero along the
# ``x``-axis. It defaults to a function that is zero on the ``x``-axis and linearly
# increases to the maximum persistence in the diagram.

# We start by splitting the diagrams into their 0 and 1 dimensional components.

dim_0 = first.(diagrams)
dim_1 = last.(diagrams)
nothing; # hide

# We feed the diagram to the `PersistenceImage` constructor which will choose ranges that
# will fit all the diagrams. We set the sigma value to 0.1, since all persistence pairs are
# in the ``[0,1]×[0,1]`` square and the default sigma of 1 would be too wide. We will use
# the default image size, which is 5×5.

image_0 = PersistenceImage(dim_0; sigma=0.1)

#

image_1 = PersistenceImage(dim_1; sigma=0.1)

# Let's see how some of the images look like.

plot(plot(dim_0[end]; persistence=true), heatmap(image_0(dim_0[end]); aspect_ratio=1))

#

plot(plot(dim_1[end]; persistence=true), heatmap(image_1(dim_1[end]); aspect_ratio=1))

# Next, we convert all diagrams to images and use `vec` to turn them into flat vectors. We
# then concatenate the zero and one-dimensional images. The result is a vector of length 50
# for each diagram.

persims = [[vec(image_0(dim_0[i])); vec(image_1(dim_1[i]))] for i in 1:length(diagrams)]

# ## Fitting A Model

# Now it's time to fit our model. We will use
# [GLMNet.jl](https://github.com/JuliaStats/GLMNet.jl) to fit a regularized linear model.

using GLMNet

# Convert the image vectors to a matrix that will be understood by `glmnet`.

X = reduce(hcat, persims)'
y = classes
nothing; # hide

# Start by randomly splitting the data into two sets, a training and a testing set.

perm = shuffle(1:200)
train_x = X[perm[1:100], :]
train_y = y[perm[1:100]]
test_x = X[perm[101:end], :]
test_y = y[perm[101:end]]
nothing; # hide

# Fit the model and predict.

path = glmnet(train_x, train_y)
cv = glmnetcv(train_x, train_y)

λ = path.lambda[argmin(cv.meanloss)]
path = glmnet(train_x, train_y; lambda=[λ])

predictions = .!iszero.(round.(GLMNet.predict(path, test_x)))
nothing; # hide

# Get the classification accuracy.

count(predictions .== test_y) / length(test_y)

# Not half bad considering we haven't touched the images and we left pretty much all
# settings on default.

# Now let's look at the misclassified examples.

missed = findall(predictions .!= test_y)
label = ("Healthy", "Infected")
plts = [plot(images[i]; title="$(label[test_y[i] + 1])", ticks=nothing) for i in missed]
plot(plts...)

# Finally, let's look at which parts of the persistence images `glmnet` considered important.

plot(
    heatmap(reshape(path.betas[1:25], (5, 5)); title="H₀ coefficients"),
    heatmap(reshape(path.betas[26:50], (5, 5)); title="H₁ coefficients"),
)

# These correspond to the area we identified at the beginning. Also note that in this case,
# the classifier does not care about ``H_1`` at all.

# ## Using MLJ

# Another, more straightforward way to execute a similar pipeline is to use Ripserer's
# [MLJ.jl](https://github.com/alan-turing-institute/MLJ.jl) integration. We will use a
# random forest classifier for this example.

# We start by loading MLJ and the classifier. Not that
# [MLJDecisionTreeInterface.jl](https://github.com/bensadeghi/DecisionTree.jl) needs to be
# installed for this to work.

using MLJ
tree = @load RandomForestClassifier pkg = "DecisionTree" verbosity = 0

# We create a pipeline of `CubicalPersistentHomology` followed by the classifier. In this
# case, `CubicalPersistentHomology` takes care of both the homology computation and the
# conversion to persistence images.

pipe = @pipeline(CubicalPersistentHomology(), tree)

# We train the pipeline the same way you would fit any other MLJ model. Remember, we need to
# use grayscale versions of images stored in `inputs`.

classes = coerce(classes, Binary)
train, test = partition(eachindex(classes), 0.7; shuffle=true, rng=1337)
mach = machine(pipe, inputs, classes)
fit!(mach; rows=train)

# Next, we predict the classes on the test data and print out the classification accuracy.

yhat = predict_mode(mach, inputs[test])
accuracy(yhat, classes[test])

# The result is quite a bit worse than before. We can try mitigating that by using a
# different vectorizer.

pipe.cubical_persistent_homology.vectorizer = PersistenceCurveVectorizer()
mach = machine(pipe, inputs, classes)
fit!(mach; rows=train)

yhat = predict_mode(mach, inputs[test])
accuracy(yhat, classes[test])

# The result could be improved further by choosing a different model and
# vectorizer. However, this is just a short introduction. Please see the [MLJ.jl
# documentation](https://alan-turing-institute.github.io/MLJ.jl/dev/) for more information
# on model tuning and selection, and the [PersistenceDiagrams.jl
# documentation](https://mtsch.github.io/PersistenceDiagrams.jl/dev/mlj/) for a list of
# vectorizers and their options.


text:
using .Settings: from_yaml
using .Calculators: SingleCalculator, SamePhDOSCalculator, DifferentPhDOSCalculator

export runcode

function runcode(cfgfile)
    config = from_yaml(cfgfile)

    if !isdir(config["output_directory"])
        mkpath(config["output_directory"])
    end
    config["qha_output"] = joinpath(config["output_directory"], "output.txt")
    if isfile(config["qha_output"])
        rm(config["qha_output"])
    end
    save_to_output(config["qha_output"], make_starting_string())

    str = lowercase(config["calculation"])
    if str == "single"
        calculator = SingleCalculator(config)
        println("You have single-configuration calculation assumed.")
    elseif str == "same phonon dos"
        calculator = SamePhDOSCalculator(config)
        println(
            "You have multi-configuration calculation with the same phonon DOS assumed.",
        )
    elseif str == "different phonon dos"
        calculator = DifferentPhDOSCalculator(config)
        println(
            "You have multi-configuration calculation with different phonon DOS assumed.",
        )
    else
        error("""
              The 'calculation' in your settings in not recognized! It must be one of:
              'single', 'same phonon dos', 'different phonon dos'!
              """)
    end

    calculator.read_input()
    println("Caution: If negative frequencies found, they are currently treated as 0!")

    calculator.refine_grid()
    calculator.desired_pressure_status()
    temperature_array = calculator.temperature_array
    desired_pressures_gpa = calculator.desired_pressures_gpa
    temperature_sample = calculator.temperature_sample_array
    p_sample_gpa = calculator.pressure_sample_array

    results_folder = expanduser(config["output_directory"])
    calculation_option = Dict(
        "F" => "f_tp",
        "G" => "g_tp",
        "H" => "h_tp",
        "U" => "u_tp",
        "V" => "v_tp",
        "Cv" => "cv_tp_jmolk",
        "Cp" => "cp_tp_jmolk",
        "Bt" => "bt_tp_gpa",
        "Btp" => "btp_tp",
        "Bs" => "bs_tp_gpa",
        "alpha" => "alpha_tp",
        "gamma" => "gamma_tp",
    )

    file_ftv_fitted = joinpath(results_folder, "f_tv_fitted_ev_ang3.txt")
    save_x_tv(
        calculator.f_tv_ev,
        temperature_array,
        calculator.finer_volumes_ang3,
        temperature_sample,
        file_ftv_fitted,
    )

    file_ftv_non_fitted = joinpath(results_folder, "f_tv_nonfitted_ev_ang3.txt")
    save_x_tv(
        calculator.vib_ev,
        temperature_array,
        calculator.volumes_ang3,
        temperature_sample,
        file_ftv_non_fitted,
    )

    file_ptv_gpa = joinpath(results_folder, "p_tv_gpa.txt")
    save_x_tv(
        calculator.p_tv_gpa,
        temperature_array,
        calculator.finer_volumes_ang3,
        temperature_sample,
        file_ptv_gpa,
    )

    file_stv_j = joinpath(results_folder, "s_tv_j.txt")
    save_x_tv(
        calculator.s_tv_j,
        temperature_array,
        calculator.finer_volumes_ang3,
        temperature_sample,
        file_stv_j,
    )

    for idx in config["thermodynamic_properties"]
        if idx in ("F", "G", "H", "U")
            attr_name = calculation_option[idx] * "_" * config["energy_unit"]
            file_name = attr_name * ".txt"
            file_dir = joinpath(results_folder, file_name)
            save_x_tp(
                getproperty(calculator, attr_name),
                temperature_array,
                desired_pressures_gpa,
                p_sample_gpa,
                file_dir,
            )
        elseif idx == "V"
            v_bohr3 = calculation_option[idx] * "_bohr3"
            file_name_bohr3 = v_bohr3 * ".txt"
            file_dir_au = joinpath(results_folder, file_name_bohr3)
            v_ang3 = calculation_option[idx] * "_ang3"
            file_name_ang3 = v_ang3 * ".txt"
            file_dir_ang3 = joinpath(results_folder, file_name_ang3)
            save_x_tp(
                getproperty(calculator, v_bohr3),
                temperature_array,
                desired_pressures_gpa,
                p_sample_gpa,
                file_dir_au,
            )
            save_x_tp(
                getproperty(calculator, v_ang3),
                temperature_array,
                desired_pressures_gpa,
                p_sample_gpa,
                file_dir_ang3,
            )
        elseif idx in ("Cv", "Cp", "Bt", "Btp", "Bs", "alpha", "gamma")
            attr_name = calculation_option[idx]
            file_name = attr_name * ".txt"
            file_dir = joinpath(results_folder, file_name)
            save_x_tp(
                getproperty(calculator, attr_name),
                temperature_array,
                desired_pressures_gpa,
                p_sample_gpa,
                file_dir,
            )
        end
    end
end

ext: jl


text:
import Base: angle

_dot2(v::Vector, w::Vector) = v[1] * w[1] + v[2] * w[2]

_no_intersect_error = "The two clines do not intersect; no angle can be calculated."

function angle(L::Line, M::Line)::Float64

    if slope(L) == slope(M)
        return 0
    end

    v = collect(reim(L.a - L.b))
    w = collect(reim(M.a - M.b))

    top = _dot2(v, w)
    bot = sqrt(_dot2(v, v) * _dot2(w, w))
    θ = abs(acos(top / bot))

    if θ > π / 2
        θ = π - θ
    end
    return θ

end

"""
    angle(C::Cline, D::Cline)
Return the angle between two clines. If `C` and `D` are equal, or are parallel lines, return 0.
Otherwise, if `C` and `D` do not intersect, an error is thrown.
"""
function angle(C::Cline, D::Cline)
    if C == D
        return 0
    end

    X = C ∩ D

    if length(X) == 0
        error(_no_intersect_error)
    end

    z = first(X)
    F = LFT(0, 1, 1, -z)

    CC = F(C)
    DD = F(D)
    if !(isa(CC, Line) && isa(DD, Line))
        error("This can't happen, but it did.")
    end
    return angle(CC, DD)
end


text:
@testset "decompose" begin
    @testset "x^2+y^2+z^2" begin
        S = semi_elementary_monomial
        x = S((2,0,0),(0,0,0))
        poly = semi_elementary_polynomial(3)
        push!(poly,S((0,0,0),(2,0,0)),1)
        push!(poly,S((0,0,0),(0,1,0)),-2)
        @test decompose(x) == poly
    end

    @testset "x^3+y^3+z^3" begin
        S = semi_elementary_monomial
        x = S((3,0,0),(0,0,0))
        poly = semi_elementary_polynomial(3)
        push!(poly,S((0,0,0),(3,0,0)),1)
        push!(poly,S((0,0,0),(1,1,0)),-3)
        push!(poly,S((0,0,0),(0,0,1)),3)
        @test decompose(x) == poly
    end

    @testset "x^4+y^4+z^4" begin
        S = semi_elementary_monomial
        x = S((4,0,0),(0,0,0))
        poly = semi_elementary_polynomial(3)
        push!(poly,S((0,0,0),(4,0,0)),1)
        push!(poly,S((0,0,0),(2,1,0)),-4)
        push!(poly,S((0,0,0),(1,0,1)),4)
        push!(poly,S((0,0,0),(0,2,0)),2)
        @test decompose(x) == poly
    end

    @testset "x^5+y^5+z^5" begin
        S = semi_elementary_monomial
        x = S((5,0,0),(0,0,0))
        poly = semi_elementary_polynomial(3)
        push!(poly,S((0,0,0),(5,0,0)),1)
        push!(poly,S((0,0,0),(3,1,0)),-5)
        push!(poly,S((0,0,0),(2,0,1)),5)
        push!(poly,S((0,0,0),(1,2,0)),5)
        push!(poly,S((0,0,0),(0,1,1)),-5)
        @test decompose(x) == poly
    end

    @testset "w^5+x^5+y^5+z^5" begin
        S = semi_elementary_monomial
        x = S((5,0,0,0),(0,0,0,0))
        poly = semi_elementary_polynomial(4)
        push!(poly,S((0,0,0,0),(5,0,0,0)),1)
        push!(poly,S((0,0,0,0),(3,1,0,0)),-5)
        push!(poly,S((0,0,0,0),(2,0,1,0)),5)
        push!(poly,S((0,0,0,0),(1,2,0,0)),5)
        push!(poly,S((0,0,0,0),(1,0,0,1)),-5)
        push!(poly,S((0,0,0,0),(0,1,1,0)),-5)
        @test decompose(x) == poly
    end
end

max_stars_repo_name: yushoteke/SymmetricPolynomials

max_issues_repo_name: yushoteke/SymmetricPolynomials

max_forks_repo_name: yushoteke/SymmetricPolynomials


text:
#  ------------------------------------------------------------------------------------------
#  Copyright (c) Microsoft Corporation. All rights reserved.
#  Licensed under the MIT License (MIT). See LICENSE in the repo root for license information.
#  ------------------------------------------------------------------------------------------

# Gather parameters from parameter file
const default_parameters = [
    ["_POOL_ID", "BatchPool"],
    ["_POOL_COUNT", "1"],
    ["_NODE_COUNT_PER_POOL", "1"],
    ["_POOL_VM_SIZE", "Standard_E2s_v3"],
    ["_JOB_ID" "BatchJob"],
    ["_STANDARD_OUT_FILE_NAME", "stdout.txt"],
    ["_NODE_OS_PUBLISHER", "Canonical"],
    ["_NODE_OS_OFFER", "UbuntuServer"],
    ["_NODE_OS_SKU", "18.04"],
    ["_BLOB_CONTAINER", "azureclusterlesshpctemp"],
    ["_INTER_NODE_CONNECTION", "0"],
    ["_NUM_RETRYS", "0"],
    ["_MPI_RUN", "0"],
    ["_CONTAINER", "None"],
    ["_NUM_NODES_PER_TASK", "1"],
    ["_NUM_PROCS_PER_NODE", "1"],
    ["_OMP_NUM_THREADS", "1"],
    ["_JULIA_DEPOT_PATH", "/mnt/batch/tasks/startup/wd/.julia"],
    ["_PYTHONPATH", "/mnt/batch/tasks/startup/wd/.local/lib/python3.6/site-packages"],
    ["_VERBOSE", "1"]
]

function create_parameter_dict(params, default_parameters)

    for parameter in default_parameters
        if haskey(params, parameter[1])
            parameter[2] = params[parameter[1]]
        end
    end

    default_parameters = Dict(default_parameters)
    if default_parameters["_INTER_NODE_CONNECTION"] == "0" && parse(Int, default_parameters["_NUM_NODES_PER_TASK"]) > 1
        @warn "Cannot set _NUM_NODES_PER_TASK > 1 if _INTER_NODE_CONNECTION=0. Defaulting to _NUM_NODES_PER_TASK = 1."
        default_parameters["_NUM_NODES_PER_TASK"] = "1"
    end

    return default_parameters
end

# Manage batch state
"""
    batch_show  ()

 Print the collected expressions that have been priorly tagged via `@batchdef` and `@bcast`.

 See also: [`batch_clear`](@ref)
"""
function batch_show()
    print("__resources__: \n$__resources__\n")
    print("__expressions__: \n$__expressions__\n")
    print("__packages__: \n$__packages__\n")
end

# Clear global vars
"""
    batch_clear()

 Clear all expressions from the AzureClusterlessHPC call stack. All expressions priorly tagged with `@batchdef` or `@bcast` will be cleared.

 See also: [`batch_show`](@ref)
"""
function batch_clear()
    global __expressions__ = nothing    # collect generic expressions tagged with @batchdef
    global __packages__ = nothing   # collect "using" statements tagged with @batchdef
    if ~isnothing(__credentials__)
        global __resources__ = Array{Any}(undef, length(__credentials__))
        for i=1:length(__credentials__)
            __resources__[i] = []
        end
    else
        global __resources__ = Array{Any}(undef, 0)
    end
    if length(__active_pools__) > 0
        global __active_pools__ = Array{Dict}(undef, 0)
    end
end

# Delete pool
"""
    delete_pool(; pool_id=nothing)

 Delete the batch pool that is specified in the parameter.json file (or delete the pool with the default name).

 *Optional input:*

 - `pool_id=nothing`: If no pool id is provided, the function deletes the pool name specified in
    the paramter json file, or if no file is provided, the default pool (`BatchPool`).

 *Output*

 - `Nothing`

 See also: [`delete_all_jobs`](@ref), [`delete_container!`](@ref)
"""
function delete_pool(; pool_id=nothing)
    try
        for pool in __active_pools__
            isnothing(pool_id) ? (current_pool = pool["pool_id"]) : (current_pool = pool_id)
            pool["clients"]["batch_client"].pool.delete(current_pool)
        end
    catch
        print("Pool does not exist.")
    end
end

# Shortcuts to blob functions
"""
    delete_container(; blobcontainer=nothing)

 Delete the blob container that contains temporary blob files as specified in the parameter.json file.

 *Optional input:*

 - `blobcontainer=nothing`: Blob container name. If no name is provided, the function deletes the container specified in
    the paramter json file, or if no file is provided, the default container (`azureclusterlesshpctemp`).

 *Output*

 - `Bool`

 See also: [`delete_all_jobs`](@ref), [`delete_pool!`](@ref)
"""
function delete_container(; blobcontainer=nothing)
    isnothing(blobcontainer) && (blobcontainer = __params__["_BLOB_CONTAINER"])
    for client in __clients__
        client["blob_client"].delete_container(blobcontainer)
    end
end

"""
    delete_all_jobs()

 Delete all jobs containing the current (base) job ID specified in the parameter file (or the default name).

 *Output*

 - `Nothing`

 See also: [`delete_container`](@ref), [`delete_pool!`](@ref)
"""
function delete_all_jobs(;base_id=nothing)
	isnothing(base_id) && (base_id = __params__["_JOB_ID"])
    for client in __clients__
        job_list = client["batch_client"].job.list(raw=true)
        while true
            try
                job = job_list.next()
                job_id = job.id
                if job_id[1:length(base_id)] == base_id
                    client["batch_client"].job.delete(job_id)
                end
            catch
                break
            end
        end
    end
end

"""
    Options(; job_name="batchjob_", task_name="task_", priority=0)

 Specify job options for batch jobs.

 *Input"

 - `job_name`: Name of the batch job.

 - `task_name`: Base string for task names.

 - `priority` (Int): Job priority

 *Output*

 - `Options` data structure.

"""
struct Options
    job_name::String
    task_name::String
    priority::Integer
    pool::Union{Nothing, Integer}
    reset_mpi::Bool
end

Options(; job_name="batchjob_", task_name="task_", priority=0, pool=nothing, reset_mpi=false) = Options(job_name, task_name, priority, pool, reset_mpi)

# Include generic text files (e.g. python files) with task
function fileinclude(s::String)
    return nothing
end


text:
using DataStructures
world_age() = ccall(:jl_get_world_counter, Int, ())
const global pull_queue = Queue{Signal}()
const global push_queue = Queue{Tuple{Signal, Any}}()
const global eventloop_cond = Condition()

function empty_queues()
    empty!(pull_queue.store)
    empty!(push_queue.store)
end

function __init__()
    @async eventloop()
end

function run_till_now()
    while !isempty(pull_queue)
        foreach(pull!, pull_queue)
        empty!(pull_queue.store)
        if !isempty(push_queue)
            (s, val) = dequeue!(push_queue)
            soft_push!(s, val)
        end
    end
end

function eventloop(eventloop_world_age = world_age())
    try
        while true
            if !isempty(pull_queue)
                if world_age() > eventloop_world_age
                    debug_mode() && println("restarting Signals eventloop")
                    @async eventloop()
                    break
                end
                debug_mode() && println("pull-queue length: $(length(pull_queue))")
                run_till_now()
            end
            wait(eventloop_cond)
        end
    catch e
        st = stacktrace(catch_backtrace())
        empty_queues()
        @async eventloop()
        showerror(stderr, SignalException(e,st))
    end
end


text:
# PGFPlots.save("logo.tex", fig.o, include_preamble=true)
# add \begin{axis}[ticks=none,height...
# add begin{tikzpicture}[scale=0.5]
# Font: https://github.com/JuliaGraphics/julia-logo-graphics/blob/master/font/TamilMN-Bold.ttf
# pdftopng test.pdf logo.png
# convert logo.pdf logo.svg
using MonteCarloMeasurements, Plots, Random, KernelDensity, Colors
darker_blue = RGB(0.251, 0.388, 0.847)
lighter_blue = RGB(0.4, 0.51, 0.878)
darker_purple = RGB(0.584, 0.345, 0.698)
lighter_purple  = RGB(0.667, 0.475, 0.757)
darker_green  = RGB(0.22, 0.596, 0.149)
lighter_green  = RGB(0.376, 0.678, 0.318)
darker_red  = RGB(0.796, 0.235, 0.2)
lighter_red  = RGB(0.835, 0.388, 0.361)
N = 30
Random.seed!(3)
a,b,c = [0.2randn(N,2) for _ in 1:3]

a .+= [0 0]
b .+= [1 0]
c .+= [0.5 1]

pa,pb,pc = Particles.((a,b,c))

opts = (markersize=9, markerstrokewidth=2, markerstrokealpha=0.8, markeralpha=0.8, size=(300,300))
scatter(eachcol(a)...; c=lighter_red, markerstrokecolor=darker_red, opts..., axis=false, grid=false, legend=false)
scatter!(eachcol(b)...; c=lighter_purple, markerstrokecolor=darker_purple, opts...)
# scatter!(eachcol(c)...; c=lighter_green, markerstrokecolor=darker_, opts...)

ls = (linewidth=10,markerstrokewidth=2)
plot!(pa[1:1],pa[2:2]; c=lighter_red, markerstrokecolor=darker_red, ls...)
plot!(pb[1:1],pb[2:2]; c=lighter_purple, markerstrokecolor=darker_purple, ls...)
# plot!(pc[1:1],pc[2:2]; c=lighter_green, markerstrokecolor=darker_, ls...)

##
x = LinRange(-0.5,1.45,30)
f(x) = (x-0.3)^2 + 0.5
y = f.(x)
plot!(x,y, l=(darker_blue,))

bi = 1:3:N
plot!([b[bi,1] b[bi,1]]', [b[bi,2] f.(b[bi,1])]', l=(lighter_purple, :dash, 0.2))
plot!([b[bi,1] fill(0.5,length(bi))]', [f.(b[bi,1]) f.(b[bi,1])]', l=(lighter_green, :dash, 0.2))

scatter!(fill(0.5,N), f.(b[:,1]); c=lighter_green, markerstrokecolor=darker_green, opts...)

kd = kde(f.(b[:,1]), npoints=200, bandwidth=0.09, boundary=(0.5,1.8))
fig = plot!(0.5 .- 0.2kd.density, kd.x, c=lighter_green, markerstrokecolor=darker_green, fill=true, fillalpha=0.2)
# PGFPlots.save("test.tex", fig.o, include_preamble=true)


text:
using ClassicalOrthogonalPolynomials, LazyArrays, QuasiArrays, BandedMatrices, ForwardDiff, Test
import ClassicalOrthogonalPolynomials: recurrencecoefficients, jacobimatrix, Clenshaw
import QuasiArrays: MulQuasiArray

@testset "Legendre" begin
    @testset "LegendreWeight" begin
        w = LegendreWeight()
        @test w.^2 isa LegendreWeight
        @test sqrt.(w) isa LegendreWeight
        @test w .* w isa LegendreWeight
    end

    @testset "basics" begin
        P = Legendre()
        @test axes(P) == (Inclusion(ChebyshevInterval()),oneto(∞))
        @test P == P == Legendre{Float32}()
        A,B,C = recurrencecoefficients(P)
        @test B isa Zeros
        P = Jacobi(0.0,0.0)
        A,B,C = recurrencecoefficients(P)
        @test B[1] == 0.0
        x = axes(P,1)
        X = P \ (x .* P)
        @testset "recurrence coefficient and jacobimatrix" begin
            @test 1/A[1] ≈ X[2,1]
            @test -B[1]/A[1] ≈ X[1,1]
            @test C[2]/A[2] ≈ X[1,2]
            @test 1/A[2] ≈ X[3,2]
            @test -B[2]/A[2] ≈ X[2,2]
            @test C[3]/A[3] ≈ X[2,3]

            @test A[1] ≈ 1/X[2,1]
            @test B[1] ≈ -X[1,1]/X[2,1]
            @test C[2] ≈ X[1,2]/X[3,2]
            @test A[2] ≈ 1/X[3,2]
            @test B[2] ≈ X[2,2]/X[3,2]
            @test C[3] ≈ X[2,3]/X[4,3]
        end
    end

    @testset "operators" begin
        P = Legendre()
        P̃ = Jacobi(0.0,0.0)
        P̄ = Ultraspherical(1/2)

        @test jacobimatrix(P̃)[1,1] == 0.0
        @test jacobimatrix(P̃)[1:10,1:10] == jacobimatrix(P)[1:10,1:10] == jacobimatrix(P̄)[1:10,1:10]

        @test P[0.1,Base.OneTo(4)] ≈ [1,0.1,-0.485,-0.1475]
        @test P̃[0.1,1:10] ≈ P[0.1,1:10] ≈ P̄[0.1,1:10]

        @test Ultraspherical(P) == P̄
        @test Jacobi(P) == P̃

        @test P̃\P === P\P̃ === P̄\P === P\P̄ === Eye(∞)
        @test_broken P̄\P̃ === P̃\P̄ === Eye(∞)

        D = Derivative(axes(P,1))
        @test Ultraspherical(3/2)\(D*P) isa BandedMatrix{Float64,<:Ones}
        @test (Ultraspherical(5/2) \ (D^2*P))[1:10,1:10] == BandedMatrix(2 => Fill(3,8))

        P = Legendre()
        x = axes(P,1)
        w = x .+ x.^2 .+ 1
        W = P \ (w .* P)
        @test W isa Clenshaw
        @test W * [1; 2; zeros(∞)] ≈ P \ (w .* (P[:,1:2] * [1,2]))
    end

    @testset "test on functions" begin
        P = Legendre()
        D = Derivative(axes(P,1))
        f = P*Vcat(randn(10), Zeros(∞))
        P̃ = Jacobi(0.0,0.0)
        P̄ = Ultraspherical(1/2)
        @test (P̃*(P̃\f))[0.1] ≈ (P̄*(P̄\f))[0.1] ≈ f[0.1]
        C = Ultraspherical(3/2)
        @test (C*(C\f))[0.1] ≈ f[0.1]

        @test (D*f)[0.1] ≈ ForwardDiff.derivative(x -> (Legendre{eltype(x)}()*f.args[2])[x],0.1)
    end

    @testset "poly broadcast" begin
        P = Legendre()
        x = axes(P,1)

        @test (x .* P) isa MulQuasiArray

        J = P \ (x .* P)
        @test (P \ (   (1 .+ x) .* P))[1:10,1:10] ≈ (I + J)[1:10,1:10]

        x = Inclusion(0..1)
        Q = P[2x.-1,:]
        @test x .* Q isa MulQuasiArray
        @test Q \ (x .* Q) isa LazyBandedMatrices.Tridiagonal
    end

    @testset "sum" begin
        P = Legendre()
        x = axes(P,1)
        w = P * (P \ exp.(x))
        @test sum(w) ≈ ℯ - inv(ℯ)
    end

    @testset "Mapped" begin
        P = legendre(0..1)
        x = axes(P,1)
        X = jacobimatrix(P)
        @test X[1:10,1:10] ≈ (P \ (x .* P))[1:10,1:10]
        @test X[band(1)][1:10] ≈ [X[k,k+1] for k=1:10]
        @test X[band(0)][1:10] ≈ [X[k,k] for k=1:10]
        @test X[band(-1)][1:10] ≈ [X[k+1,k] for k=1:10]
        A,B,C = recurrencecoefficients(P)
        @test P[0.1,1] == 1
        @test P[0.1,2] ≈ A[1]*0.1 + B[1]
        @test P[0.1,3] ≈ (A[2]*0.1 + B[2])*P[0.1,2] - C[2]
        @test P[0.1,4] ≈ (A[3]*0.1 + B[3])*P[0.1,3] - C[3]*P[0.1,2]
    end

    @testset "Christoffel–Darboux" begin
        P = Legendre()
        X = P\ (axes(P,1) .* P)
        Mi = inv(P'P)
        x,y = 0.1,0.2
        n = 10
        Pn = Diagonal([Ones(n); Zeros(∞)])
        Min = Pn * Mi
        @test (X*Min - Min*X')[1:n,1:n] ≈ zeros(n,n)
        @test (x-y) * P[x,1:n]'Mi[1:n,1:n]*P[y,1:n] ≈ P[x,n:n+1]' * (X*Min - Min*X')[n:n+1,n:n+1] * P[y,n:n+1]
        β = X[n,n+1]*Mi[n+1,n+1]
        @test (x-y) * P[x,1:n]'Mi[1:n,1:n]*P[y,1:n] ≈ P[x,n:n+1]' * [0 -β; β 0] * P[y,n:n+1]
    end

    @testset "special syntax" begin
        @test legendrep.(0:5, 0.3) == Legendre()[0.3, 1:6]
    end
end


text:
"""
    StatsBase.predict(m::LinearMixedModel, newdata;
                    new_re_levels=:missing)
    StatsBase.predict(m::GeneralizedLinearMixedModel, newdata;
                    new_re_levels=:missing, type=:response)

Predict response for new data.

!!! note
    Currently, no in-place methods are provided because these methods
    internally construct a new model and therefore allocate not just a
    response vector but also many other matrices.

!!! warning
    Models are assumed to be full rank.

!!! warning
    These methods construct an entire MixedModel behind the scenes and
    as such may use a large amount of memory when `newdata` is large.

The keyword argument `new_re_levels` specifies how previously unobserved
values of the grouping variable are handled. Possible values are
`:population` (return population values, i.e. fixed-effects only),
`:missing` (return `missing`), `:error` (error on this condition;
the error type is an implementation detail). If you want simulated values
for unobserved levels of the grouping variable, consider the
[`simulate!`](@ref) and `simulate` methods.

Predictions based purely on the fixed effects can be obtained by
specifying previously unobserved levels of the random effects and setting
`new_re_levels=:population`. In the future, it may be possible to specify
a subset of the grouping variables or overall random-effects structure to use,
but not at this time.

For `GeneralizedLinearMixedModel`, the `type` parameter specifies
whether the predictions should be returned on the scale of linear predictor
(`:linpred`) or on the response scale (`:response`). If you don't know the
difference between these terms, then you probably want `type=:response`.

Regression weights are not yet supported in prediction.
Similarly, offsets are also not supported for `GeneralizedLinearMixedModel`.
"""
function StatsBase.predict(
    m::LinearMixedModel, newdata::Tables.ColumnTable; new_re_levels=:population
)
    return _predict(m, newdata, m.β; new_re_levels)
end

function StatsBase.predict(
    m::GeneralizedLinearMixedModel,
    newdata::Tables.ColumnTable;
    new_re_levels=:population,
    type=:response,
)
    type in (:linpred, :response) || throw(ArgumentError("Invalid value for type: $(type)"))

    y = _predict(m.LMM, newdata, m.β; new_re_levels)

    return type == :linpred ? y : broadcast!(Base.Fix1(linkinv, Link(m)), y, y)
end

# β is separated out here because m.β != m.LMM.β depending on how β is estimated for GLMM
function _predict(m::MixedModel{T}, newdata, β; new_re_levels) where {T}
    new_re_levels in (:population, :missing, :error) ||
        throw(ArgumentError("Invalid value for new_re_levels: $(new_re_levels)"))

    # if we ever support simulation, here some old bits from the docstrings
    # `new_re_levels=:simulate` is also not yet available for `GeneralizedLinearMixedModel`.
    # , `:simulate` (simulate new values).
    # For `:simulate`, the values are determined by solving for their values by
    # using the existing model's estimates for the new data. (These are in general
    # *not* the same values as the estimates computed on the new data.)

    # the easiest thing here is to just assemble a new model and
    # pass that to the other predict methods....
    # this can probably be made much more efficient

    # note that the contrasts don't matter for prediction purposes
    # (at least for the response)

    # add a response column
    # we get type stability via constant propogation on `new_re_levels`
    y, mnew = let ytemp = ones(T, length(first(newdata)))
        f, contr = _abstractify_grouping(m.formula)
        lmm = LinearMixedModel(f, newdata; contrasts=contr)
        ytemp =
            new_re_levels == :missing ? convert(Vector{Union{T,Missing}}, ytemp) : ytemp

        ytemp, lmm
    end

    grps = fnames(m)
    mul!(y, mnew.X, β)
    # mnew.reterms for the correct Z matrices
    # ranef(m) for the BLUPs from the original fit

    # because the reterms are sorted during model construction by
    # number of levels and that number may not be the same for the
    # new data, we need to permute the reterms from both models to be
    # in the same order
    newreperm = sort(1:length(mnew.reterms); by=x -> string(mnew.reterms[x].trm))
    oldreperm = sort(1:length(m.reterms); by=x -> string(m.reterms[x].trm))
    newre = mnew.reterms[newreperm]
    oldre = m.reterms[oldreperm]

    if new_re_levels == :error
        for (grp, known_levels, data_levels) in
            zip(grps, levels.(m.reterms), levels.(mnew.reterms))
            if sort!(known_levels) != sort!(data_levels)
                throw(ArgumentError("New level enountered in $grp"))
            end
        end

        # we don't have to worry about the BLUP ordering within a given
        # grouping variable because we are in the :error branch
        blups = ranef(m)[oldreperm]
    elseif new_re_levels == :population
        blups = ranef(mnew)[newreperm]
        blupsold = ranef(m)[oldreperm]

        for (idx, B) in enumerate(blups)
            oldlevels = levels(oldre[idx])
            for (lidx, ll) in enumerate(levels(newre[idx]))
                oldloc = findfirst(isequal(ll), oldlevels)
                if oldloc === nothing
                    # setting a BLUP to zero gives you the population value
                    B[:, lidx] .= zero(T)
                else
                    B[:, lidx] .= @view blupsold[idx][:, oldloc]
                end
            end
        end
    elseif new_re_levels == :missing
        # we can't quite use ranef! because we need
        # Union{T, Missing} and not just T
        blups = Vector{Matrix{Union{T,Missing}}}(undef, length(m.reterms))
        copyto!(blups, ranef(mnew)[newreperm])
        blupsold = ranef(m)[oldreperm]
        for (idx, B) in enumerate(blups)
            oldlevels = levels(oldre[idx])
            for (lidx, ll) in enumerate(levels(newre[idx]))
                oldloc = findfirst(isequal(ll), oldlevels)
                if oldloc === nothing
                    # missing is poisonous so propogates
                    B[:, lidx] .= missing
                else
                    B[:, lidx] .= @view blupsold[idx][:, oldloc]
                end
            end
        end
        # elseif new_re_levels == :simulate
        #     @show m.θ
        #     updateL!(setθ!(mnew, m.θ))
        #     blups = ranef(mnew)[newreperm]
        #     blupsold = ranef(m)[oldreperm]
        #     for (idx, B) in enumerate(blups)
        #         oldlevels = levels(oldre[idx])
        #         for (lidx, ll) in enumerate(levels(newre[idx]))
        #             oldloc = findfirst(isequal(ll), oldlevels)
        #             if oldloc === nothing
        #                 # keep the new value
        #             else
        #                 B[:, lidx] = @view blupsold[idx][:, oldloc]
        #             end
        #         end
        #     end
    else
        throw(ErrorException("Impossible branch reached. Please report an issue on GitHub"))
    end

    for (rt, bb) in zip(newre, blups)
        mul!(y, rt, bb, one(T), one(T))
    end

    return y
end

# yup, I got lazy on this one -- let the dispatched method handle kwarg checking
function StatsBase.predict(m::MixedModel, newdata; kwargs...)
    return predict(m, columntable(newdata); kwargs...)
end

## should we add in code for prediction intervals?
# we don't directly implement (Wald) confidence intervals, so directly
# supporting (Wald) prediction intervals seems a step too far right now


text:
import FEM.getmat
import FEM.getvec

facts("FEM.Tools") do

context("FEM.Tools.FloatArrayPool") do

matpool = MatPool()
vecpool = VecPool()

A = getmat(2, 2, "J", matpool)
B = getmat(2, 2, "J", matpool)
C = getmat(2, 2, "K", matpool)
D = getmat(3, 2, "J", matpool)

@fact is(A, B) => true
@fact is(A, D) => false
@fact is(B, C) => false

a = getvec(4, "N", vecpool)
b = getvec(4, "N", vecpool)
c = getvec(6, "N", vecpool)
d = getvec(4, "K", vecpool)

@fact is(a, b) => true
@fact is(b, c) => false
@fact is(b, d) => false

function function_creating_matrix(matpool::MatPool)
    return getmat(2,2, "a", matpool)
end

f = getmat(2,2, "a", matpool)
v = function_creating_matrix(matpool)
@fact is(f, v) => true

end # context

end # facts


text:
using Test

using KalmanFilters.UnscentedKalmanFilter
using KalmanFilters.SigmaPoints

@testset "MerweScaled" begin
include("test_merwe.jl")
end

@testset "Unscented Kalman Filter" begin
include("test_ukf.jl")
end

max_stars_repo_name: pxl-th/KalmanFilters

max_issues_repo_name: pxl-th/KalmanFilters

max_forks_repo_name: pxl-th/KalmanFilters


text:
Base.isequal(F::C, G::C) where {C <: ComponentData} =
    all(f -> isequal(getfield(F, f), getfield(G, f)), 1:nfields(F))::Bool

Base.:(==)(F::C, G::C) where {C <: ComponentData} =
    all(f -> getfield(F, f)== getfield(G, f), 1:nfields(F))::Bool

@inline function Base.hash(c::C, h::UInt) where {C <: ComponentData}
    for f in nfields(c)
        h = hash(getfield(c, f), h)
    end
    return h
end

"Can be used to specify the type of component storage to be used for a given `ComponentData`."
component_type(::Type{<:ComponentData}) = Component

@inline indices_iterator(a::AbstractComponent) = a.indices

"""
The most basic Component type.

Indexing into a component with an `Entity` will return the data linked to that entity,
indexing with a regular `Int` will return directly the `ComponentData` that is stored in the data
vector at that index, i.e. generally not the storage linked to the `Entity` with that `Int` as id.
"""
struct Component{T <: ComponentData} <: AbstractComponent{T}
    indices::Indices
    data   ::Vector{T}
end

Component{T}() where {T} = Component(Indices(), T[])

##### BASE Extensions ####
Base.eltype(::Type{<:AbstractComponent{T}}) where T = T

Base.length(c::AbstractComponent) = length(c.data)

Base.in(i::Integer, c::AbstractComponent) = in(i, c.indices)
Base.in(e::AbstractEntity, c::AbstractComponent)  = in(e.id, c)

Base.isempty(c::AbstractComponent) = isempty(c.data)

function Base.delete!(c::AbstractComponent, es::Vector{Entity})
    for e in es
        if e in c
            pop!(c, e)
        end
    end
end

function Base.permute!(c::AbstractComponent, permvec)
    permute!(c.data, permvec)
    permute!(c.indices, permvec)
end

Base.@propagate_inbounds @inline Base.getindex(c::Component,       e::AbstractEntity)  = c.data[c.indices[e.id]]
Base.@propagate_inbounds @inline Base.getindex(c::Component,       i::Integer) = c.data[i]

@inline function Base.setindex!(c::Component{T}, v::T, e::AbstractEntity) where {T}
    eid = e.id
    @boundscheck if !in(e, c)
        push!(c.indices, eid)
        push!(c.data, v)
        return v
    end
    @inbounds c.data[c.indices[eid]] = v
    return v
end

#TODO This is necessary  to support || in @entities_in, do we need that?
# See below for the code without
@inline function Base.pointer(c::Component{T}, e::AbstractEntity) where {T}
    if e.id in c.indices
        if T.mutable
            unsafe_load(reinterpret(Ptr{Ptr{T}}, pointer(c.data, @inbounds c.indices[e.id])))
        else
            return pointer(c.data, @inbounds c.indices[e.id])
        end
    else
        return Ptr{T}()
    end
end

# @inline Base.pointer(c::Component{T}, e::Entity) where {T}=
#     pointer(c.data, @inbounds c.indices[e.id])

function Base.empty!(c::Component)
    empty!(c.indices)
    empty!(c.data)
    return c
end

function swap_order!(c::AbstractComponent, e1::AbstractEntity, e2::AbstractEntity)
    @boundscheck if !in(e1, c)
        throw(BoundsError(c, e1))
    elseif !in(e2, c)
        throw(BoundsError(c, e2))
    end
    @inbounds begin
        id1, id2 = swap_order!(c.indices, e1.id, e2.id)
        c.data[id1], c.data[id2] = c.data[id2], c.data[id1]
    end
end

function pop_indices_data!(c::AbstractComponent, e::AbstractEntity)
    @boundscheck if !in(e, c)
        throw(BoundsError(c, e))
    end
    n = length(c)
    @inbounds begin
        id = c.indices[e.id]
        v = c.data[id]
        c.data[id] = c.data[end]
        pop!(c.data)
        pop!(c.indices, e.id)
        return v
    end
end

Base.pop!(c::Component, e::AbstractEntity) = pop_indices_data!(c, e)

@inline Base.iterate(c::Component, args...) = iterate(c.data, args...)

function ensure_entity_id!(c::AbstractComponent, e::Int, id::Int)
    indices = c.indices
    @inbounds packed_id = indices[e]
    if packed_id != id
        @inbounds id_to_swap = indices.packed[id]
        swap_order!(indices, e, id_to_swap)
        c.data[id], c.data[packed_id] = c.data[packed_id], c.data[id]
    end
    return true
end

function shared_entity_ids(cs)
    l, id = findmin(map(length, cs))
    shortest = cs[id]
    shared_entity_ids = Int[]
    for (i, e) in enumerate(shortest.indices)
        if all(x->in(e, x.indices), cs)
            push!(shared_entity_ids, e)
        end
    end
    return shared_entity_ids
end

@inline function Base.hash(c::C, h::UInt) where {C <: AbstractComponent}
    for f in nfields(c)
        h = hash(getfield(c,f), h)
    end
    return h
end

########################################
#                                      #
#            Iteration                 #
#                                      #
########################################
struct EntityIterator{T <: Union{IndicesIterator,Indices,AbstractGroup}, TT <: Tuple}
    it::T
    components::TT
end

Base.IteratorSize(::EntityIterator) = Base.SizeUnknown()
Base.IteratorEltype(::EntityIterator) = Base.HasEltype()
Base.eltype(::EntityIterator{T, TT}) where {T, TT} = EntityState{TT}
Base.length(i::EntityIterator) = length(i.it)

Base.in(e::AbstractEntity, i::EntityIterator) = in(e.id, i.it)

struct EntityState{TT<:Tuple} <: AbstractEntity
    e::Entity
    components::TT
end

# TODO: Cleanup, can these two be merged?
@generated function Base.getproperty(e::EntityState{TT}, f::Symbol) where {TT}
    fn_to_DT = Dict{Symbol, DataType}()
    ex = :(getfield(e, f))
    ex = Expr(:elseif, :(f === :id), :(return getfield(getfield(e, :e), :id)::Int), ex)
    for PDT in TT.parameters
        DT = eltype(PDT)
        for (fn, ft) in zip(fieldnames(DT), fieldtypes(DT))
            if haskey(fn_to_DT, fn)
                fnq = QuoteNode(fn)
                DT_ = fn_to_DT[fn]
                ft_ = fieldtype(DT_, fn)
                ex = MacroTools.postwalk(ex) do x
                    if @capture(x, return getfield(e[$DT_], $fnq)::$ft_)
                        return quote
                            throw(error("Field $f found in multiple components in $e.\nPlease use entity_state[$($DT)].$f instead."))
                        end
                    else
                        return x
                    end
                end
            else
                fn_to_DT[fn] = DT
                fnq = QuoteNode(fn)
                ex = Expr(:elseif, :(f === $fnq), :(return getfield(e[$DT], $fnq)::$ft), ex)
            end
        end
    end
    ex.head = :if
    return quote
        $(Expr(:meta, :inline))
        $(Expr(:meta, :propagate_inbounds))
        $ex
    end
end

@generated function Base.setproperty!(e::EntityState{TT}, f::Symbol, val) where {TT}
    fn_to_DT = Dict{Symbol, DataType}()
    ex = :(setfield!(e, f))
    ex = Expr(:elseif, :(f === :id), :(return setfield!(getfield(e, :e), :id)::Int, val), ex)
    for PDT in TT.parameters
        DT = eltype(PDT)
        for (fn, ft) in zip(fieldnames(DT), fieldtypes(DT))
            if haskey(fn_to_DT, fn)
                fnq = QuoteNode(fn)
                DT_ = fn_to_DT[fn]
                ft_ = fieldtype(DT_, fn)
                ex = MacroTools.postwalk(ex) do x
                    if @capture(x, setfield!(e[$DT_], $fnq, val))
                        return quote
                            throw(error("Field $f found in multiple components in $e.\nPlease use entity_state[$($DT)].$f instead."))
                        end
                    else
                        return x
                    end
                end
            else
                fn_to_DT[fn] = DT
                fnq = QuoteNode(fn)
                ex = Expr(:elseif, :(f === $fnq), :(return setfield!(e[$DT], $fnq, val)), ex)
            end
        end
    end
    ex.head = :if
    return quote
        $(Expr(:meta, :inline))
        $(Expr(:meta, :propagate_inbounds))
        $ex
    end
end

@generated function component(e::EntityState{TT}, ::Type{T}) where {TT<:Tuple, T<:ComponentData}
    id = findfirst(x -> eltype(x) == T, TT.parameters)
    return quote
        $(Expr(:meta, :inline))
        getfield(e,:components)[$id]
    end
end

@inline Base.@propagate_inbounds Base.getindex(e::EntityState, ::Type{T}) where {T<:ComponentData} = component(e, T)[e.e]

@inline Base.@propagate_inbounds Base.setindex!(e::EntityState, x::T, ::Type{T}) where {T<:ComponentData} =
    component(e, T)[e] = x

@inline Base.@propagate_inbounds Base.length(::EntityState{TT}) where {TT} = length(TT.parameters)

@inline function Base.iterate(i::EntityIterator, state = 1)
    n = iterate(i.it, state)
    n === nothing && return n
    e = Entity(n[1])
    return EntityState(e, i.components), n[2]
end

Base.getindex(iterator::EntityIterator, i) = Entity(iterator.it.shortest.packed[i])

macro entities_in(indices_expr)
    expr, t_sets, t_orsets = expand_indices_bool(indices_expr)
    if length(t_sets) == 1 && isempty(t_orsets) && expr.args[2] isa Symbol
        return esc(:(Overseer.EntityIterator(Overseer.indices_iterator($(t_sets[1])), ($(t_sets[1]),))))
    else
        return esc(quote
            t_comps = $(Expr(:tuple, t_sets...))
            t_or_comps = $(Expr(:tuple, t_orsets...))
            sets = map(Overseer.indices_iterator, t_comps)
            orsets = map(Overseer.indices_iterator, t_or_comps)
            if isempty(sets)
                minlen, minid = findmin(map(length, orsets))
                t_shortest = orsets[minid]
            else
                minlen, minid = findmin(map(length, sets))
                t_shortest = sets[minid]
            end
            if $(!isempty(t_orsets))
                shortest = deepcopy(t_shortest)
                for s in orsets
                    union!(shortest, s)
                end
            else
                shortest = t_shortest
            end
            Overseer.EntityIterator(Overseer.IndicesIterator(shortest, x->$expr), (t_comps..., t_or_comps...,))
        end)
    end
end

macro entities_in(ledger, indices_expr)
    expr, t_sets, t_orsets = expand_indices_bool(indices_expr)
    t_comp_defs = quote
    end
    comp_sym_map = Dict()
    for s in [t_sets; t_orsets]
        sym = gensym()
        t_comp_defs = quote
            $t_comp_defs
            $sym = $ledger[$s]
        end
        comp_sym_map[s] = sym
    end
    t_comp_defs = MacroTools.rmlines(MacroTools.flatten(t_comp_defs))

    expr = MacroTools.postwalk(expr) do x
        if x in keys(comp_sym_map)
            return comp_sym_map[x]
        else
            return x
        end
    end
    t_sets = map(x -> comp_sym_map[x], t_sets)
    t_orsets = map(x -> comp_sym_map[x], t_orsets)

    if length(t_sets) == 1 && isempty(t_orsets) && expr.args[2] isa Symbol
        return esc(quote
            $t_comp_defs
            Overseer.EntityIterator(Overseer.indices_iterator($(t_sets[1])), ($(t_sets[1]),))
        end)
    else
        return esc(quote
            $t_comp_defs
            t_comps = $(Expr(:tuple, t_sets...))
            t_or_comps = $(Expr(:tuple, t_orsets...))
            sets = map(Overseer.indices_iterator, t_comps)
            orsets = map(Overseer.indices_iterator, t_or_comps)
            if isempty(sets)
                minlen, minid = findmin(map(length, orsets))
                t_shortest = orsets[minid]
            else
                minlen, minid = findmin(map(length, sets))
                t_shortest = sets[minid]
            end
            if $(!isempty(t_orsets))
                shortest = deepcopy(t_shortest)
                for s in orsets
                    union!(shortest, s)
                end
            else
                shortest = t_shortest
            end
            Overseer.EntityIterator(Overseer.IndicesIterator(shortest, x->$expr), (t_comps..., t_or_comps...,))
        end)
    end
end

function Base.:(==)(c1::C1, c2::C2) where {C1 <: AbstractComponent, C2 <: AbstractComponent}
    if eltype(C1) != eltype(C2) ||length(c1) != length(c2)
        return false
    elseif length(c1) > 20 && hash(c1) != hash(c2)
        return false
    else
        return all(e -> (e in c2) && (@inbounds c2[e] == c1[e]), @entities_in(c1))
    end
end

##############################
#                            #
#     Component Macros       #
#                            #
##############################

function process_typedef(typedef, mod)
    global td = nothing
    MacroTools.postwalk(typedef) do x
        if @capture(x, struct T_ fields__ end | mutable struct T_ fields__ end)
            global td = T
        end
        x
    end
    tn = MacroTools.namify(td)
    if @capture(td, T_ <: V_)
        if !Base.eval(mod, :($V <: Overseer.ComponentData))
            error("Components can only have supertypes which are subtypes of ComponentData.")
        else
            return typedef, tn
        end
    else
        typedef_ = MacroTools.postwalk(typedef) do x
            if MacroTools.isexpr(x) && x.head == :struct
                x.args[2] = :($(x.args[2]) <: Overseer.ComponentData)
            end
            x
        end
        return typedef_, tn
    end
end

macro component(typedef)
	return esc(Overseer._component(typedef, __module__))
end
function _component(typedef, mod)
    t = process_typedef(typedef, mod)
	t1, tn = t
	return quote
	    $t1
        Overseer.component_type(::Type{$tn}) = Overseer.Component
    end
end

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

struct ApplyToPool
    e::Entity
end

Base.parent(e::Entity) = ApplyToPool(e)

struct PooledComponent{T <: ComponentData} <: AbstractComponent{T}
    indices::Indices
    pool::Vector{Int}
    pool_size::Vector{Int}
    data::Vector{T}
end

PooledComponent{T}() where {T <: ComponentData} = PooledComponent{T}(Indices(), Int[], Int[], T[])

Base.@propagate_inbounds @inline Base.getindex(c::PooledComponent, e::AbstractEntity) = c.data[pool(c, e)]
Base.@propagate_inbounds @inline Base.getindex(c::PooledComponent, i::Integer) = c.data[c.pool[i]]
Base.@propagate_inbounds @inline pool(c::PooledComponent, e::AbstractEntity) = c.pool[c.indices[e.id]]
Base.@propagate_inbounds @inline pool(c::PooledComponent, e::Int) = c.pool[c.indices[e]]

Base.@propagate_inbounds @inline function Base.parent(c::PooledComponent, i::Int)
    @boundscheck if i > length(c.data)
        throw(BoundsError(c, i))
    end
    return @inbounds Entity(c.indices.packed[findfirst(isequal(i), c.pool)])
end
Base.@propagate_inbounds @inline Base.parent(c::PooledComponent, e::Entity) = parent(c, pool(c, e))

function is_unique_in(value, collection)
    count = 0
    for element in collection
        count += element == value
    end
    return count == 1
end

# c[entity] = value
# set value of <only> this entity
@inline function Base.setindex!(c::PooledComponent{T}, v::T, e::AbstractEntity) where {T}
    eid = e.id
    @inbounds if in(e, c)
        pid = c.indices[eid]
        g = c.pool[pid]
        if c.pool_size[g] == 1
            # the entity already has its own (otherwise empty) pool - adjust value
            c.data[g] = v
        else
            # the entity is part of a larger pool - create a new one
            c.pool_size[g] -= 1
            push!(c.data, v)
            push!(c.pool_size, 1)
            c.pool[pid] = length(c.data)
        end
    else
        # the entity is not in the component - add it
        push!(c.indices, eid)
        push!(c.pool, length(c.data)+1)
        push!(c.pool_size, 1)
        push!(c.data, v)
    end
    return v
end

# c[entity] = parent
# set the value of this entity to that of parent
@inline function Base.setindex!(c::PooledComponent, p::AbstractEntity, e::AbstractEntity)
    @boundscheck if !in(p, c)
        throw(BoundsError(c, p))
    end
    @inbounds begin
        pg = pool(c, p)
        if in(e, c)
            eg = pool(c, e)
            if c.pool_size[eg] == 1
                # if this entity is the only one holding onto a value, remove
                # that value and cleanup pool indices
                deleteat!(c.data, eg)
                deleteat!(c.pool_size, eg)
                for i in eachindex(c.pool)
                    c.pool[i] = c.pool[i] - (c.pool[i] > eg)
                end
            else
                c.pool_size[eg] -= 1
            end
            # adjust pool index either way
            c.pool[c.indices[e.id]] = pg
        else
            # if the entity is not in there we have to add it
            push!(c.indices, e.id)
            push!(c.pool, pg)
        end
        c.pool_size[pg] += 1

        return c[p]
    end
end

# c[ParentGroup(entity)] = value
# set the value for all entities pooled with entity
@inline function Base.setindex!(c::PooledComponent{T}, v::T, x::ApplyToPool) where {T}
    e = x.e
    @boundscheck if !in(e, c)
        throw(BoundsError(c, e))
    end
    @inbounds c.data[pool(c, e.id)] = v
    return v
end

Base.length(c::PooledComponent) = length(c.pool)

function Base.empty!(c::PooledComponent)
    empty!(c.indices)
    empty!(c.pool)
    empty!(c.pool_size)
    empty!(c.data)
    return c
end

function Base.pop!(c::PooledComponent, e::AbstractEntity)
    @boundscheck if !in(e, c)
        throw(BoundsError(c, e))
    end

    @inbounds begin
        id = c.indices[e.id]
        g = c.pool[id]

        c.pool[id] = c.pool[end]
        c.pool_size[g] -= 1
        pop!(c.pool)
        pop!(c.indices, e.id)

        val = c.data[g]

        if c.pool_size[g] == 0
            deleteat!(c.data, g)
            deleteat!(c.pool_size, g)
            for i in eachindex(c.pool)
                if c.pool[i] > g
                    c.pool[i] -= 1
                end
            end
        end
    end

    return val
end

@inline Base.iterate(c::PooledComponent, args...) = iterate(c.data, args...)
Base.sortperm(c::PooledComponent) = sortperm(c.pool)

macro pooled_component(typedef)
    return esc(Overseer._pooled_component(typedef, __module__))
end

function _pooled_component(typedef, mod)
    t = process_typedef(typedef, mod)
    t1, tn = t
    return quote
        $t1
       	Overseer.component_type(::Type{$tn}) = Overseer.PooledComponent
    end
end

function make_unique!(c::PooledComponent)
    @inbounds begin
        # Go through all pools and check if a later pool has the same value.
        # If it does mark the later pool empty.
        for i in eachindex(c.pool)
            g0 = c.pool[i]
            if c.pool_size[g0] > 0
                v0 = c.data[g0]
                for j in i+1:length(c.pool)
                    g = c.pool[j]
                    if c.pool_size[g] > 0 && c.data[g] == v0
                        c.pool_size[g] -= 1
                        c.pool_size[g0] += 1
                        c.pool[j] = g0
                    end
                end
            end
        end

        # Go through all pools again and replace empty pools (i.e. its data
        # and pool_size) by the current last pool. Adjust pool indices
        # accordingly.
        i = 1
        while i <= length(c.pool_size)
            if c.pool_size[i] == 0
                if i == length(c.pool_size)
                    pop!(c.pool_size)
                    pop!(c.data)
                    break
                else
                    N = length(c.pool_size)
                    c.pool_size[i] = pop!(c.pool_size)
                    c.data[i] = pop!(c.data)
                    for j in eachindex(c.pool)
                        c.pool[j] = c.pool[j] == N ? i : c.pool[j]
                    end
                end
            else
                i += 1
            end
        end
    end

    return
end

struct PooledEntityIterator{T}
    c::PooledComponent{T}
    pool_id::Int
end

indices_iterator(g::PooledEntityIterator) = g
indices(g::PooledEntityIterator) = g

function entity_pool(c::PooledComponent, pool_id::Int)
    @boundscheck if length(c.data) < pool_id
        throw(BoundsError(c, pool_id))
    end
    return PooledEntityIterator(c, pool_id)
end

@inline entity_index(it::PooledEntityIterator, i::Int) = @inbounds it.c.indices.packed[findnext(isequal(it.pool_id), it.c.pool, i)]

@inline function Base.iterate(i::PooledEntityIterator, state = 1)
    state > i.c.pool_size[i.pool_id] && return nothing
    n = findnext(isequal(i.pool_id), i.c.pool, state)
    n === nothing && return n
    return Entity(i.c.indices.packed[n]), n + 1
end

Base.getindex(iterator::PooledEntityIterator, i) = iterate(iterator, i)[1]

Base.IteratorSize(::Type{PooledEntityIterator}) = Base.HasLength()
Base.IteratorEltype(::Type{PooledEntityIterator}) = Base.HasEltype()
Base.eltype(::PooledEntityIterator) = Entity
Base.eltype(::Type{PooledEntityIterator{T}}) where {T} = T
Base.length(i::PooledEntityIterator) = i.c.pool_size[i.pool_id]

Base.@propagate_inbounds @inline Base.in(e::Entity, it::PooledEntityIterator) =
    in(e, it.c) && @inbounds pool(it.c, e) == it.pool_id
Base.@propagate_inbounds @inline Base.in(e::Int, it::PooledEntityIterator) =
    in(e, it.c.indices) && @inbounds pool(it.c, e) == it.pool_id


text:
@testset "35.search-insert-position.jl" begin
    @test search_insert([1, 3, 5, 6], 5) == 2
    @test search_insert([1, 3, 5, 6], 2) == 1
    @test search_insert([1, 3, 5, 6], 7) == 4
    @test search_insert([1, 3, 5, 6], 0) == 0
end


text:
using GridArrays, DomainSets, StaticArrays, Test
@testset "boundingbox" begin
    @test boundingbox(0.0..1.0) == 0.0..1.0
    @test boundingbox((0.0..1.0)^3) == (0.0..1.0)^3
    @test boundingbox(UnitDisk())≈ChebyshevInterval()^2
    @test boundingbox(UnitBall())≈ChebyshevInterval()^3

    @test boundingbox(3UnitDisk())≈(-3.0..3.0)^2
    @test boundingbox(3UnitBall())≈(-3.0..3.0)^3
    @test boundingbox(3*((0.0..1.0)^3))≈ (0.0..3.0)^3
    @test boundingbox(3*(-1.0..2.0))≈ (-3.0..6.0)

    @test boundingbox(3UnitDisk().+SVector(.4,.5))≈(-2.6..3.4)×(-2.5..3.5)
    @test boundingbox(3UnitBall().+SVector(-.4,.5,0.))≈(-3.4..2.6)×(-2.5..3.5)×(-3.0..3.0)

    uniondomain = UnionDomain(3UnitDisk() .+ SVector(.4,.5)  , 3UnitDisk() .+ SVector(-.4,.5))
    @test boundingbox(uniondomain)≈(-3.4..3.4)×(-2.5..3.5)

    intersectdomain = IntersectDomain(3UnitDisk() .+ SVector(.4,.5)  , 3UnitDisk() .+ SVector(-.4,.5))
    @test boundingbox(intersectdomain)≈(-2.6..2.6)×(-2.5..3.5)

    diffdomain = SetdiffDomain(3UnitDisk() .+ SVector(.4,.5)  , 3UnitDisk() .+ SVector(-.4,.5))
    @test boundingbox(diffdomain)≈(-2.6..3.4)×(-2.5..3.5)
end


text:
@static if VERSION < v"0.7-"
    srand(12345)
else
    Random.seed!(12345)
end

@testset "gemm!" begin
    for elty in eltypes
        # A_mul_B!
        A = rand(elty, m_L3, k_L3)
        A_cl = cl.CLArray(queue, A)
        B = rand(elty, k_L3, n_L3)
        B_cl = cl.CLArray(queue, B)
        C = rand(elty, m_L3, n_L3)
        C_cl = cl.CLArray(queue, C)
        α = rand(elty)
        β = rand(elty)

        CLBlast.gemm!('N', 'N', α, A_cl, B_cl, β, C_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.gemm!('N', 'N', α, A, B, β, C)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(B_cl, queue=queue) ≈ B
        @test cl.to_host(C_cl, queue=queue) ≈ C

        @test_throws DimensionMismatch CLBlast.gemm!('T', 'N', α, A_cl, B_cl, β, C_cl, queue=queue)
        @test_throws DimensionMismatch CLBlast.gemm!('C', 'N', α, A_cl, B_cl, β, C_cl, queue=queue)
        @test_throws DimensionMismatch CLBlast.gemm!('N', 'T', α, A_cl, B_cl, β, C_cl, queue=queue)
        @test_throws DimensionMismatch CLBlast.gemm!('N', 'C', α, A_cl, B_cl, β, C_cl, queue=queue)

        @test_throws ArgumentError CLBlast.gemm!('A', 'N', α, A_cl, B_cl, β, C_cl, queue=queue)
        @test_throws ArgumentError CLBlast.gemm!('N', 'A', α, A_cl, B_cl, β, C_cl, queue=queue)

        # At_mul_B!, Ac_mul_B!
        A = rand(elty, k_L3, m_L3)
        A_cl = cl.CLArray(queue, A)
        B = rand(elty, k_L3, n_L3)
        B_cl = cl.CLArray(queue, B)
        C = rand(elty, m_L3, n_L3)
        C_cl = cl.CLArray(queue, C)
        α = rand(elty)
        β = rand(elty)

        CLBlast.gemm!('T', 'N', α, A_cl, B_cl, β, C_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.gemm!('T', 'N', α, A, B, β, C)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(B_cl, queue=queue) ≈ B
        @test cl.to_host(C_cl, queue=queue) ≈ C

        CLBlast.gemm!('C', 'N', α, A_cl, B_cl, β, C_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.gemm!('C', 'N', α, A, B, β, C)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(B_cl, queue=queue) ≈ B
        @test cl.to_host(C_cl, queue=queue) ≈ C

        @test_throws DimensionMismatch CLBlast.gemm!('N', 'N', α, A_cl, B_cl, β, C_cl, queue=queue)
        @test_throws DimensionMismatch CLBlast.gemm!('N', 'T', α, A_cl, B_cl, β, C_cl, queue=queue)
        @test_throws DimensionMismatch CLBlast.gemm!('N', 'C', α, A_cl, B_cl, β, C_cl, queue=queue)

        # A_mul_Bt!, A_mul_Bc!
        A = rand(elty, m_L3, k_L3)
        A_cl = cl.CLArray(queue, A)
        B = rand(elty, n_L3, k_L3)
        B_cl = cl.CLArray(queue, B)
        C = rand(elty, m_L3, n_L3)
        C_cl = cl.CLArray(queue, C)
        α = rand(elty)
        β = rand(elty)

        CLBlast.gemm!('N', 'T', α, A_cl, B_cl, β, C_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.gemm!('N', 'T', α, A, B, β, C)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(B_cl, queue=queue) ≈ B
        @test cl.to_host(C_cl, queue=queue) ≈ C

        CLBlast.gemm!('N', 'C', α, A_cl, B_cl, β, C_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.gemm!('N', 'C', α, A, B, β, C)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(B_cl, queue=queue) ≈ B
        @test cl.to_host(C_cl, queue=queue) ≈ C

        @test_throws DimensionMismatch CLBlast.gemm!('N', 'N', α, A_cl, B_cl, β, C_cl, queue=queue)
        @test_throws DimensionMismatch CLBlast.gemm!('T', 'N', α, A_cl, B_cl, β, C_cl, queue=queue)
        @test_throws DimensionMismatch CLBlast.gemm!('C', 'N', α, A_cl, B_cl, β, C_cl, queue=queue)

        # At_mul_Bt!, Ac_mul_Bt!, At_mul_Bt!, Ac_mul_Bc!
        A = rand(elty, k_L3, m_L3)
        A_cl = cl.CLArray(queue, A)
        B = rand(elty, n_L3, k_L3)
        B_cl = cl.CLArray(queue, B)
        C = rand(elty, m_L3, n_L3)
        C_cl = cl.CLArray(queue, C)
        α = rand(elty)
        β = rand(elty)

        CLBlast.gemm!('T', 'T', α, A_cl, B_cl, β, C_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.gemm!('T', 'T', α, A, B, β, C)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(B_cl, queue=queue) ≈ B
        @test cl.to_host(C_cl, queue=queue) ≈ C

        CLBlast.gemm!('C', 'T', α, A_cl, B_cl, β, C_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.gemm!('C', 'T', α, A, B, β, C)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(B_cl, queue=queue) ≈ B
        @test cl.to_host(C_cl, queue=queue) ≈ C

        CLBlast.gemm!('T', 'C', α, A_cl, B_cl, β, C_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.gemm!('T', 'C', α, A, B, β, C)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(B_cl, queue=queue) ≈ B
        @test cl.to_host(C_cl, queue=queue) ≈ C

        CLBlast.gemm!('C', 'C', α, A_cl, B_cl, β, C_cl, queue=queue)
        Compat.LinearAlgebra.BLAS.gemm!('C', 'C', α, A, B, β, C)
        @test cl.to_host(A_cl, queue=queue) ≈ A
        @test cl.to_host(B_cl, queue=queue) ≈ B
        @test cl.to_host(C_cl, queue=queue) ≈ C

        @test_throws DimensionMismatch CLBlast.gemm!('N', 'N', α, A_cl, B_cl, β, C_cl, queue=queue)
        @test_throws DimensionMismatch CLBlast.gemm!('T', 'N', α, A_cl, B_cl, β, C_cl, queue=queue)
        @test_throws DimensionMismatch CLBlast.gemm!('C', 'N', α, A_cl, B_cl, β, C_cl, queue=queue)
        @test_throws DimensionMismatch CLBlast.gemm!('N', 'T', α, A_cl, B_cl, β, C_cl, queue=queue)
        @test_throws DimensionMismatch CLBlast.gemm!('N', 'C', α, A_cl, B_cl, β, C_cl, queue=queue)
    end
end

@testset "symm!" begin
    for elty in eltypes
        # multiply from the left
        A = rand(elty, m_L3, m_L3)
        A_cl = cl.CLArray(queue, A)
        B = rand(elty, m_L3, n_L3)
        B_cl = cl.CLArray(queue, B)
        C = rand(elty, m_L3, n_L3)
        C_cl = cl.CLArray(queue, C)
        α = rand(elty)
        β = rand(elty)

        for uplo in ['U','L']
            CLBlast.symm!('L', uplo, α, A_cl, B_cl, β, C_cl, queue=queue)
            Compat.LinearAlgebra.BLAS.symm!('L', uplo, α, A, B, β, C)
            @test cl.to_host(A_cl, queue=queue) ≈ A
            @test cl.to_host(B_cl, queue=queue) ≈ B
            @test cl.to_host(C_cl, queue=queue) ≈ C

            @test_throws DimensionMismatch CLBlast.symm!('L', uplo, α, B_cl, A_cl, β, C_cl, queue=queue)
            @test_throws DimensionMismatch CLBlast.symm!('R', uplo, α, A_cl, B_cl, β, C_cl, queue=queue)
        end

        @test_throws ArgumentError CLBlast.symm!('A', 'U', α, A_cl, B_cl, β, C_cl, queue=queue)
        @test_throws ArgumentError CLBlast.symm!('U', 'A', α, A_cl, B_cl, β, C_cl, queue=queue)

        # multiply from the right
        A = rand(elty, n_L3, n_L3)
        A_cl = cl.CLArray(queue, A)
        B = rand(elty, m_L3, n_L3)
        B_cl = cl.CLArray(queue, B)
        C = rand(elty, m_L3, n_L3)
        C_cl = cl.CLArray(queue, C)
        α = rand(elty)
        β = rand(elty)

        for uplo in ['U','L']
            CLBlast.symm!('R', uplo, α, A_cl, B_cl, β, C_cl, queue=queue)
            Compat.LinearAlgebra.BLAS.symm!('R', uplo, α, A, B, β, C)
            @test cl.to_host(A_cl, queue=queue) ≈ A
            @test cl.to_host(B_cl, queue=queue) ≈ B
            @test cl.to_host(C_cl, queue=queue) ≈ C

            @test_throws DimensionMismatch CLBlast.symm!('L', uplo, α, A_cl, B_cl, β, C_cl, queue=queue)
            @test_throws DimensionMismatch CLBlast.symm!('R', uplo, α, B_cl, A_cl, β, C_cl, queue=queue)
        end
    end
end

@testset "hemm!" begin
    for elty in eltypes
        elty <: Complex || continue

        # multiply from the left
        A = rand(elty, m_L3, m_L3)
        A_cl = cl.CLArray(queue, A)
        B = rand(elty, m_L3, n_L3)
        B_cl = cl.CLArray(queue, B)
        C = rand(elty, m_L3, n_L3)
        C_cl = cl.CLArray(queue, C)
        α = rand(elty)
        β = rand(elty)

        for uplo in ['U','L']
            CLBlast.hemm!('L', uplo, α, A_cl, B_cl, β, C_cl, queue=queue)
            Compat.LinearAlgebra.BLAS.hemm!('L', uplo, α, A, B, β, C)
            @test cl.to_host(A_cl, queue=queue) ≈ A
            @test cl.to_host(B_cl, queue=queue) ≈ B
            @test cl.to_host(C_cl, queue=queue) ≈ C

            @test_throws DimensionMismatch CLBlast.hemm!('L', uplo, α, B_cl, A_cl, β, C_cl, queue=queue)
            @test_throws DimensionMismatch CLBlast.hemm!('R', uplo, α, A_cl, B_cl, β, C_cl, queue=queue)
        end

        @test_throws ArgumentError CLBlast.hemm!('A', 'U', α, A_cl, B_cl, β, C_cl, queue=queue)
        @test_throws ArgumentError CLBlast.hemm!('U', 'A', α, A_cl, B_cl, β, C_cl, queue=queue)

        # multiply from the right
        A = rand(elty, n_L3, n_L3)
        A_cl = cl.CLArray(queue, A)
        B = rand(elty, m_L3, n_L3)
        B_cl = cl.CLArray(queue, B)
        C = rand(elty, m_L3, n_L3)
        C_cl = cl.CLArray(queue, C)
        α = rand(elty)
        β = rand(elty)

        for uplo in ['U','L']
            CLBlast.hemm!('R', uplo, α, A_cl, B_cl, β, C_cl, queue=queue)
            Compat.LinearAlgebra.BLAS.hemm!('R', uplo, α, A, B, β, C)
            @test cl.to_host(A_cl, queue=queue) ≈ A
            @test cl.to_host(B_cl, queue=queue) ≈ B
            @test cl.to_host(C_cl, queue=queue) ≈ C

            @test_throws DimensionMismatch CLBlast.hemm!('L', uplo, α, A_cl, B_cl, β, C_cl, queue=queue)
            @test_throws DimensionMismatch CLBlast.hemm!('R', uplo, α, B_cl, A_cl, β, C_cl, queue=queue)
        end
    end
end

@testset "syrk!" begin
    for elty in eltypes
        # A*A'
        A = rand(elty, n_L3, k_L3)
        A_cl = cl.CLArray(queue, A)
        C = rand(elty, n_L3, n_L3)
        C_cl = cl.CLArray(queue, C)
        α = rand(elty)
        β = rand(elty)

        for uplo in ['U','L']
            CLBlast.syrk!(uplo, 'N', α, A_cl, β, C_cl, queue=queue)
            Compat.LinearAlgebra.BLAS.syrk!(uplo, 'N', α, A, β, C)
            @test cl.to_host(A_cl, queue=queue) ≈ A
            @test cl.to_host(C_cl, queue=queue) ≈ C

            @test_throws DimensionMismatch CLBlast.syrk!(uplo, 'T', α, A_cl, β, C_cl, queue=queue)
        end

        @test_throws ArgumentError CLBlast.syrk!('A', 'N', α, A_cl, β, C_cl, queue=queue)
        @test_throws ArgumentError CLBlast.syrk!('U', 'A', α, A_cl, β, C_cl, queue=queue)

        # A'*A
        A = rand(elty, k_L3, n_L3)
        A_cl = cl.CLArray(queue, A)
        C = rand(elty, n_L3, n_L3)
        C_cl = cl.CLArray(queue, C)
        α = rand(elty)
        β = rand(elty)

        for uplo in ['U','L']
            CLBlast.syrk!(uplo, 'T', α, A_cl, β, C_cl, queue=queue)
            Compat.LinearAlgebra.BLAS.syrk!(uplo, 'T', α, A, β, C)
            @test cl.to_host(A_cl, queue=queue) ≈ A
            @test cl.to_host(C_cl, queue=queue) ≈ C

            @test_throws DimensionMismatch CLBlast.syrk!(uplo, 'N', α, A_cl, β, C_cl, queue=queue)
        end

        @test_throws ArgumentError CLBlast.syrk!('U', 'C', α, A_cl, β, C_cl, queue=queue)
    end
end

@testset "herk!" begin
    for elty in eltypes
        elty <: Complex || continue

        # A*A'
        A = rand(elty, n_L3, k_L3)
        A_cl = cl.CLArray(queue, A)
        C = rand(elty, n_L3, n_L3)
        for i in 1:n_L3
            C[i,i] = real(C[i,i])
        end
        C_cl = cl.CLArray(queue, C)
        α = real(rand(elty))
        β = real(rand(elty))

        for uplo in ['U','L']
            CLBlast.herk!(uplo, 'N', α, A_cl, β, C_cl, queue=queue)
            Compat.LinearAlgebra.BLAS.herk!(uplo, 'N', α, A, β, C)
            @test cl.to_host(A_cl, queue=queue) ≈ A
            @test_skip cl.to_host(C_cl, queue=queue) ≈ C

            @test_throws DimensionMismatch CLBlast.herk!(uplo, 'C', α, A_cl, β, C_cl, queue=queue)
        end

        @test_throws ArgumentError CLBlast.herk!('A', 'N', α, A_cl, β, C_cl, queue=queue)
        @test_throws ArgumentError CLBlast.herk!('U', 'A', α, A_cl, β, C_cl, queue=queue)

        # A'*A
        A = rand(elty, k_L3, n_L3)
        A_cl = cl.CLArray(queue, A)
        C = rand(elty, n_L3, n_L3)
        for i in 1:n_L3
            C[i,i] = real(C[i,i])
        end
        C_cl = cl.CLArray(queue, C)
        α = real(rand(elty))
        β = real(rand(elty))

        for uplo in ['U','L']
            CLBlast.herk!(uplo, 'C', α, A_cl, β, C_cl, queue=queue)
            Compat.LinearAlgebra.BLAS.herk!(uplo, 'C', α, A, β, C)
            @test cl.to_host(A_cl, queue=queue) ≈ A
            @test_skip cl.to_host(C_cl, queue=queue) ≈ C

            @test_throws DimensionMismatch CLBlast.herk!(uplo, 'N', α, A_cl, β, C_cl, queue=queue)
        end

        @test_throws ArgumentError CLBlast.herk!('U', 'T', α, A_cl, β, C_cl, queue=queue)
    end
end

@testset "syr2k!" begin
    for elty in eltypes
        # A*B'
        A = rand(elty, n_L3, k_L3)
        A_cl = cl.CLArray(queue, A)
        B = rand(elty, n_L3, k_L3)
        B_cl = cl.CLArray(queue, B)
        C = rand(elty, n_L3, n_L3)
        C_cl = cl.CLArray(queue, C)
        α = rand(elty)
        β = rand(elty)

        for uplo in ['U','L']
            CLBlast.syr2k!(uplo, 'N', α, A_cl, B_cl, β, C_cl, queue=queue)
            Compat.LinearAlgebra.BLAS.syr2k!(uplo, 'N', α, A, B, β, C)
            @test cl.to_host(A_cl, queue=queue) ≈ A
            @test cl.to_host(B_cl, queue=queue) ≈ B
            @test cl.to_host(C_cl, queue=queue) ≈ C

            @test_throws DimensionMismatch CLBlast.syr2k!(uplo, 'T', α, A_cl, B_cl, β, C_cl, queue=queue)
        end

        @test_throws ArgumentError CLBlast.syr2k!('A', 'N', α, A_cl, B_cl, β, C_cl, queue=queue)
        @test_throws ArgumentError CLBlast.syr2k!('U', 'A', α, A_cl, B_cl, β, C_cl, queue=queue)

        # A'*B
        A = rand(elty, k_L3, n_L3)
        A_cl = cl.CLArray(queue, A)
        B = rand(elty, k_L3, n_L3)
        B_cl = cl.CLArray(queue, B)
        C = rand(elty, n_L3, n_L3)
        C_cl = cl.CLArray(queue, C)
        α = rand(elty)
        β = rand(elty)

        for uplo in ['U','L']
            CLBlast.syr2k!(uplo, 'T', α, A_cl, B_cl, β, C_cl, queue=queue)
            Compat.LinearAlgebra.BLAS.syr2k!(uplo, 'T', α, A, B, β, C)
            @test cl.to_host(A_cl, queue=queue) ≈ A
            @test cl.to_host(B_cl, queue=queue) ≈ B
            @test cl.to_host(C_cl, queue=queue) ≈ C

            @test_throws DimensionMismatch CLBlast.syr2k!(uplo, 'N', α, A_cl, B_cl, β, C_cl, queue=queue)
        end
    end
end

@testset "her2k!" begin
    for elty in eltypes
        elty <: Complex || continue

        # A*B'
        A = rand(elty, n_L3, k_L3)
        A_cl = cl.CLArray(queue, A)
        B = rand(elty, n_L3, k_L3)
        B_cl = cl.CLArray(queue, B)
        C = rand(elty, n_L3, n_L3)
        C_cl = cl.CLArray(queue, C)
        for i in 1:n_L3
            C[i,i] = real(C[i,i])
        end
        α = rand(elty)
        β = real(rand(elty))

        for uplo in ['U','L']
            CLBlast.her2k!(uplo, 'N', α, A_cl, B_cl, β, C_cl, queue=queue)
            Compat.LinearAlgebra.BLAS.her2k!(uplo, 'N', α, A, B, β, C)
            @test cl.to_host(A_cl, queue=queue) ≈ A
            @test cl.to_host(B_cl, queue=queue) ≈ B
            @test_skip cl.to_host(C_cl, queue=queue) ≈ C

            @test_throws DimensionMismatch CLBlast.her2k!(uplo, 'C', α, A_cl, B_cl, β, C_cl, queue=queue)
        end

        @test_throws ArgumentError CLBlast.her2k!('A', 'N', α, A_cl, B_cl, β, C_cl, queue=queue)
        @test_throws ArgumentError CLBlast.her2k!('U', 'A', α, A_cl, B_cl, β, C_cl, queue=queue)

        # A'*B
        A = rand(elty, k_L3, n_L3)
        A_cl = cl.CLArray(queue, A)
        B = rand(elty, k_L3, n_L3)
        B_cl = cl.CLArray(queue, B)
        C = rand(elty, n_L3, n_L3)
        for i in 1:n_L3
            C[i,i] = real(C[i,i])
        end
        C_cl = cl.CLArray(queue, C)
        α = rand(elty)
        β = real(rand(elty))

        for uplo in ['U','L']
            CLBlast.her2k!(uplo, 'C', α, A_cl, B_cl, β, C_cl, queue=queue)
            Compat.LinearAlgebra.BLAS.her2k!(uplo, 'C', α, A, B, β, C)
            @test cl.to_host(A_cl, queue=queue) ≈ A
            @test cl.to_host(B_cl, queue=queue) ≈ B
            @test_skip cl.to_host(C_cl, queue=queue) ≈ C

            @test_throws DimensionMismatch CLBlast.her2k!(uplo, 'N', α, A_cl, B_cl, β, C_cl, queue=queue)
        end

        @test_throws ArgumentError CLBlast.her2k!('U', 'T', α, A_cl, B_cl, β, C_cl, queue=queue)
    end
end

@testset "trmm!" begin
    for elty in eltypes
        # multiply from the left
        A = rand(elty, m_L3, m_L3)
        A_cl = cl.CLArray(queue, A)
        B = rand(elty, m_L3, n_L3)
        B_cl = cl.CLArray(queue, B)
        α = rand(elty)

        for uplo in ['U','L'], transA in ['N','T','C'], diag in ['N','U']
            CLBlast.trmm!('L', uplo, transA, diag, α, A_cl, B_cl, queue=queue)
            Compat.LinearAlgebra.BLAS.trmm!('L', uplo, transA, diag, α, A, B)
            @test cl.to_host(A_cl, queue=queue) ≈ A
            @test cl.to_host(B_cl, queue=queue) ≈ B

            @test_throws DimensionMismatch CLBlast.trmm!('L', uplo, transA, diag, α, B_cl, A_cl, queue=queue)
            @test_throws DimensionMismatch CLBlast.trmm!('R', uplo, transA, diag, α, A_cl, B_cl, queue=queue)
        end

        @test_throws ArgumentError CLBlast.trmm!('A', 'U', 'N', 'N', α, A_cl, B_cl, queue=queue)
        @test_throws ArgumentError CLBlast.trmm!('U', 'A', 'N', 'N', α, A_cl, B_cl, queue=queue)

        # multiply from the right
        A = rand(elty, n_L3, n_L3)
        A_cl = cl.CLArray(queue, A)
        B = rand(elty, m_L3, n_L3)
        B_cl = cl.CLArray(queue, B)
        α = rand(elty)
        β = rand(elty)

        for uplo in ['U','L'], transA in ['N','T','C'], diag in ['N','U']
            CLBlast.trmm!('R', uplo, transA, diag, α, A_cl, B_cl, queue=queue)
            Compat.LinearAlgebra.BLAS.trmm!('R', uplo, transA, diag, α, A, B)
            @test cl.to_host(A_cl, queue=queue) ≈ A
            @test cl.to_host(B_cl, queue=queue) ≈ B

            @test_throws DimensionMismatch CLBlast.trmm!('L', uplo, transA, diag, α, A_cl, B_cl, queue=queue)
            @test_throws DimensionMismatch CLBlast.trmm!('R', uplo, transA, diag, α, B_cl, A_cl, queue=queue)
        end
    end
end

@testset "trsm!" begin
    for elty in eltypes
        # On Travis, there is some strange error
        # https://travis-ci.org/JuliaGPU/CLBlast.jl/jobs/414395954#L347
        @static if VERSION < v"0.7-"
            is_apple() && break
        else
            Sys.isapple() && break
        end

        # multiply from the left
        A = rand(elty, m_L3, m_L3)
        for i in 1:m_L3
            A[i,i] = i
        end
        A_cl = cl.CLArray(queue, A)
        B = rand(elty, m_L3, n_L3)
        B_cl = cl.CLArray(queue, B)
        α = rand(elty)

        for uplo in ['U','L'], transA in ['N','T','C'], diag in ['N','U']
            CLBlast.trsm!('L', uplo, transA, diag, α, A_cl, B_cl, queue=queue)
            Compat.LinearAlgebra.BLAS.trsm!('L', uplo, transA, diag, α, A, B)
            @test cl.to_host(A_cl, queue=queue) ≈ A
            @test cl.to_host(B_cl, queue=queue) ≈ B

            @test_throws DimensionMismatch CLBlast.trsm!('L', uplo, transA, diag, α, B_cl, A_cl, queue=queue)
            @test_throws DimensionMismatch CLBlast.trsm!('R', uplo, transA, diag, α, A_cl, B_cl, queue=queue)
        end

        @test_throws ArgumentError CLBlast.trsm!('A', 'U', 'N', 'N', α, A_cl, B_cl, queue=queue)
        @test_throws ArgumentError CLBlast.trsm!('U', 'A', 'N', 'N', α, A_cl, B_cl, queue=queue)

        # multiply from the right
        A = rand(elty, n_L3, n_L3)
        for i in 1:n_L3
            A[i,i] = i
        end
        A_cl = cl.CLArray(queue, A)
        B = rand(elty, m_L3, n_L3)
        B_cl = cl.CLArray(queue, B)
        α = rand(elty)
        β = rand(elty)

        for uplo in ['U','L'], transA in ['N','T','C'], diag in ['N','U']
            CLBlast.trsm!('R', uplo, transA, diag, α, A_cl, B_cl, queue=queue)
            Compat.LinearAlgebra.BLAS.trsm!('R', uplo, transA, diag, α, A, B)
            @test cl.to_host(A_cl, queue=queue) ≈ A
            @test cl.to_host(B_cl, queue=queue) ≈ B

            @test_throws DimensionMismatch CLBlast.trsm!('L', uplo, transA, diag, α, A_cl, B_cl, queue=queue)
            @test_throws DimensionMismatch CLBlast.trsm!('R', uplo, transA, diag, α, B_cl, A_cl, queue=queue)
        end
    end
end


text:
input = read(joinpath(@__DIR__, "input"), String)
coords, instructions = split(input, "\n\n")

paper = Set()
for coord in split(coords, "\n")
    x, y = parse.(Int, match(r"(\d+),(\d+)", coord).captures)
    push!(paper, (x, y))
end

folds = []
for i in split(instructions, "\n", keepempty=false)
    dim, n = match(r".+([xy])=(\d+)$", i).captures
    push!(folds, (dim, parse(Int, n)))
end

function translatex(coord, n)
    x, y = coord
    x > n && return (n - (x - n), y)
    return (x, y)
end

function translatey(coord, n)
    x, y = coord
    y > n && return (x, n - (y - n))
    return (x, y)
end

function fold(paper, folds)
    paper = copy(paper)
    for (dim, n) in folds
        if dim == "x"
            paper = map(idx -> translatex(idx, n), collect(paper))
            paper = filter(idx -> first(idx) <= n, paper)
            paper = unique(paper)
        else
            paper = map(idx -> translatey(idx, n), collect(paper))
            paper = filter(idx -> last(idx) <= n, paper)
            paper = unique(paper)
        end
    end
    return paper
end

firstfold = first(folds)
p1 = length(fold(paper, [firstfold]))

paper = fold(paper, folds)
x = maximum(map(idx -> first(idx), collect(paper)))
y = maximum(map(idx -> last(idx), collect(paper)))

println("-----------------------------------------------------------------------")
println("transparent origami -- part one :: $p1")
println("transparent origami -- part two ::")
println()
for j in 0:y
    print("    ")
    for i in 0:x
        (i, j) in paper ? print("█") : print(" ")
    end
    println()
end
println()
println("-----------------------------------------------------------------------")

@assert(p1 == 759)


text:
export @D, @DD, @GRAD, @pot

# ===========================================================================
#     implement some fun little macros for easier access
#     to the potentials
# ===========================================================================

# REMARK on @pot:
# ---------------
# Julia 0.4 version:
# call(pp::Potential, varargs...) = evaluate(pp, varargs...)
# call(pp::Potential, ::Type{Val{:D}}, varargs...) = evaluate_d(pp, varargs...)
# call(pp::Potential, ::Type{Val{:DD}}, varargs...) = evaluate_dd(pp, varargs...)
# call(pp::Potential, ::Type{Val{:GRAD}}, varargs...) = grad(pp, varargs...)
# unfortunately, in 0.5 `call` doesn't take anabstractargument anymore,
# which means that we need to specify for every potential how to
# create this syntactic sugar. This is what `@pot` is for.

"""
Annotate a type with `@pot` to setup the syntax sugar
for `evaluate, evaluate_d, evaluate_dd, grad`.

## Usage:

For example, the declaration
```julia
"documentation for `LennardJones`"
mutable struct LennardJones <: PairPotential
   r0::Float64
end
@pot LennardJones
```
creates the following aliases:
```julia
lj = LennardJones(1.0)
lj(args...) = evaluate(lj, args...)
@D lj(args...) = evaluate_d(lj, args...)
@DD lj(args...) = evaluate_dd(lj, args...)
@GRAD lj(args...) = grad(lj, args...)
```

Usage of `@pot` is not restricted to pair potentials, but can be applied to
*any* type.
"""
macro pot(fsig)
   @assert fsig isa Symbol
   sym = esc(:x)
   tsym = esc(fsig)
   quote
      @inline ($sym::$tsym)(args...) = evaluate($sym, args...)
      @inline ($sym::$tsym)(::Type{Val{:D}}, args...) = evaluate_d($sym, args...)
      @inline ($sym::$tsym)(::Type{Val{:DD}}, args...) = evaluate_dd($sym, args...)
      @inline ($sym::$tsym)(::Type{Val{:GRAD}}, args...) = grad($sym, args...)
   end
end

# --------------------------------------------------------------------------

# next create macros that translate
"""
`@D`: Use to evaluate the derivative of a potential. E.g., to compute the
Lennard-Jones potential,
```julia
lj = LennardJones()
r = 1.0 + rand(10)
ϕ = lj(r)
ϕ' = @D lj(r)
```
see also `@DD`.
"""
macro D(fsig::Expr)
    @assert fsig.head == :call
    insert!(fsig.args, 2, Val{:D})
    for n = 1:length(fsig.args)
        fsig.args[n] = esc(fsig.args[n])
    end
    return fsig
end

"`@DD` : analogous to `@D`"
macro DD(fsig::Expr)
    @assert fsig.head == :call
    for n = 1:length(fsig.args)
        fsig.args[n] = esc(fsig.args[n])
    end
    insert!(fsig.args, 2, Val{:DD})
    return fsig
end

"`@GRAD` : analogous to `@D`, but escapes to `grad`"
macro GRAD(fsig::Expr)
    @assert fsig.head == :call
    for n = 1:length(fsig.args)
        fsig.args[n] = esc(fsig.args[n])
    end
    insert!(fsig.args, 2, Val{:GRAD})
    return fsig
end


text:
module ArgParse2

using OrderedCollections: LittleDict
using TextWrap: println_wrapped, wrap

export ArgumentParser, add_argument!, argument_adder, parse_args, show_help, show_usage
export add_argument  # TODO: Remove add_argument before registering

const Optional{T} = Union{T,Nothing}

include("utils.jl")
include("argument.jl")
include("variable.jl")
include("argument_parser.jl")
include("help.jl")
include("parse.jl")

include("precompile.jl")

_precompile_()

end # module


text:
# this section includes:

# 1. methods to do simulation
# 2.

using Dierckx  # use 1D B-spline (linear) interpolations to approximate gradients, Dierckx is the julia wrapper of Fortran library Dierckx
using PyPlot

# ----------------------------------
# SECTION: quick methods
# ----------------------------------
# 1. expr to refresh env & workspace data
macro macro_RefreshEnvData( PopulationCSVFile::String )
   # usage: @macro_RefreshEnvData( "./data/Demography_base.csv" )
   return :(
   env = (
      T = 400, S = 80, Sr = 40, START_AGE = 20, START_YEAR = 1945,
      PATH_DEMOGCSV = $PopulationCSVFile,
      PATH_WAGEPROFILE = "./data/WageProfileCoef.csv", PATH_MA2MB = "./data/MA2MBCoef.csv",
      PATH_M2C = "./data/M2C.csv", PATH_TFPGROW = "./data/tfpGrowthProfile.csv",
   );
   include("$(pwd())\\src\\proc_VarsDeclare.jl"); # use pwd() to avoid possible ambiguity of working directory
   include("$(pwd())\\src\\proc_InitPars.jl");
   )
end # macro_RefreshEnvData

# 2. expr to solve a model
expr_RunModel = :(begin;
   include("$(pwd())\\src\\proc_VarsDeclare.jl"); # use pwd() to avoid possible ambiguity of working directory
   include("$(pwd())\\src\\proc_InitPars.jl");
   println("solving......")
   Guess = ( r = 0.08, L = 0.2 );EasySearch.SteadyState!( 1, Guess, Dt, Dst, Pt, Ps, Pc, env, atol = 1E-8, MaxIter = 1000, PrintMode = "silent", MagicNum = 2.0, StepLen = 0.5 ); Guess = ( r = 0.12, L = 0.75 );
   EasySearch.SteadyState!( env.T, Guess, Dt, Dst, Pt, Ps, Pc, env, atol = 1E-6, MaxIter = 1000, PrintMode = "silent", MagicNum = 2.0, StepLen = 0.5 );
   PerfLog = EasySearch.Transition!( Dt, Dst, Pt, Ps, Pc, env, atol = 1E-3, MaxIter = 500, PrintMode = "silent", MagicNum = 2.0, StepLen = 0.7, ReturnLog = true );
   EasySearch.ProcAfterTransition!( Dt, Dst, Pt, Ps, Pc, env );
   end;
)

# 3. small function to generate a two-stage-flat cross-sectional demography whose population size is 1
"""
   afunc_FlatDemog( Nretire2Nrate::Float64 ,S::Int, Sr::Int )

生成一个given max age S, given retirement age Sr, given working population/total population rate Nwork2Nrate,
的cross-sectional demography (a vector with length S);
# NOTE: 固定总人口为1，工作期flat，退休期flat，调整两期总人口比例

returns a vector with length S.
"""
afunc_FlatDemog( Nretire2Nrate::Float64 ,S::Int, Sr::Int ) = begin
   @assert( Nretire2Nrate <= 0.501, "non stable demography!" ) # allows a small ex to compute gradient (well, strictly we should use left-side gradient but ... to save time)
   return cat( fill( (1.0 - Nretire2Nrate) / Sr, Sr), fill( Nretire2Nrate / (S - Sr) , S - Sr)  , dims = 1 )::Vector{Float64}
end # afunc_FlatDemog

# 4. method to get the collection which saves required results
# NOTE: also returns c̄, w̄, l̄, Ã, B̃, LI/PoolExp, LI/PoolIn, LI/Y, LI/TaxRev in the same tuple;
#       used to compute contributions, GE effects etc.
# NOTE: EasySearch.ProcAfterTransition! is not essential
function afunc_getRes( t, Dt, Dst, Pt, Ps, Pc, env )
   # compute temp vars
   local AggMB = sum( Ps[:N][t,1:env.S] .* Dst[:MB][t,1:env.S] )
   local AggPoolExp = AggMB * ( 1.0 - Pt[:cpB][t] )
   local AggPoolIn = AggPoolExp - Dt[:LI][t]
   local NinSS = 1.0
   # mean q (well ... in fact, they are equal among ages)
   local qmean = sum(Ps[:q][t,:]) / env.S

   # collect
   tmpRes = (
      LIpercapita = Dt[:LI][t] / NinSS,  # gap per capita
      N = NinSS,  # total population
      ρ = 1.0 - sum(Ps[:N][t,1:env.Sr]) / NinSS  , # the retired population ratio
      c̄ = Dt[:C][t] / NinSS ,
      w̄ = Dt[:w̄][t],
      l̄ = Dt[:L][t] / sum(Ps[:N][t,1:env.Sr]),  # labor per capita is calculated on working people
      Ã = qmean * (1.0 - Pt[:cpB][t]) / ( 1.0 + sum(Ps[:p]) / env.S ) ,  # multiplier of benefit policy
      B̃ = (1.0 - Pt[:𝕒][t] - Pt[:𝕓][t]) * Pt[:ζ][t] / (1.0 + Pt[:z][t] * Pt[:η][t] + Pt[:ζ][t] ) ,  # multiplier of contribution policy
      gap2exp = Dt[:LI][t] / AggPoolExp,
      gap2in = Dt[:LI][t] / AggPoolIn,
      gap2gdp = Dt[:LI][t] / Dt[:Y][t],
      gap2taxrev = Dt[:LI][t] / (Dt[:TRw][t] + Dt[:TRc][t]),
      AggPoolExp = AggPoolExp,
      AggPoolIn = AggPoolIn,
      Y = Dt[:Y][t],
      TaxRev = Dt[:TRw][t] + Dt[:TRc][t],
      TRw = Dt[:TRw][t],
      TRc = Dt[:TRw][t],
      L = Dt[:L][t],
      K  =Dt[:K][t],
      G = Dt[:G][t],
      r = Dt[:r][t],
      C = Dt[:C][t],
      I = Dt[:I][t],
   )
   return tmpRes::NamedTuple
end # afunc_getRes

# 6. funtions from ρ to: 𝕃, c̄, w̄, l̄; and their gradient at specific ρ level
# NOTE: ρ is defined as the retired population share
# NOTE: the following functions are defined IN ORDER, they refer defined functions
# NOTE: recommend to select the whole 6. then execute
# NOTE: we use Interpolations.gradients to compute gradient which avoids (visible) floating errors
# --------------------
   function solveSSwithRho( ρ::Real; t::Int = 2010 - 1945 + 1 , Dt = Dt, Dst = Dst, Pt = Pt, Ps = Ps, Pc = Pc, env = env )
      # create two-stage flat demography
      Ps[:N][t,:] = afunc_FlatDemog( ρ, env.S, env.Sr )
      # Ps[:N][t,:] = Ps[:N][end,:]  # natural exponential mortality case
      # do not forget to adjust the survival rate (more important than population size)
      for s in 1:(env.S-1)
         Ps[:F][t,s] = 1.0 - Ps[:N][t,s+1] / Ps[:N][t,s]
      end; Ps[:F][t,env.S] = 0.0
      # use p̄ to replace p_s
      Ps[:p][:] .= sum(Ps[:p]) / env.S
      # search SS at that year
      EasySearch.SteadyState!( t, ( r = 0.08, L = 0.2 ), Dt, Dst, Pt, Ps, Pc, env, atol = 1E-8, MaxIter = 1000, PrintMode = "silent", MagicNum = 2.0, StepLen = 0.5 ); Guess = ( r = 0.12, L = 0.75 );
      # collect data and return
      return afunc_getRes( t, Dt, Dst, Pt, Ps, Pc, env )
   end
   f_𝕃percapita(ρ::Float64 ; Year::Int = 2010 ) = solveSSwithRho(ρ, t = Year - 1945 + 1 ).LIpercapita
   f_c̄(ρ::Float64 ; Year::Int = 2010 ) = solveSSwithRho(ρ, t = Year - 1945 + 1 ).c̄
   f_w̄(ρ::Float64 ; Year::Int = 2010 ) = solveSSwithRho(ρ, t = Year - 1945 + 1 ).w̄
   f_l̄(ρ::Float64 ; Year::Int = 2010 ) = solveSSwithRho(ρ, t = Year - 1945 + 1 ).l̄

# solve dense-enough curves(e.g. ρ -> 𝕃) and do interpolations; they are used to evaluate gradients
# NOTE: use the linear B-spline of Dierckx package
# NOTE: itp_ vars are used to plot 0-order; grad_ means 1st-order
grid_ρ = Array(LinRange( 0.08, 0.5, 300 ) ) # ρ lattice
grid_Res = [ solveSSwithRho(ρ) for ρ in grid_ρ ]

# interpolations & gradient
   # 1) gap per capita (𝕃)
   itp_𝕃percapita = Spline1D(grid_ρ,[ grid_Res[x].LIpercapita for x in 1:length(grid_Res) ])
   # 2) the 1st-order derivatives of 𝕃 (i.e. d𝕃/dρ = 𝔾 + 𝔻)
   grad_𝕃percapita = derivative(itp_𝕃percapita, grid_ρ)
   # 3) consumption per capita (c̄)
   itp_c̄ = Spline1D(grid_ρ, [ grid_Res[x].c̄ for x in 1:length(grid_Res) ])
   # 4) dc̄/dρ
   grad_c̄ = derivative(itp_c̄, grid_ρ)
   # 5) average wage (w̄)
   itp_w̄ = Spline1D(grid_ρ, [ grid_Res[x].w̄ for x in 1:length(grid_Res) ])
   # 6) dw̄/dρ
   grad_w̄ = derivative(itp_w̄, grid_ρ)
   # 7) labor supply per capita (l̄)
   itp_l̄ = Spline1D(grid_ρ, [ grid_Res[x].l̄ for x in 1:length(grid_Res) ])
   # 8) dl̄/dρ
   grad_l̄ = derivative(itp_l̄, grid_ρ)

# extract vectors used later
   itp_Ã = Spline1D(grid_ρ, [ grid_Res[x].Ã for x in 1:length(grid_Res) ])
   itp_B̃ = Spline1D(grid_ρ, [ grid_Res[x].B̃ for x in 1:length(grid_Res) ])
   itp_gap2exp = Spline1D(grid_ρ, [ grid_Res[x].gap2exp for x in 1:length(grid_Res) ])
   itp_gap2in = Spline1D(grid_ρ, [ grid_Res[x].gap2in for x in 1:length(grid_Res) ])
   itp_gap2gdp = Spline1D(grid_ρ, [ grid_Res[x].gap2gdp for x in 1:length(grid_Res) ])
   itp_gap2taxrev = Spline1D(grid_ρ, [ grid_Res[x].gap2taxrev for x in 1:length(grid_Res) ])

# compute lattice vectors of GE effect 𝔾, and cross effect 𝔻
   # 1) GE effect 𝔾
   itp_𝔾 = itp_Ã(grid_ρ) .* grad_c̄ .+ itp_B̃(grid_ρ) .* itp_w̄(grid_ρ) .* itp_l̄(grid_ρ)
   itp_𝔾 = Spline1D(grid_ρ,itp_𝔾) # to unify denotations
   # 2) cross effect 𝔻
   itp_𝔻 = (grid_ρ .- 1.0) .* itp_B̃(grid_ρ) .* ( grad_w̄ .* itp_l̄(grid_ρ) + grad_l̄ .* itp_w̄(grid_ρ) )
   itp_𝔻 = Spline1D(grid_ρ,itp_𝔻)

# 7. reset datasets under benchmark fertility/demography projection
#    then tell & record the real(projected) ρ in 2010
@macro_RefreshEnvData( "./data/Demography_base.csv" ) # refresh it!
tmp = Ps[:N][2010 - 1945 + 1, :]
realρ2010 = 1.0 - sum(tmp[1:env.Sr]) / sum(tmp)
println("The real/projected 2010's rho is: ", realρ2010 )

# ----------------------------------
# SECTION: relation: ρ -> LI/N (ρ is retired population ratio)
# NOTE: collects gap per capita, N, ρ, c̄, w̄, l̄, Ã, B̃, LI/PoolExp, LI/PoolIn, LI/Y, LI/TaxRev in the same tuple
# NOTE: using two-stage flat cross-sectional demography
# NOTE: ρ <= 0.5 to ensure it is a stable demography
# ----------------------------------
tmpexpr = :(
   grid(true); xlabel(L"$\rho$");
)
figure(figsize=(12,8))
   subplot(3,1,1)  # lines: 𝔾, 𝔻
      plot(grid_ρ,itp_𝔾(grid_ρ),"-.r")
      plot(grid_ρ,itp_𝔻(grid_ρ),"--b")
      axvline(realρ2010, linestyle = ":" )  # add a vertical line indicating real ρ in 2010
      text( realρ2010, 0.003, L"Real 2010 ρ = " * string(realρ2010) )
      eval(tmpexpr)
      legend([
              L"GE effect ($\mathbb{G}$)",
              L"Cross term ($\mathbb{D}$)"],
              loc="center left", fontsize = 14)
  subplot(3,1,2)  # lines: 𝔾 = Ã(dc̄/dρ) + B̃w̄l̄
     plot(grid_ρ, grad_c̄ .+ itp_B̃(grid_ρ) ,"-.r")
     plot(grid_ρ, itp_B̃(grid_ρ) .* itp_w̄(grid_ρ) .* itp_l̄(grid_ρ) ,"--g")
     axvline(realρ2010, linestyle = ":" )  # add a vertical line indicating real ρ in 2010
     eval(tmpexpr)
     legend([
             L"Medical demand influence: $\tilde{A}(d\bar{c}/d\bar{\rho})$",
             L"Income level influence: $\tilde{B}\bar{w}\bar{l}$"],
             loc="center left", fontsize = 14)
   subplot(3,1,3)  # lines: 𝔻 = GE part * ρ
      plot(grid_ρ,grid_ρ .- 1.0,"-.r")
      plot(grid_ρ,itp_𝔻(grid_ρ) ./ (grid_ρ .- 1.0),"--g")
      axvline(realρ2010, linestyle = ":" )  # add a vertical line indicating real ρ in 2010
      eval(tmpexpr)
      legend([
              L"Linear form of $\rho$: $\rho-1$",
              L"Sensitivity to incomes: $ \tilde{B} ( \frac{\partial \bar{w}}{\partial \rho} \bar{l} + \frac{\partial \bar{l}}{\partial \rho} \bar{w} )  $  "],
              loc="center left", fontsize = 14)
   tight_layout()
   # save figure
   savefig("$(pwd())/output/Channel_DecomposeInSState.pdf", format = "pdf")

# ----------------------------------
# SECTION: extract other economic variables & interpolations
# NOTE: please run the programs above to get grid_Res and grid_ρ
# ----------------------------------
tmpEconVars = (
   itp_Y = Spline1D(grid_ρ,[ grid_Res[x].Y for x in 1:length(grid_Res) ]),
   itp_K = Spline1D(grid_ρ,[ grid_Res[x].K for x in 1:length(grid_Res) ]),
   itp_L = Spline1D(grid_ρ,[ grid_Res[x].L for x in 1:length(grid_Res) ]),
   itp_G = Spline1D(grid_ρ,[ grid_Res[x].G for x in 1:length(grid_Res) ]),
   itp_C = Spline1D(grid_ρ,[ grid_Res[x].C for x in 1:length(grid_Res) ]),
   itp_I = Spline1D(grid_ρ,[ grid_Res[x].I for x in 1:length(grid_Res) ]),
   itp_TRw = Spline1D(grid_ρ,[ grid_Res[x].TRw for x in 1:length(grid_Res) ]),
   itp_TRc = Spline1D(grid_ρ,[ grid_Res[x].TRc for x in 1:length(grid_Res) ]),
   itp_r = Spline1D(grid_ρ,[ grid_Res[x].r for x in 1:length(grid_Res) ]),
)

#


text:
"""
    st_as_stars(data::AbstractArray{T, 2}, mask::AbstractArray{Bool, 2})

    st_as_stars(data::AbstractArray{T, 2}, r_mask::GeoArray; outfile = nothing)
    st_as_stars(data::AbstractArray{T, 2}, mask::AbstractArray{Bool, 2},
        b::bbox; kw...)
    st_as_stars(data::AbstractArray{T, 2}, mask::AbstractArray{Bool, 2},
        b::AbstractArray{<:Real, 1}; kw...)
"""
function st_as_stars(data::AbstractArray{T, 2}, mask::AbstractArray{Bool, 2}) where T <: Real
    ntime = size(data, 2)
    arr = zeros(T, size(mask)..., ntime)

    ind = findall(mask)
    @views for i = 1:ntime
        x = arr[:, :, i]
        x[ind] = data[:, i]
    end
    arr
end

function st_as_stars(data::AbstractArray{T, 2}, r_mask::AbstractGeoArray; outfile = nothing) where T <: Real

    arr = st_as_stars(data, r_mask.A[:, :, 1])
    @show size(arr)

    r = GeoArray(r_mask, vals = arr)
    if outfile !== nothing; st_write(r, outfile); end
    r
end

function st_as_stars(data::AbstractArray{T, 2}, mask::AbstractArray{Bool, 2},
    b::bbox; kw...) where T <: Real

    st_as_stars(data, rast(mask, b); kw...)
end

function st_as_stars(data::AbstractArray{T, 2}, mask::AbstractArray{Bool, 2},
    b::AbstractArray{<:Real, 1}; kw...) where T <: Real

    st_as_stars(data, rast(mask, bbox(b...)); kw...)
end

export st_as_stars