text:
#! ~/bin/julia

# gen_hclust.jl

#=
    written circa Julia version 1.5
    If starting from nothing....

using Pkg;
Pkg.add("Clustering")
Pkg.add("StatsPlots")

*** can also enter the builtin package manager from the REPL with ']'
*** then for example just write `add DelimitedFiles`
=#

using DelimitedFiles
using Clustering
# using StatsPlots  # too expensive on a headless server

include("hclust2newick.jl")

function main(ARGS)

    # read in the sparse delimited representation of a distance matrix
    sdm_dir = ARGS[1]
    sdm_file = ARGS[2]

    sparse = readdlm(string(sdm_dir, "/", sdm_file))

    # find max index for observations  (in this case anyway)
    obs = maximum(sparse[1:size(sparse,1), 2]);

    # init a distance matrix
    dist = zeros(obs+1, obs+1); # experiment: +1 for zeroa in outgroup?

    # populate the upper and lower triangles with Jaccard distances
    for i in 1:size(sparse,1)
        dist[sparse[i,1],sparse[i,2]] = dist[sparse[i,2],sparse[i,1]] = sparse[i,7]
    end

    # save a copy of the symmetric matrix
    open(sdm_dir * "/hct-jaccard-dist-mat-sym.tsv", "w") do io
        writedlm(io, dist)
    end

    #=
        hclust struct info gleened

     propertynames(hc)
     (:merges, :heights, :order, :linkage, :height, :labels, :merge, :method)

    linkage === method one of

        :single (the default): use the minimum distance between any of the cluster members
        :average: use the mean distance between any of the cluster members
        :complete: use the maximum distance between any of the members
        :ward

    uplo:
        upper :U or lower :L   if unspecified expects symmetric

        ---------------
    branchorder (order?) one of:
        :r (the default as R does): based on the node heights & original elements order
        :barjoseph (or :optimal):

	merge === merges

	height ===  heights

    =#

    # generate the hierarchical cluster
    hc = hclust(dist, linkage=:complete, uplo=:L, branchorder=:optimal)

    # lookey lou (from repl any way)
    # plot(hc, xticks=false)

    #=
        oooo. grab a screen shot

        looks plausible. no clue how to export in useful format
        answer: have to make one!
        write scripts/hclust2newick.jl
    =#

    open(sdm_dir * "/hcdm.tree", "w") do io
        write_newick(io, hc)
    end

    # partition the observations into a set of clusters(features)
    # appropriate to the number of  cultivars (observations) (10:1 per folklore)
    # make the h-clust assignments

    hca = cutree(hc; k=Integer(ceil(size(hc.order)[1]/10)))

    # make the implicit index explicit
    idx_hca = [ia for ia in zip([1:1:length(hca);], hca)]

    # save a copy of the cluster assignments
    open(sdm_dir * "/iss_cultivar_cluster_assignment", "w") do io
        writedlm(io, idx_hca)
    end

end

main(ARGS)


text:
### A Pluto.jl notebook ###
# v0.14.8

using Markdown
using InteractiveUtils

# ╔═╡ ec4e1740-0b51-11ec-3fa2-d3c1d67f4572
using Revise

# ╔═╡ 9e67b5cd-f084-4bd4-9eb0-b94311184197
using CSV

# ╔═╡ b2a2fdb9-949a-498c-ac19-dab106582558
using DataFrames

# ╔═╡ 72553f5c-276d-4256-8354-04154017d346
using Plots

# ╔═╡ 67b93cb9-6335-448e-a09e-ead44db9d5c3
using LaTeXStrings

# ╔═╡ d20b171e-1bd3-44aa-b067-7bf1a3c0227e
Ω = real.(parse.(ComplexF64, (CSV.File("conductivity_data_6.csv") |> Tables.matrix)[3:end, 1]))

# ╔═╡ 8c909d72-7f28-49ed-9f6f-6b7f77e4e984
β = real.(parse.(ComplexF64, (CSV.File("conductivity_data_6.csv") |> Tables.matrix)[1, 2:end-1]))

# ╔═╡ 2980129c-4be6-4928-8bfe-b6b13eb0d50a
T = 1 ./ β

# ╔═╡ 29d4a2e3-692b-48b3-9ce0-fca83b55e8d3
v = real.((CSV.File("v_data.csv") |> Tables.matrix)[2:end-1, 2:end])

# ╔═╡ eab11453-2f25-493a-84e4-307231496653
v_0 = real.((CSV.File("v_data.csv") |> Tables.matrix)[end, 2:end])

# ╔═╡ 136a409d-9ea0-40ab-998a-f5bf954315ca
α = 1:12

# ╔═╡ a1bef2ba-0fb0-4a43-a09d-4e2e5ab6e069
w = real.((CSV.File("w_data.csv") |> Tables.matrix)[2:end-1, 2:end])

# ╔═╡ 671bed0e-082a-4b3a-9467-7be0ee2327df
w_0 = real.((CSV.File("w_data.csv") |> Tables.matrix)[end, 2:end])

# ╔═╡ a7431fd9-588a-41f1-b3d8-ef15837ccf75
F = real.((CSV.File("F_data.csv") |> Tables.matrix)[2:end-1, 2:end])

# ╔═╡ f7746913-7afa-4bce-b261-9532bd6f27cb
F_0 = real.((CSV.File("F_data.csv") |> Tables.matrix)[end, 2:end])

# ╔═╡ 780dd0f5-b922-4691-94b5-fa9881381c07
M = real.((CSV.File("M_data.csv") |> Tables.matrix)[2:end-1, 2:end])

# ╔═╡ d3f88ec2-6f32-4429-a34b-4db33fc55433
M_0 = real.((CSV.File("M_data.csv") |> Tables.matrix)[end, 2:end])

# ╔═╡ 39cff408-72c6-42ea-8172-d9d0fb9e94fc
σ = [parse.(ComplexF64, (CSV.File("conductivity_data_$i.csv") |> Tables.matrix)[3:end, 2:end-1]) for i in α]

# ╔═╡ 8e432e19-984e-4727-ab9c-25894e2585d8
σ_0 = hcat([parse.(ComplexF64, (CSV.File("conductivity_data_$i.csv") |> Tables.matrix)[3:end, end]) for i in α]...)

# ╔═╡ 83972098-4232-437e-9d2d-7620c17672f0
σ_0_dc = hcat([parse.(ComplexF64, (CSV.File("conductivity_data_$i.csv") |> Tables.matrix)[2, end]) for i in α]...)

# ╔═╡ d8762608-09eb-4ea0-9a84-5f46673a619d
σ_dc = hcat([parse.(ComplexF64, (CSV.File("conductivity_data_$i.csv") |> Tables.matrix)[2, 2:end-1]) for i in α]...)

# ╔═╡ 9ab20c3a-a0f5-4938-9339-720630fbd48b
dc_σ_real = Plots.contourf(α, reverse(T), reverse(log.(real.(σ_dc)), dims = 1), xlabel = "α", ylabel = "T / ω", fill = cgrad(:thermal, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, xticks = (1:12, hcat(["$i" for i in 1:12]...)), yticks = ([0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0], hcat([0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0]...)), yaxis = :log, right_margin = 4Plots.mm)

# ╔═╡ ed78a2eb-3313-4cae-b0d3-2118f5dff33f
savefig(dc_σ_real, "dc_mobility.pdf")

# ╔═╡ f524e982-6208-44bc-86df-ccfa851b0c92
athermal_theory = plot(α, [-F_0, M_0, v_0, w_0], yaxis = :log, legend = :topleft, linewidth = 3, linestyle = :auto, labels = [L"-\textrm{E_{gs}}" L"\textrm{m_p}" L"\textrm{v}" L"\textrm{w}" L"\textrm{\mu_0 (10^{6})}"], xlabel = "α", xticks = (1:12, ["$i" for i in 1:12]), ylabel = "Athermal Theory", tickfontsize = 12, labelfontsize = 12, legendfontsize = 12, size = (550, 550), yticks = ([1, 2, 4, 8, 16, 32, 64, 128, 256], ["1" "2" "4" "8" "16" "32" "64" "128" "256"]), ylims = (0.8, 256))

# ╔═╡ 9f09e259-014e-48b0-bd00-f7e2c9c897b3
plot(T, real.(σ_dc), ylims = (0.001, 10), linestyle = :dash, marker = :circle, xaxis = :log, yaxis = :log, legend = false)

# ╔═╡ ccda1394-45b4-4a9a-8c14-425cb0b81a87
log(v[end, 9])

# ╔═╡ fadad1b3-6da6-4724-be97-469d81641d8f
plot(α, real.(σ_dc)'[:, 1:10:end], linestyle = :dash, marker = :circle, legend = :outerright, yaxis = :log, label = hcat(["$(T[i])" for i in 1:10:length(T)]...))

# ╔═╡ 3918264a-f72b-4b4c-b67f-cd6fa4c6c846
savefig(athermal_theory, "athermal_theory.pdf")

# ╔═╡ c1d56910-a645-40ed-a5d8-0f7d5be19ae7
plot(Ω, hcat([real.(σ_0[:, i]) ./ maximum(real.(σ_0[:, i])) for i in 1:12]...), labels = hcat(["α = $i" for i in 1:12]...), linewidth = 2, legend = :outerright, xlabel = "Ω / ω", ylabel = "σ(Ω) [arb]", minorgrid= true)

# ╔═╡ 19065839-85a3-4768-830e-26f64bf4a278
hcat([real.(σ_0[:, i]) ./ maximum(real.(σ_0[:, i])) for i in 1:12]...)

# ╔═╡ 7c978142-5319-45f8-b854-efac81271fcc
function ridgeline(x, y, z; shift = 1, jump = 10, ylabel = :none, xlabel = :none, palette = :twilight, linewidth = 1.5, color = :black, size = (595, 842), ymirror = false, fillalpha = 1, step = 1)
	if !ymirror
		p = plot(x, z[:, 1], fillrange = z[:, 2 * jump] .- shift, fillalpha = fillalpha, ylabel = ylabel, xlabel = xlabel, xticks = (0:maximum(x)/10:maximum(x), string.(x[1:Int(floor((length(x)-1)/10)):end])), yticks = (z[1, 1:jump * step:end] .- shift .* range(0, stop = length(y[1:jump:end]) - 1, step = step), string.(y[1:jump * step:end])), legend = false, size = size, grid = false, ymirror = ymirror, palette = palette, tickfontsize = 12, labelfontsize = 12)
	else
		p = plot(x, z[:, 1], fillrange = z[:, 2 * jump] .- 2 * shift, fillalpha = fillalpha, ylabel = ylabel, xlabel = xlabel, xticks = (0:maximum(x)/11:maximum(x), string.(reverse(reverse(x)[1:Int(floor(length(x)/11)):end]))), yticks = (z[end, 1:jump * step:end] .- shift .* range(0, stop = length(y[1:jump:end]) - 1, step = step), string.(y[1:jump * step:end])), legend = false, size = size, grid = false, ymirror = ymirror, palette = palette, tickfontsize = 12, labelfontsize = 12)
	end

	plot!(x, z[:, 1], color = color, linewidth = linewidth)
	plot!(x, z[:, 2 * jump] .- shift, color = color, linewidth = linewidth)

	for i in 2:Int(floor(length(y)/jump) - 1)
	plot!(x, z[:, i * jump] .- shift * (i - 1), fillrange = z[:, (i + 1) * jump] .- shift * i, fillalpha = fillalpha, palette = palette)
	plot!(x, z[:, i * jump] .- shift * (i - 1), color = color, linewidth = linewidth)
	plot!(x, z[:, (i + 1) * jump] .- shift * i, color = color, linewidth = linewidth)
	end
	return p
end

# ╔═╡ b41bfff2-c899-491e-b382-e7e3b136edbd
ridgeline(log.(reverse(T)), α, reverse(log.(real.(σ_dc)), dims = 1), jump = 1, shift = 1, xlabel = "Ω / ω", ylabel = "T / ω", size = (550, 550), linewidth = 1, palette = :thermal, fillalpha = 0.5, step = 1)

# ╔═╡ 22a8f0c4-f15a-46eb-adbe-4ee2575233b0
zero_temp_α = ridgeline(round.(Ω, digits = 1), α, hcat([real.(σ_0[:, i]) ./ maximum(real.(σ_0[:, i])) for i in 1:12]...), jump = 1, shift = 0.2, xlabel = "Ω / ω", ylabel = "α", size = (550, 550), palette = cgrad(:tab20, 12, categorical = true), linewidth = 1.5, fillalpha = 0.35, step = 1)

# ╔═╡ a6e1cd52-3384-44a5-aa15-be0f13691eea
savefig(zero_temp_α, "zero_temp_α.pdf")

# ╔═╡ 0aa5786a-f7f4-43b1-becf-f35a848302db
plot(T, F, labels = hcat(["α = $i" for i in α]...), legend = :outerright, linewidth = 2, linestyle = :auto, xlabel = "T / ω", ylabel = "Free energy / ω", tickfontsize = 12, legendfontsize = 12, labelfontsize = 12, size = (700, 550), minorgrid = true)

# ╔═╡ c1e13990-899e-4c65-81ff-ca3d08a46851
plot(T, v .- w, labels = hcat(["α = $i" for i in α]...), legend = :outerright, linewidth = 2, linestyle = :auto, xlabel = "T / ω", ylabel = "Free energy / ω", tickfontsize = 12, legendfontsize = 12, labelfontsize = 12, size = (700, 550), minorgrid = true, xaxis = :log)

# ╔═╡ c49ef6c5-61b6-4a40-9362-2166c642638d
free_energy = Plots.contourf(α, reverse(T), reverse(F, dims = 1), xlabel = "α", ylabel = "T / ω", fill = cgrad(:thermal, rev = false, categorical = true, scale = :log), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, xticks = (1:12, hcat(["$i" for i in 1:12]...)), yticks = ([0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0], hcat([0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0]...)), yaxis = :log, right_margin = 2Plots.mm)

# ╔═╡ d81b64a3-bf3d-443b-a712-783ae6f7eef3
savefig(free_energy, "free_energy.pdf")

# ╔═╡ 4020f1da-b186-48b9-b652-a480e20707a7
effective_mass = Plots.contourf(α, reverse(T), reverse(log.(M), dims = 1), xlabel = "α", ylabel = "T / ω", fill = cgrad(:thermal, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, xticks = (1:12, hcat(["$i" for i in 1:12]...)), yticks = ([0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0], hcat([0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0]...)), yaxis = :log)

# ╔═╡ ba9cc94a-26f3-45c8-b793-cc13ad62a28a
savefig(effective_mass, "effective_mass.pdf")

# ╔═╡ 2e68e769-79ff-4509-a9ec-3917237341bb
v_parameter = Plots.contourf(α, reverse(T), reverse(log.(v), dims = 1), xlabel = "α", ylabel = "T / ω", fill = cgrad(:thermal, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, xticks = (1:12, hcat(["$i" for i in 1:12]...)), yticks = ([0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0], hcat([0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0]...)), yaxis = :log)

# ╔═╡ 89b85619-ee78-46bf-8741-76d84741cc9f
savefig(v_parameter, "v_parameter.pdf")

# ╔═╡ db726be5-6572-4781-a407-1826a56dd1da
w_parameter = Plots.contourf(α, reverse(T), reverse(log.(w), dims = 1), xlabel = "α", ylabel = "T / ω", fill = cgrad(:thermal, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, xticks = (1:12, hcat(["$i" for i in 1:12]...)), yticks = ([0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0], hcat([0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0]...)), yaxis = :log)

# ╔═╡ a9cea3b1-55d8-4ac8-b140-120761e5daae
savefig(w_parameter, "w_parameter.pdf")

# ╔═╡ e4a79f98-c29a-4b9e-81fe-e9f55da2c71f
σ_contour_real = [Plots.contourf(Ω, reverse(T), reverse(log.(abs.(real.(σ[i])))', dims = 1), yaxis = :log, xlabel = "Ω / ω", ylabel = "T / ω", fill = cgrad(:viridis, rev = false, categorical = true, scale = :log), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, xticks = (0:4:20, hcat(["$i" for i in 0:4:20]...)), yticks = ([0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0], hcat([0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0]...)), right_margin = 7Plots.mm) for i in α]

# ╔═╡ c39d711d-c97e-4c8b-9f8e-304ef2be54de
begin
	for i in α
		savefig(σ_contour_real[i], "conductivity_contour_plots/conductivity_contour_real_$i.pdf")
	end
end

# ╔═╡ f8f1085a-ce7f-4854-a273-f8f9db5f8a74
σ_contour_imag = [Plots.contourf(Ω, reverse(T), reverse(log.(abs.(imag.(σ[i])))', dims = 1), yaxis = :log, xlabel = "Ω / ω", ylabel = "T / ω", fill = cgrad(:viridis, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, xticks = (0:4:20, hcat(["$i" for i in 0:4:20]...)), yticks = ([0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0], hcat([0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0]...)), right_margin = 7Plots.mm) for i in α]

# ╔═╡ c1d2a272-e4d4-4291-b424-f66cfde52135
begin
	for i in α
		savefig(σ_contour_imag[i], "conductivity_contour_plots/conductivity_contour_imag_$i.pdf")
	end
end

# ╔═╡ 67b1453c-e81a-4e7d-a7a5-a46ed60d4217
σ_contour_abs = [Plots.contourf(Ω, reverse(T), reverse(log.(abs.(abs.(σ[i])))', dims = 1), yaxis = :log, xlabel = "Ω / ω", ylabel = "T / ω", fill = cgrad(:viridis, rev = false, categorical = true), linewidth = 1.3, color = :grey10, size = (550, 500), tickfontsize = 12, labelfontsize = 12, xticks = (0:4:20, hcat(["$i" for i in 0:4:20]...)), yticks = ([0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0], hcat([0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0]...)), right_margin = 7Plots.mm) for i in α]

# ╔═╡ 5e4f5ba2-21e7-427c-93d4-2ebf23d204d4
begin
	for i in α
		savefig(σ_contour_abs[i], "conductivity_contour_plots/conductivity_contour_abs_$i.pdf")
	end
end

# ╔═╡ 67248e4a-babc-4b69-8bf2-34d56d1c70f4
σ_plot_temp_real = [ridgeline(round.(Ω, digits = 1), round.(reverse(T), digits = 2), reverse(log.(abs.(real.(σ[i]))), dims = 2), jump = 2, shift = 0.8, xlabel = "Ω / ω", ylabel = "T / ω", size = (550, 550), linewidth = 1, palette = :thermal, fillalpha = 0.75, step = 8) for i in α]

# ╔═╡ 5559b8fd-c862-4fdd-9a8b-4c196d401b46
begin
	for i in α
		savefig(σ_plot_temp_real[i], "conductivity_contour_plots/conductivity_plot_temp_real_$i.pdf")
	end
end

# ╔═╡ 6c86267c-799c-4585-a5ad-a221e8990d16
σ_plot_freq_real = [ridgeline(round.(reverse(T), digits = 2), round.(Ω, digits = 1), reverse(log.(abs.(real.(σ[i]'))), dims = 1), jump = 4, shift = 0.8, ylabel = "Ω / ω", xlabel = "T / ω", size = (600, 600), linewidth = 1, palette = :twilight, fillalpha = 0.75, step = 8, ymirror = true) for i in α]

# ╔═╡ f872eee7-e2a8-4173-a149-3c68597d421a
begin
	for i in α
		savefig(σ_plot_freq_real[i], "conductivity_contour_plots/conductivity_plot_freq_real_$i.pdf")
	end
end

# ╔═╡ 4a31fd0c-3cc4-4680-91a0-cf52e6ca75e3
σ_plot_temp_imag = [ridgeline(round.(Ω, digits = 1), reverse(β[1:end]), reverse(log.(abs.(imag.(σ[i]))), dims = 2), jump = 2, shift = 0.8, xlabel = "Ω (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 1, palette = :thermal, fillalpha = 0.75, step = 9) for i in α]

# ╔═╡ 39bbd242-e6e9-419e-9bc4-a5f183c33e58
begin
	for i in α
		savefig(σ_plot_temp_imag[i], "conductivity_contour_plots/conductivity_plot_temp_imag_$i.pdf")
	end
end

# ╔═╡ c052d7ff-58ce-413c-8358-851180759fdb
σ_plot_freq_imag = [ridgeline(round.(reverse(T), digits = 2), round.(Ω, digits = 1), reverse(log.(abs.(imag.(σ[i]'))), dims = 1), jump = 4, shift = 0.8, ylabel = "Ω / ω", xlabel = "T / ω", size = (550, 550), linewidth = 1, palette = :twilight, fillalpha = 0.75, step = 10, ymirror = true) for i in α]

# ╔═╡ e5ebc4f2-be00-4d5b-bad1-a71c6b5cdd93
begin
	for i in α
		savefig(σ_plot_freq_imag[i], "conductivity_contour_plots/conductivity_plot_freq_imag_$i.pdf")
	end
end

# ╔═╡ f21050fa-9718-4c23-aff0-de5474d6c5f6
σ_plot_temp_abs = [ridgeline(round.(Ω, digits = 1), reverse(β[1:end]), reverse(log.(abs.(σ[i])), dims = 2), jump = 2, shift = 0.8, xlabel = "Ω (THz)", ylabel = "T (K)", size = (550, 550), linewidth = 1, palette = :thermal, fillalpha = 0.75, step = 8) for i in α]

# ╔═╡ 71d9e25c-eec6-49bb-b43f-1927fb0618c8
begin
	for i in α
		savefig(σ_plot_temp_abs[i], "conductivity_contour_plots/conductivity_plot_temp_abs_$i.pdf")
	end
end

# ╔═╡ e3924eb5-8cfb-4db4-9d79-9f5dc7194a95
σ_plot_freq_abs = [ridgeline(round.(reverse(T), digits = 2), round.(Ω, digits = 1), reverse(log.(abs.(σ[i]')), dims = 1), jump = 4, shift = 0.8, ylabel = "Ω / ω", xlabel = "T / ω", size = (550, 550), linewidth = 1, palette = :twilight, fillalpha = 0.75, step = 8, ymirror = true) for i in α]

# ╔═╡ 87dfdce3-c94c-4f2a-9055-b699d37c1848
begin
	for i in α
		savefig(σ_plot_freq_abs[i], "conductivity_contour_plots/conductivity_plot_freq_abs_$i.pdf")
	end
end

# ╔═╡ Cell order:
# ╠═ec4e1740-0b51-11ec-3fa2-d3c1d67f4572
# ╠═9e67b5cd-f084-4bd4-9eb0-b94311184197
# ╠═b2a2fdb9-949a-498c-ac19-dab106582558
# ╠═72553f5c-276d-4256-8354-04154017d346
# ╠═67b93cb9-6335-448e-a09e-ead44db9d5c3
# ╠═d20b171e-1bd3-44aa-b067-7bf1a3c0227e
# ╠═8c909d72-7f28-49ed-9f6f-6b7f77e4e984
# ╠═2980129c-4be6-4928-8bfe-b6b13eb0d50a
# ╠═29d4a2e3-692b-48b3-9ce0-fca83b55e8d3
# ╠═eab11453-2f25-493a-84e4-307231496653
# ╠═136a409d-9ea0-40ab-998a-f5bf954315ca
# ╠═a1bef2ba-0fb0-4a43-a09d-4e2e5ab6e069
# ╠═671bed0e-082a-4b3a-9467-7be0ee2327df
# ╠═a7431fd9-588a-41f1-b3d8-ef15837ccf75
# ╠═f7746913-7afa-4bce-b261-9532bd6f27cb
# ╠═780dd0f5-b922-4691-94b5-fa9881381c07
# ╠═d3f88ec2-6f32-4429-a34b-4db33fc55433
# ╠═39cff408-72c6-42ea-8172-d9d0fb9e94fc
# ╠═8e432e19-984e-4727-ab9c-25894e2585d8
# ╠═83972098-4232-437e-9d2d-7620c17672f0
# ╠═d8762608-09eb-4ea0-9a84-5f46673a619d
# ╠═9ab20c3a-a0f5-4938-9339-720630fbd48b
# ╠═ed78a2eb-3313-4cae-b0d3-2118f5dff33f
# ╠═b41bfff2-c899-491e-b382-e7e3b136edbd
# ╠═f524e982-6208-44bc-86df-ccfa851b0c92
# ╠═9f09e259-014e-48b0-bd00-f7e2c9c897b3
# ╠═ccda1394-45b4-4a9a-8c14-425cb0b81a87
# ╠═fadad1b3-6da6-4724-be97-469d81641d8f
# ╠═3918264a-f72b-4b4c-b67f-cd6fa4c6c846
# ╠═c1d56910-a645-40ed-a5d8-0f7d5be19ae7
# ╠═22a8f0c4-f15a-46eb-adbe-4ee2575233b0
# ╠═a6e1cd52-3384-44a5-aa15-be0f13691eea
# ╠═19065839-85a3-4768-830e-26f64bf4a278
# ╠═7c978142-5319-45f8-b854-efac81271fcc
# ╠═0aa5786a-f7f4-43b1-becf-f35a848302db
# ╠═c1e13990-899e-4c65-81ff-ca3d08a46851
# ╠═c49ef6c5-61b6-4a40-9362-2166c642638d
# ╠═d81b64a3-bf3d-443b-a712-783ae6f7eef3
# ╠═4020f1da-b186-48b9-b652-a480e20707a7
# ╠═ba9cc94a-26f3-45c8-b793-cc13ad62a28a
# ╠═2e68e769-79ff-4509-a9ec-3917237341bb
# ╠═89b85619-ee78-46bf-8741-76d84741cc9f
# ╠═db726be5-6572-4781-a407-1826a56dd1da
# ╠═a9cea3b1-55d8-4ac8-b140-120761e5daae
# ╠═e4a79f98-c29a-4b9e-81fe-e9f55da2c71f
# ╠═c39d711d-c97e-4c8b-9f8e-304ef2be54de
# ╠═f8f1085a-ce7f-4854-a273-f8f9db5f8a74
# ╠═c1d2a272-e4d4-4291-b424-f66cfde52135
# ╠═67b1453c-e81a-4e7d-a7a5-a46ed60d4217
# ╠═5e4f5ba2-21e7-427c-93d4-2ebf23d204d4
# ╠═67248e4a-babc-4b69-8bf2-34d56d1c70f4
# ╠═5559b8fd-c862-4fdd-9a8b-4c196d401b46
# ╠═6c86267c-799c-4585-a5ad-a221e8990d16
# ╠═f872eee7-e2a8-4173-a149-3c68597d421a
# ╠═4a31fd0c-3cc4-4680-91a0-cf52e6ca75e3
# ╠═39bbd242-e6e9-419e-9bc4-a5f183c33e58
# ╠═c052d7ff-58ce-413c-8358-851180759fdb
# ╠═e5ebc4f2-be00-4d5b-bad1-a71c6b5cdd93
# ╠═f21050fa-9718-4c23-aff0-de5474d6c5f6
# ╠═71d9e25c-eec6-49bb-b43f-1927fb0618c8
# ╠═e3924eb5-8cfb-4db4-9d79-9f5dc7194a95
# ╠═87dfdce3-c94c-4f2a-9055-b699d37c1848

ext: jl


text:
struct MockSolver end

mutable struct CurentState
    route_idx::Int
    v_set_idx::Int
    v_idx::Int # inside a set idx
    act_idx::Int # inside a route
    t::Float64 # inside a route
    c::Dict{String,Float64} # inside a route
    CurentState() = new(1, 1, 1, 1, 0.0, Dict{String,Float64}())
end

function create_action(cur_state::CurentState, data::RvrpInstance,
                       loc::Location, act_type::Int; req_id = "")
    return Action(string("act_", cur_state.act_idx),
                       loc.id, act_type, req_id, cur_state.t)
end

function init_route(cur_state::CurentState, data::RvrpInstance,
                    computed_data::RvrpComputedData)

    v_cat = get_vehicle_category(data.vehicle_sets[cur_state.v_set_idx],
                                      data, computed_data)
    cur_state.c = Dict{String,Float64}(k => 0.0 for k in keys(v_cat.capacities.of_vehicle))
    cur_state.t = data.vehicle_sets[cur_state.v_set_idx].work_periods[1].soft_range.lb
    loc = get_departure_locations(data.vehicle_sets[cur_state.v_set_idx], data, computed_data)[1]

    act = create_action(cur_state, data, loc, 0)
    return Route(string("r_", cur_state.route_idx), data.vehicle_sets[cur_state.v_set_idx].id, [act], 0)
end

function finalize_route(r::Route, data::RvrpInstance,
                        computed_data::RvrpComputedData,
                        cur_state::CurentState, travel_times::Array{Float64,2})
    depot_loc = get_arrival_locations(data.vehicle_sets[cur_state.v_set_idx],
                                           data, computed_data)[1]
    if !isempty(travel_times)
        prev_location_idx = computed_data.location_id_2_index[r.sequence[end].location_id]
        cur_state.t += travel_times[prev_location_idx,depot_loc.index]
    end
    cur_state.act_idx += 1
    depot_act = create_action(cur_state, data, depot_loc, 0)
    push!(r.sequence, depot_act)

    cur_state.v_idx += 1
    cur_state.route_idx += 1
    cur_state.act_idx = 1
    # Check if need to change v_set_idx every time before initialize the route
    if cur_state.v_idx > data.vehicle_sets[cur_state.v_set_idx].nb_of_vehicles_range.soft_range.ub
        cur_state.v_set_idx += 1
        if cur_state.v_set_idx > length(data.vehicle_sets)
            return false
        end
        cur_state.v_idx = 1
    end
    return true
end

function check_if_feasible(req::Request, cur_state::CurentState,
                           r::Route, data::RvrpInstance,
                           computed_data::RvrpComputedData,
                           travel_times::Array{Float64,2})

    cap_cons = get_capacity_consumptions(req, data.product_specification_classes, computed_data)

    v_cat = get_vehicle_category(data.vehicle_sets[cur_state.v_set_idx],
                                      data, computed_data)
    for (k,v) in v_cat.capacities.of_vehicle
        if cur_state.c[k] + cap_cons[k] > v_cat.capacities.of_vehicle[k]
            return false, -1, cap_cons
        end
    end

    prev_loc_idx = computed_data.location_id_2_index[r.sequence[end].location_id]
    t = cur_state.t
    if req.request_type == 0 || req.request_type == 1
        p_loc_idx = get_pickup_locations(req, data, computed_data)[1].index
        # Time of arrival
        if !isempty(travel_times)
            t += travel_times[prev_loc_idx,p_loc_idx]
        end
        t = max(t, req.pickup_time_windows[1].soft_range.lb)
        if (t > req.pickup_time_windows[1].soft_range.ub ||
            t > data.locations[p_loc_idx].opening_time_windows[1].ub)
            return false, -1, cap_cons
        end
        t += req.pickup_service_time # Leaving time
        prev_loc_idx = p_loc_idx
    end

    if req.request_type == 0 || req.request_type == 2
        d_loc_idx = get_delivery_locations(req, data, computed_data)[1].index
        if !isempty(travel_times)
            t += travel_times[prev_loc_idx,d_loc_idx]
        end
        t = max(t, req.delivery_time_windows[1].soft_range.lb)
        if (t > req.delivery_time_windows[1].soft_range.ub ||
            t > data.locations[d_loc_idx].opening_time_windows[1].ub)
            return false, -1, cap_cons
        end
        t += req.delivery_service_time
        prev_loc_idx = d_loc_idx
    end

    depot_loc_idx = get_arrival_locations(data.vehicle_sets[cur_state.v_set_idx],
                                               data, computed_data)[1].index
    return_to_base_time = 0
    if !isempty(travel_times)
        return_to_base_time = travel_times[prev_loc_idx,depot_loc_idx]
    end
    if t + return_to_base_time > data.vehicle_sets[cur_state.v_set_idx].work_periods[1].soft_range.ub
        return false, -1, cap_cons
    end

    return true, t, cap_cons

end

function build_routes(data::RvrpInstance, solver::MockSolver,
                      computed_data::RvrpComputedData)
    unassigned = String[]
    travel_times = data.travel_specifications[1].travel_time_matrix
    routes = Route[]
    cur_state = CurentState()
    route = init_route(cur_state, data, computed_data)
    req_idx = 1
    while req_idx <= length(data.requests)
        req = data.requests[req_idx]
        feas, t, cap_cons = check_if_feasible(req, cur_state, route, data, computed_data, travel_times)
        if feas
            cur_state.t = t
            if req.request_type != 0
                for k in keys(cap_cons)
                    cur_state.c[k] += cap_cons[k]
                end
            end
            if req.request_type == 0 || req.request_type == 1
                p_loc = get_pickup_locations(req, data, computed_data)[1]
                cur_state.act_idx += 1
                p_act = create_action(cur_state, data, p_loc, 1, req_id = req.id)
                push!(route.sequence, p_act)
            end

            if req.request_type == 0 || req.request_type == 2
                d_loc = get_delivery_locations(req, data, computed_data)[1]
                cur_state.act_idx += 1
                d_act = create_action(cur_state, data, d_loc, 2, req_id = req.id)
                push!(route.sequence, d_act)
            end
            req_idx += 1
        end
        # If not feasible or the last req was just inserted
        if !feas || req_idx == length(data.requests) + 1
            feas = finalize_route(route, data, computed_data, cur_state, travel_times)
            push!(routes, route)
            if !feas
                unassigned = [req.id for req in data.requests[req_idx:end]]
                break
            end
            if req_idx <= length(data.requests) # If not all requests were done
                route = init_route(cur_state, data, computed_data)
            else
                break
            end
        end
    end
    return routes, unassigned
end

function solve(data::RvrpInstance, solver::MockSolver)

    id = string(data.id, "_mock_SOL_", rand(1:10000))
    if !check_initial_feasibility(data)
        return build_empty_sol(data, id)
    end
    computed_data = build_computed_data(data)
    routes, unassigned =  build_routes(data, solver, computed_data)
    problem_id = data.id
    sol = RvrpSolution(id, problem_id, 0.0, routes, unassigned)
    compute_cost_and_scheduled_times(sol, data, computed_data)
    return sol

end

mutable struct SolutionTimes
    current::Float64
    in_service::Float64
    in_travel::Float64
    in_waiting::Float64
    total_distance::Float64
end

function update_times_waiting(solution_times::SolutionTimes,
                              req::Request, act::Action,
                              waiting_time::Float64)
    # Waiting time based on the request
    if act.operation_type == 1
        tws = [_tw.soft_range for _tw in
               req.pickup_time_windows]
    elseif act.operation_type == 2
        tws = [_tw.soft_range for _tw in
               req.delivery_time_windows]
    end
    tw = get_closest_tw(tws, solution_times.current)
    waiting_time = max(0, max(waiting_time, tw.lb - solution_times.current))
    solution_times.in_waiting += waiting_time
    solution_times.current += waiting_time
    # Scheduled start time
    act.scheduled_start_time = solution_times.current
end

function update_times_service(solution_times::SolutionTimes,
                              req::Request, operation_type::Int)
    if operation_type == 1
        service_time = req.pickup_service_time
    elseif operation_type == 2
        service_time = req.delivery_service_time
    end
    solution_times.current += service_time
    solution_times.in_service += service_time
end

function compute_route_cost(solution_times::SolutionTimes,
                           v_set::HomogeneousVehicleSet)
    return (
        + v_set.cost_periods[1].travel_distance_unit_cost * solution_times.total_distance
        + v_set.cost_periods[1].travel_time_unit_cost * solution_times.in_travel
        + v_set.cost_periods[1].service_time_unit_cost * solution_times.in_service
        + v_set.cost_periods[1].waiting_time_unit_cost * solution_times.in_waiting
        + v_set.cost_periods[1].fixed_cost
    )
end

function compute_cost_and_scheduled_times(sol::RvrpSolution,
                                          data::RvrpInstance,
                                          computed_data::RvrpComputedData)

    cost = 0.0
    distances = data.travel_specifications[1].travel_distance_matrix
    times = data.travel_specifications[1].travel_time_matrix
    has_distances = !isempty(distances)
    has_times = !isempty(times)

    for r in sol.routes
        v_set = data.vehicle_sets[computed_data.vehicle_set_id_2_index[r.vehicle_set_id]]
        solution_times = SolutionTimes(v_set.work_periods[1].soft_range.lb,
                                       0.0, 0.0, 0.0, 0.0)
        prev_act = r.sequence[1]

        for act in r.sequence[2:end]

            prev_loc_idx = computed_data.location_id_2_index[prev_act.location_id]
            curr_loc_idx = computed_data.location_id_2_index[act.location_id]

            # Distance
            if has_distances
                solution_times.total_distance += distances[prev_loc_idx,curr_loc_idx]
            end

            # Travel time
            if has_times
                travel_time = times[prev_loc_idx,curr_loc_idx]
                solution_times.current += travel_time
                solution_times.in_travel += travel_time
            end

            # Waiting time based on the location
            tw = get_closest_tw(data.locations[prev_loc_idx].opening_time_windows,
                                solution_times.current)
            waiting_time = max(0, tw.lb - solution_times.current)

            req_id = act.request_id
            # If action is from a request
            if haskey(computed_data.request_id_2_index, req_id)
                req = data.requests[computed_data.request_id_2_index[req_id]]
                update_times_waiting(solution_times, req, act, waiting_time)
                update_times_service(solution_times, req, act.operation_type)
            end

            prev_act = act
        end
        # Update costs
        cost += compute_route_cost(solution_times, v_set)
    end
    sol.cost = cost
end

function supported_features(::Type{MockSolver})
    features = BitSet()

    # Location based features  #
    union!(features, HAS_OPENING_TIME_WINDOWS)

    # Product based features #
    union!(features, HAS_PRODUCT_CAPACITY_CONSUMPTIONS)
    union!(features, HAS_MULTIPLE_PRODUCT_CAPACITY_CONSUMPTIONS)

    # Request based features
    union!(features, HAS_SHIPMENT_REQUESTS)
    union!(features, HAS_PICKUPONLY_REQUESTS)
    union!(features, HAS_DELIVERYONLY_REQUESTS)
    union!(features, HAS_PICKUP_TIME_WINDOWS)
    union!(features, HAS_DELIVERY_TIME_WINDOWS)

    # VehicleCategory based features
    union!(features, HAS_VEHICLE_CAPACITIES)
    union!(features, HAS_MULTIPLE_VEHICLE_CAPACITIES)

    # HomogeneousVehicleSet based features
    union!(features, HAS_TRAVEL_TIME_UNIT_COST)
    union!(features, HAS_SERVICE_TIME_UNIT_COST)
    union!(features, HAS_WAITING_TIME_UNIT_COST)
    union!(features, HAS_TRAVEL_DISTANCE_UNIT_COST)
    union!(features, HAS_MAX_NB_VEHICLES)
    union!(features, HAS_WORKING_TIME_WINDOW)
    union!(features, HAS_FIXED_COST_PER_VEHICLE)
    union!(features, HAS_ARRIVAL_DIFFERENT_FROM_DEPARTURE)

    # Instance based features
    union!(features, HAS_MULTIPLE_VEHICLE_CATEGORIES)
    union!(features, HAS_MULTIPLE_VEHICLE_SETS)

    return [features]
end


text:
using CmeSIM.FluorophoreMarkers
using CmeSIM.Acquisitions
using Test
using Unitful: s, nm
using Images

FM = FluorophoreMarkers

@testset "Acquisition" begin

    @testset "fluorophore_marker.jl" begin

        @test begin
            marker = FluorophoreMarker(FM.dTomato)
            marker == FluorophoreMarker(FM.dTomato, 554nm, 581nm)
        end

        @test begin
            all(
                [typeof(FluorophoreMarker(ins)) for ins in instances(Fluorophore)] .== FluorophoreMarker,
            )
        end

        global marker = FluorophoreMarker(FM.EGFP)

    end

    @testset "acquisition.jl" begin
        images = unreconstructedimages
        images_rgb = unreconstructedimages_rgb

        @test begin
            acquisitiontimes = range(0s, 4s; length=18)
            typeof(Acquisition(images, acquisitiontimes, NA, M, marker, pixelsize)) ==
            Acquisition
        end

        @test begin
            typeof(Acquisition(images_rgb, 0.1s, NA, M, marker, pixelsize)) == Acquisition
        end

        global acquisition = Acquisition(images, 0.1s, NA, M, marker, pixelsize)

        @test begin
            all(times(acquisition) .== range(0s; length=18, step=0.1s))
        end

    end

    @testset "simacquisition.jl" begin
        @test begin
            acquisitiontimes = range(0s, 4s; length=18)
            typeof(SIMAcquisition(
                reconstructedimages_rgb,
                unreconstructedimages_rgb,
                acquisitiontimes,
                NA,
                M,
                marker,
                pixelsize,
            )) <: Acquisitions.AbstractAcquisition
        end

        @test begin
            typeof(SIMAcquisition(
                reconstructedimages, unreconstructedimages, 0.1s, NA, M, marker, pixelsize
            )) == SIMAcquisition
        end

        global simacquisition = SIMAcquisition(
            reconstructedimages, unreconstructedimages, 0.1s, NA, M, marker, pixelsize
        )

        @test begin
            all(times(simacquisition) .== range(0s; length=18, step=0.1s))
        end

    end

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "CMake"
version = v"3.22.2"

# Collection of sources required to build CMake
sources = [
    ArchiveSource("https://github.com/Kitware/CMake/releases/download/v$(version)/cmake-$(version).tar.gz",
                  "3c1c478b9650b107d452c5bd545c72e2fad4e37c09b89a1984b9a2f46df6aced"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/cmake-*/

cmake \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_BUILD_TYPE:STRING=Release \
    -DCMAKE_TOOLCHAIN_FILE=$CMAKE_TARGET_TOOLCHAIN

make -j${nproc}
make install
"""

# Build for all supported platforms.
platforms = expand_cxxstring_abis(supported_platforms())

# The products that we will ensure are always built
products = [
    ExecutableProduct("cmake", :cmake),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("OpenSSL_jll")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
#conjugate gradient solver
#Note on line searches- HZ much better than backtracking
function solveA!(loss, Θ::AbstractVolumePartsIter{TM,TV,T},
    Xv::AbstractXYIter{TM, TV, T},
    SolveType::Val{:bootstrap};
    RegressionType::Val = Val(PARAM[:iterbootstrapregressiontype])) where {TM, TV, T}

  dims = (K=size(Θ.A,1), T=size(Θ.A,2))

  if (Θ.G[1,1] .≈ Θ.G) |> all
    throw("All G values are the same!!! Implies G was not sourced
    bootstrapped correctly from the source")
  end

  #compute the loss and the values from A
  λ = updateA₁!(loss, Θ, Xv, RegressionType)
  updateAfromGAₜ!(Θ,Xv, 1)
  @assert λ ≈ lossfromA(Θ, Xv, RegressionType)

  #lower::Vector{T} = fill(lowerbound, (dims.T-1) * dims.K)
  #upper::Vector{T} = fill(upperbound, (dims.T-1) * dims.K)

  b = IOBuffer()
  write(b, "\n\n***Boostrap additional simulation parameters:\n" *
    "SolveType: $SolveType" *
    "measuretype: $(PARAM[:measuretype])\n" *
    "regressionmethod: $(PARAM[:regressionmethod])\n" *
    "controlsbykt: $(PARAM[:controlsbykt])\n" *
    "controlsbyt: $(PARAM[:controlsbyt])\n" *
    "controlsglobal: $(PARAM[:controlsglobal])\n" *
    "growth type: N/A (level since bootstrap)\n" *
    "momentum weight spec: $(PARAM[:momentumweightspec])\n" *
    "sorted default thresholds: $(PARAM[:sorteddefaultthresholds])\n" *
    "regression method (bootstrap): $(RegressionType)\n" *
    "Types (Matrix, Vector, T): $TM, $TV, $T\n" *
    "bounds: disabled\n"
    )
  write(b, "\n$(formatparameters())\n")

  iterinfo = String(take!(b))
  #println(iterinfo)
  open("$(PARAM[:outputpath])\\$(iterresultsname(λ))_info(bs).txt", "w+") do f
    write(f, iterinfo)
  end
  return λ
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module CommonSubexpressions

export @cse, cse

struct Cache
    args_to_symbol::Dict{Symbol, Symbol}
    disqualified_symbols::Set{Symbol}
    setup::Vector{Expr}
end

Cache() = Cache(Dict{Symbol,Symbol}(), Set{Symbol}(), Vector{Expr}())

function add_element!(cache::Cache, name, expr::Expr)
    sym = gensym(expr.args[1])
    cache.args_to_symbol[name] = sym
    push!(cache.setup, :($sym = $(expr)))
    sym
end

disqualify!(cache::Cache, x) = nothing
disqualify!(cache::Cache, s::Symbol) = push!(cache.disqualified_symbols, s)
disqualify!(cache::Cache, expr::Expr) = foreach(arg -> disqualify!(cache, arg), expr.args)

# fallback for non-Expr arguments
combine_subexprs!(setup, expr, warn_enabled::Bool) = expr

const standard_expression_forms = Set{Symbol}(
    (:call,
     :block,
     :comprehension,
     :(=>),
     :(:),
     :(&),
     :(&&),
     :(|),
     :(||),
     :tuple,
     :for,
     :ref,
     :macrocall,
     Symbol("'")))

const assignment_expression_forms = Set{Symbol}(
    (:(=),
     :(+=),
     :(-=),
     :(*=),
     :(/=)))

function combine_subexprs!(cache::Cache, expr::Expr, warn_enabled::Bool)
    if expr.head == :function
        # We can't continue CSE through a function definition, but we can
        # start over inside the body of the function:
        for i in 2:length(expr.args)
            expr.args[i] = combine_subexprs!(expr.args[i], warn_enabled)
        end
    elseif expr.head == :line
        # nothing
    elseif expr.head in assignment_expression_forms
        disqualify!(cache, expr.args[1])
        for i in 2:length(expr.args)
            expr.args[i] = combine_subexprs!(cache, expr.args[i], warn_enabled)
        end
    elseif expr.head == :generator
        for i in vcat(2:length(expr.args), 1)
            expr.args[i] = combine_subexprs!(cache, expr.args[i], warn_enabled)
        end
    elseif expr.head in standard_expression_forms
        for (i, child) in enumerate(expr.args)
            expr.args[i] = combine_subexprs!(cache, child, warn_enabled)
        end
        if expr.head == :call
            for (i, child) in enumerate(expr.args)
                expr.args[i] = combine_subexprs!(cache, child, warn_enabled)
            end
            if all(!isa(arg, Expr) && !(arg in cache.disqualified_symbols) for arg in expr.args)
                combined_args = Symbol(expr.args...)
                if !haskey(cache.args_to_symbol, combined_args)
                    sym = add_element!(cache, combined_args, expr)
                else
                    sym = cache.args_to_symbol[combined_args]
                end
                return sym
            else
            end
        end
    else
        warn_enabled && warn("CommonSubexpressions can't yet handle expressions of this form: $(expr.head)")
    end
    return expr
end

combine_subexprs!(x, warn_enabled::Bool = true) = x

function combine_subexprs!(expr::Expr, warn_enabled::Bool)
    cache = Cache()
    expr = combine_subexprs!(cache, expr, warn_enabled)
    Expr(:block, cache.setup..., expr)
end

macro cse(expr, warn_enabled::Bool = true)
    result = combine_subexprs!(expr, warn_enabled)
    # println(result)
    esc(result)
end

cse(expr, warn_enabled::Bool = true) = combine_subexprs!(copy(expr), warn_enabled)

end


text:
### Vector Policy ###
# maintained by @zsunberg and @etotheipluspi

"""
    VectorPolicy{S,A}
A generic MDP policy that consists of a vector of actions. The entry at `stateindex(mdp, s)` is the action that will be taken in state `s`.

# Fields
- `mdp::MDP{S,A}` the MDP problem
- `act::Vector{A}` a vector of size |S| mapping state indices to actions
"""
mutable struct VectorPolicy{S,A} <: Policy
    mdp::MDP{S,A}
    act::Vector{A}
end

action(p::VectorPolicy, s) = p.act[stateindex(p.mdp, s)]

"""
    VectorSolver{A}
Solver for VectorPolicy. Doesn't do any computation - just sets the action vector.

# Fields
- `act::Vector{A}` the action vector
"""
mutable struct VectorSolver{A}
    act::Vector{A}
end

function solve(s::VectorSolver{A}, mdp::MDP{S,A}) where {S,A}
    return VectorPolicy{S,A}(mdp, s.act)
end

function Base.show(io::IO, mime::MIME"text/plain", p::VectorPolicy)
    summary(io, p)
    println(io, ':')
    ds = get(io, :displaysize, displaysize(io))
    ioc = IOContext(io, :displaysize=>(first(ds)-1, last(ds)))
    showpolicy(ioc, mime, p.mdp, p)
end

"""
     ValuePolicy{P<:Union{POMDP,MDP}, T<:AbstractMatrix{Float64}, A}
A generic MDP policy that consists of a value table. The entry at `stateindex(mdp, s)` is the action that will be taken in state `s`.
It is expected that the order of the actions in the value table is consistent with the order of the actions in `act`.
If `act` is not explicitly set in the construction, `act` is ordered according to `actionindex`.

# Fields
- `mdp::P` the MDP problem
- `value_table::T` the value table as a |S|x|A| matrix
- `act::Vector{A}` the possible actions
"""
struct ValuePolicy{P<:Union{POMDP,MDP}, T<:AbstractMatrix{Float64}, A} <: Policy
    mdp::P
    value_table::T
    act::Vector{A}
end
function ValuePolicy(mdp::Union{MDP,POMDP}, value_table=zeros(length(states(mdp)), length(actions(mdp))))
    return ValuePolicy(mdp, value_table, ordered_actions(mdp))
end

action(p::ValuePolicy, s) = p.act[argmax(p.value_table[stateindex(p.mdp, s),:])]

actionvalues(p::ValuePolicy, s) = p.value_table[stateindex(p.mdp, s), :]

function Base.show(io::IO, mime::MIME"text/plain", p::ValuePolicy{M}) where M <: MDP
    summary(io, p)
    println(io, ':')
    ds = get(io, :displaysize, displaysize(io))
    ioc = IOContext(io, :displaysize=>(first(ds)-1, last(ds)))
    showpolicy(io, mime, p.mdp, p)
end

max_stars_repo_name: JuliaPOMDP/POMDPPolicies

max_issues_repo_name: JuliaPOMDP/POMDPPolicies

max_forks_repo_name: JuliaPOMDP/POMDPPolicies


text:
# ---
# title: 1373. Maximum Sum BST in Binary Tree
# id: problem1373
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming, Binary Search Tree
# link: <https://leetcode.com/problems/maximum-sum-bst-in-binary-tree/description/>
# hidden: true
# ---
#
# Given a **binary tree** `root`, the task is to return the maximum sum of all
# keys of **any**  sub-tree which is also a Binary Search Tree (BST).
#
# Assume a BST is defined as follows:
#
#   * The left subtree of a node contains only nodes with keys  **less than**  the node's key.
#   * The right subtree of a node contains only nodes with keys  **greater than**  the node's key.
#   * Both the left and right subtrees must also be binary search trees.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/01/30/sample_1_1709.png)
#
#
#
#     Input: root = [1,4,3,2,4,2,5,null,null,null,null,null,null,4,6]
#     Output: 20
#     Explanation: Maximum sum in a valid Binary search tree is obtained in root node with key equal to 3.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/01/30/sample_2_1709.png)
#
#
#
#     Input: root = [4,3,null,1,2]
#     Output: 2
#     Explanation: Maximum sum in a valid Binary search tree is obtained in a single root node with key equal to 2.
#
#
# **Example 3:**
#
#
#
#     Input: root = [-4,-2,-5]
#     Output: 0
#     Explanation: All values are negatives. Return an empty BST.
#
#
# **Example 4:**
#
#
#
#     Input: root = [2,1,3]
#     Output: 6
#
#
# **Example 5:**
#
#
#
#     Input: root = [5,4,8,3,null,6,3]
#     Output: 7
#
#
#
#
# **Constraints:**
#
#   * The given binary tree will have between `1` and `40000` nodes.
#   * Each node's value is between `[-4 * 10^4 , 4 * 10^4]`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# ================================================================
# Taylor model reach set
# ================================================================

using TaylorModels: TaylorModel1, TaylorN, fp_rpa

"""
    TaylorModelReachSet{N} <: AbstractTaylorModelReachSet{N}

Reach-set representation consisting of a vector of taylor models in one variable
(the "time" variable) whose coefficients are multivariate polynomials
(the "space" variables).

### Notes

The parameter `N` refers to the numerical type of the representation.
The space variables are assumed to be normalized to the interval `[-1, 1]`.
It is assumed that the time domain is the same for all components.
"""
struct TaylorModelReachSet{N} <: AbstractTaylorModelReachSet{N}
    X::Vector{TaylorModel1{TaylorN{N}, N}}
    Δt::TimeInterval
end

# interface functions
@inline set(R::TaylorModelReachSet) = R.X
@inline setrep(::TaylorModelReachSet{N}) where {N} = Vector{TaylorModel1{TaylorN{N}, N}}
@inline setrep(::Type{TaylorModelReachSet{N}}) where {N} = Vector{TaylorModel1{TaylorN{N}, N}}
@inline tstart(R::TaylorModelReachSet) = tstart(R.Δt) # t0 - δt + Δt0.lo
@inline tend(R::TaylorModelReachSet) = tend(R.Δt) # t0 + Δt0.hi
@inline tspan(R::TaylorModelReachSet) = R.Δt # # Interval(tstart(R), tend(R))
@inline dim(R::TaylorModelReachSet) = get_numvars()
@inline vars(R::TaylorModelReachSet) = Tuple(Base.OneTo(length(R.X)),)

# overload getter functions for the taylor model
# we assume that the first element is representative
domain(R::TaylorModelReachSet) = domain(first(R.X)) # normalized time domain
remainder(R::TaylorModelReachSet) = remainder.(R.X)
polynomial(R::TaylorModelReachSet) = polynomial.(R.X)
get_order(R::TaylorModelReachSet) = get_order.(R.X)
expansion_point(R::TaylorModelReachSet) = [Xi.x0 for Xi in R.X]

function shift(R::TaylorModelReachSet, t0::Number)
    return TaylorModelReachSet(set(R), tspan(R) + t0)
end

function reconstruct(R::TaylorModelReachSet, X)
    return TaylorModelReachSet(X, tspan(R))
end

# TODO remove but add comment in docs
function project(R::TaylorModelReachSet, vars::NTuple{D, M}) where {D, M<:Integer}
    throw(ArgumentError("the concrete projection of Taylor model reach-set is not " *
            "available; try first to overapproximate the Taylor model and the  project"))
end

# constructor with a time point
function TaylorModelReachSet(X::Vector{TaylorModel1{TaylorN{N}, N}}, t::Real) where {N}
    return TaylorModelReachSet(X, interval(t))
end

# ======================
# Intersection methods
# ======================

# for TMs, we overapproximate with a zonotope
function _intersection(R::TaylorModelReachSet, Y::LazySet)
    _intersection(R, Y, FallbackIntersection())
end

function _intersection(R::TaylorModelReachSet, Y::LazySet, ::FallbackIntersection)
    X = set(overapproximate(R, Zonotope))
    intersection(X, Y)
end

# TODO overapprox Y with a box if it's bounded?
function _intersection(R::TaylorModelReachSet, Y::LazySet, ::BoxIntersection)
    X = set(overapproximate(R, Hyperrectangle))
    out = intersection(X, Y)
    return isempty(out) ? EmptySet(dim(out)) : overapproximate(out, Hyperrectangle)
end

# TODO refactor? See LazySets#2158
function _intersection(R::TaylorModelReachSet, Y::UnionSetArray{N, HT}, ::BoxIntersection) where {N, HT<:HalfSpace{N}}
    X = set(overapproximate(R, Hyperrectangle))

    # find first non-empty intersection
    m = length(Y.array) # can't use array(Y) ?
    i = 1
    local W
    @inbounds while i <= m
        W = intersection(X, Y.array[i])
        !isempty(W) && break
        i += 1
    end
    if i == m+1
        return EmptySet(dim(Y))
    end

    # add all other non-empty intersections
    out = [W]
    @inbounds for j in i+1:m
        W = intersection(X, Y.array[j])
        !isempty(W) && push!(out, W)
    end
    return isempty(out) ? EmptySet(dim(out)) : overapproximate(UnionSetArray(out), Hyperrectangle)
end

# ======================
# Inclusion checks
# ======================

function Base.:⊆(R::TaylorModelReachSet, X::LazySet)
    Rz = overapproximate(R, Zonotope)
    return set(Rz) ⊆ X
end

# ======================
# Disjointness checks
# ======================

_is_intersection_empty(R::TaylorModelReachSet, Y::LazySet, method::FallbackDisjointness) = is_intersection_empty(set(overapproximate(R, Zonotope)), Y)
_is_intersection_empty(R::TaylorModelReachSet, Y::LazySet, method::ZonotopeEnclosure) = is_intersection_empty(set(overapproximate(R, Zonotope)), Y)

# FIXME when UnionSet, UnionSetArray <: LazySet
_is_intersection_empty(R::TaylorModelReachSet, Y::Union{UnionSet, UnionSetArray}, method::FallbackDisjointness) = is_intersection_empty(set(overapproximate(R, Zonotope)), Y)
_is_intersection_empty(R::TaylorModelReachSet, Y::Union{UnionSet, UnionSetArray}, method::ZonotopeEnclosure) = is_intersection_empty(set(overapproximate(R, Zonotope)), Y)
_is_intersection_empty(R::TaylorModelReachSet, Y::Union{UnionSet, UnionSetArray}, method::BoxEnclosure) = is_intersection_empty(set(overapproximate(R, Hyperrectangle)), Y)

_is_intersection_empty(R1::TaylorModelReachSet, R2::TaylorModelReachSet, method::FallbackDisjointness) = is_intersection_empty(set(overapproximate(R1, Zonotope)), set(overapproximate(R2, Zonotope)))
_is_intersection_empty(R1::TaylorModelReachSet, R2::TaylorModelReachSet, method::ZonotopeEnclosure) = is_intersection_empty(set(overapproximate(R1, Zonotope)), set(overapproximate(R2, Zonotope)))

@commutative function _is_intersection_empty(R::TaylorModelReachSet, Y::Union{LazySet, UnionSet, UnionSetArray}, method::ZonotopeEnclosure)
    Z = overapproximate(R, Zonotope)
    return is_intersection_empty(set(Z), Y)
end

@commutative function _is_intersection_empty(R::TaylorModelReachSet, X::Universe, method::AbstractDisjointnessMethod)
    return false
end

@commutative function _is_intersection_empty(R::TaylorModelReachSet, X::Universe, method::ZonotopeEnclosure)
    return false
end

@commutative function _is_intersection_empty(R::TaylorModelReachSet, X::Universe, method::BoxEnclosure)
    return false
end

# ==============================================================
# Conversion and overapproximation of Taylor model reach-sets
# ==============================================================

# no-op
function overapproximate(R::TaylorModelReachSet{N}, ::Type{<:TaylorModelReachSet}) where {N}
    return R
end

function overapproximate(R::TaylorModelReachSet{N}, ::BoxEnclosure) where {N}
    return overapproximate(R, Hyperrectangle)
end

function overapproximate(R::TaylorModelReachSet{N}, ::ZonotopeEnclosure) where {N}
    return overapproximate(R, Zonotope)
end

# tdom defined in the methods below is the same as Δt - t0, but the domain inclusion check
# in TM.evauate may fail, so we unpack the domain again here; also note that
# by construction the TMs in time are centered at zero

# overapproximate taylor model reachset with one hyperrectangle
function _overapproximate(R::TaylorModelReachSet{N}, ::Type{<:Hyperrectangle}) where {N}
    # dimension of the reachset
    D = dim(R)

    # normalized time domain
    tdom = domain(R)

    # evaluate the Taylor model in time
    # X_Δt is a vector of TaylorN (spatial variables) whose coefficients are intervals
    X_Δt = evaluate(set(R), tdom)

    # evaluate the spatial variables in the symmetric box
    Bn = symBox(D)
    X̂ib = IntervalBox([evaluate(X_Δt[i], Bn) for i in 1:D]...)
    X̂ = convert(Hyperrectangle, X̂ib)

    Δt = tspan(R)
    return ReachSet(X̂, Δt)
end

LazySets.box_approximation(R::TaylorModelReachSet) = _overapproximate(R, Hyperrectangle)

# overapproximate taylor model reachset with several hyperrectangles,
# by splitting in time (ntdiv), and/or space (with either nsdiv for uniform partitions
# or by specifying a partition (vector of integers))
# e.g. if the partition is uniform for each dimension, the number of returned sets
# is nsdiv^D * ntdiv
function overapproximate(R::TaylorModelReachSet{N}, ::Type{<:Hyperrectangle};
                         partition=nothing, nsdiv=1, ntdiv=1) where {N}

    # no splitting
    if isnothing(partition) && nsdiv == 1 && ntdiv == 1
        return _overapproximate(R, Hyperrectangle)
    end

    D = dim(R)
    X = set(R)

    # time domain splittig
    tspdiv = IntervalArithmetic.mince(tspan(R), ntdiv)
    tdom = domain(R)
    Δtdiv = IntervalArithmetic.mince(tdom, ntdiv)

    # evaluate the Taylor model in (normalized) time
    # X_Δt is a vector of TaylorN (spatial variables) whose coefficients are intervals
    X_Δt = [evaluate(X, Δti) for Δti in Δtdiv]

    # partition the symmetric box
    S = symBox(D)
    if isnothing(partition)
        if nsdiv == 1
            partition = [S]

        else
            # FIXME (also below) may use IntervalArithmetic.mince directly
            partition = fill(nsdiv, D)
            Sdiv = LazySets.split(convert(Hyperrectangle, S), partition)
            partition = convert.(IntervalBox, Sdiv)
        end

    else
        Sdiv = LazySets.split(convert(Hyperrectangle, S), partition)
        partition = convert.(IntervalBox, Sdiv)
    end
    nparts = length(partition)

    # preallocate output vectors
    ST = Hyperrectangle{N, SVector{D, N}, SVector{D, N}}
    X̂out = Vector{ST}(undef, nparts * ntdiv)
    R̂out = Vector{ReachSet{N, ST}}(undef, nparts * ntdiv)

    # evaluate the spatial variables in the symmetric box
    @inbounds for k in 1:ntdiv
        for j in 1:nparts
            Bj = partition[j]
            Xk = X_Δt[k]
            X̂ib = IntervalBox([evaluate(Xk[i], Bj) for i in 1:D])

            idx = j + (k-1)*nparts
            X̂out[idx] = convert(Hyperrectangle, X̂ib)

            R̂out[idx] = ReachSet(X̂out[idx], tspdiv[k])
        end
    end

    return R̂out
end

# overapproximate taylor model reachset with one zonotope
function overapproximate(R::TaylorModelReachSet{N}, ::Type{<:Zonotope}) where {N}
    # dimension of the reachset
    n = dim(R)

    # normalized time domain
    tdom = domain(R)

    # evaluate the Taylor model in time
    # X_Δt is a vector of TaylorN (spatial variables) whose coefficients are intervals
    X = set(R)
    X_Δt = evaluate(X, tdom)

    # builds the associated taylor model for each coordinate j = 1...n
    #  X̂ is a TaylorModelN whose coefficients are intervals
    X̂ = [TaylorModelN(X_Δt[j], zeroI, zeroBox(n), symBox(n)) for j in 1:n]

    # compute floating point rigorous polynomial approximation
    # fX̂ is a TaylorModelN whose coefficients are floats
    fX̂ = fp_rpa.(X̂)

    Δt = tspan(R)
    # LazySets can overapproximate a Taylor model with a Zonotope
    Zi = overapproximate(fX̂, Zonotope)
    return ReachSet(Zi, Δt)
end

# overapproximate taylor model reachset with several zonotopes given the partition,
# that specifies the way in which the symmetric box [-1, 1]^n is split:
#
# - if `partition` is an integer, makes as uniform partition for each coordinate
# - if `partition` is a vector of integers, split the domain according to each element in partition
function overapproximate(R::TaylorModelReachSet{N}, ::Type{<:Zonotope}, partition::Union{Int, Vector{Int}}) where {N}
    # dimension of the reachsets
    D = dim(R)

    # normalized time domain
    tdom = domain(R)

    # evaluate the Taylor model in time
    # X_Δt is a vector of TaylorN (spatial variables) whose coefficients are intervals
    X = set(R)
    X_Δt = evaluate(X, tdom)

    # evaluate the spatial variables in the symmetric box
    part = _split_symmetric_box(D, partition)
    fX̂ = Vector{Vector{TaylorModelN{length(X_Δt), N, N}}}(undef, length(part))
    @inbounds for (i, Bi) in enumerate(part)
        x0 = IntervalBox(mid.(Bi))
        X̂ib = [TaylorModelN(X_Δt[j], zeroI, x0, Bi) for j in 1:D]
        fX̂[i] = fp_rpa.(X̂ib)
    end
    Z = overapproximate.(fX̂, Zonotope)
    Δt = tspan(R)
    #return ReachSet(UnionSetArray(Z), Δt) # but UnionSetArray is not yet a lazyset
    return ReachSet(ConvexHullArray(Z), Δt)
end

# evaluate at a given time and overapproximate the resulting set with a zonotope
function overapproximate(R::TaylorModelReachSet{N}, ::Type{<:Zonotope}, t::AbstractFloat;
                         remove_zero_generators=true) where {N}
    @assert t ∈ tspan(R) "the given time point $t does not belong to the reach-set's time span, $(tspan(R))"

    X = set(R)
    Δtn = domain(R)
    if t == tstart(R)
        tn = zero(N)
    elseif t == tend(R)
        tn = Δtn.hi
    else
        tn = t - tstart(R)
    end
    X_Δt = evaluate(X, tn)
    n = dim(R)
    X̂ = [TaylorModelN(X_Δt[j], zeroI, zeroBox(n), symBox(n)) for j in 1:n]
    fX̂ = fp_rpa.(X̂)
    Zi = overapproximate(fX̂, Zonotope, remove_zero_generators=remove_zero_generators)
    Δt = TimeInterval(t, t)
    return ReachSet(Zi, Δt)
end

function overapproximate(R::TaylorModelReachSet{N}, ::Type{<:Zonotope}, Δt::TimeInterval) where {N}
    @assert Δt ⊆ tspan(R) "the given time point $t does not belong to the reach-set's time span, $(tspan(R))"

    X = set(R)
    Δtn = domain(R)

    if Δt.lo == tstart(R)
        dtn_lo = zero(N)
    elseif Δt.lo == tend(R)
        dtn_lo = Δtn.hi
    else
        dtn_lo = Δt.lo - tstart(R)
    end

    if Δt.hi == tstart(R)
        dtn_hi = zero(N)
    elseif Δt.hi == tend(R)
        dtn_hi = Δtn.hi
    else
        dtn_hi = Δt - tstart(R)
    end

    dtn = IA.Interval(dtn_lo, dtn_hi)
    X_Δt = evaluate(X, dtn)
    n = dim(R)
    X̂ = [TaylorModelN(X_Δt[j], zeroI, zeroBox(n), symBox(n)) for j in 1:n]
    fX̂ = fp_rpa.(X̂)
    Zi = overapproximate(fX̂, Zonotope)

    return ReachSet(Zi, Δt)
end

# convert a hyperrectangular set to a taylor model reachset
function convert(::Type{<:TaylorModelReachSet}, H::AbstractHyperrectangle{N};
                 orderQ::Integer=2, orderT::Integer=8, Δt::TimeInterval=zeroI) where {N}

    n = dim(H)
    x = set_variables("x", numvars=n, order=orderQ)

    # preallocations
    vTM = Vector{TaylorModel1{TaylorN{N}, N}}(undef, n)

    # normalized time domain
    Δtn = TimeInterval(zero(N), diam(Δt))

    # for each variable i = 1, .., n, compute the linear polynomial that covers
    # the line segment corresponding to the i-th edge of H
    @inbounds for i in 1:n
        α = radius_hyperrectangle(H, i)
        β = center(H, i)
        pi = α * x[i] + β
        vTM[i] = TaylorModel1(Taylor1(pi, orderT), zeroI, zeroI, Δtn)
    end

    return TaylorModelReachSet(vTM, Δt)
end

# overapproximate a hyperrectangular set with a taylor model reachset, fallback to convert
function overapproximate(H::AbstractHyperrectangle{N}, T::Type{<:TaylorModelReachSet};
                         orderQ::Integer=2, orderT::Integer=8, Δt::TimeInterval=zeroI) where {N}
    convert(T, H; orderQ=orderQ, orderT=orderT, Δt=Δt)
end

# overapproximate a zonotopic set with a taylor model reachset
# FIXME pass algorithm option to choose between using parallelotope oa or  order reduction
function overapproximate(Z::AbstractZonotope{N}, ::Type{<:TaylorModelReachSet};
                         orderQ::Integer=2, orderT::Integer=8, Δt::TimeInterval=zeroI,
                         indices=1:dim(Z), box_reduction=false) where {N}

    n = dim(Z)
    x = set_variables("x", numvars=n, order=orderQ)

    if order(Z) > 1
        if box_reduction
            # diagonal generators matrix
            Z = _reduce_order(Z, 1)
        else
            # indices selects the indices that we want to keep
            Z = _overapproximate_hparallelotope(Z, indices)
        end
    end
    c = center(Z)
    G = genmat(Z)

    # preallocations
    vTM = Vector{TaylorModel1{TaylorN{N}, N}}(undef, n)

    # normalized time domain
    Δtn = TimeInterval(zero(N), diam(Δt))

    # for each variable i = 1, .., n, compute the linear polynomial that covers
    # the line segment corresponding to the i-th edge of Z
    @inbounds for i in 1:n
        pi = c[i] + sum(view(G, i, :) .* x)
        vTM[i] = TaylorModel1(Taylor1(pi, orderT), zeroI, zeroI, Δtn)
    end

    return TaylorModelReachSet(vTM, Δt)
end

# convert a zonotopic set to a taylor model reachset (only if it has order 1)
function convert(TM::Type{<:TaylorModelReachSet}, Z::AbstractZonotope; kwargs...)
    if order(Z) == 1
        return overapproximate(Z, TM; kwargs...)

    else
        throw(ArgumentError("can't convert a zonotope of order $(order(Z)) to a `TaylorModelReachSet`; " *
                            "use `overapproximate(Z, TaylorModelReachSet)` instead"))
    end
end

function convert(TM::Type{<:TaylorModelReachSet}, R::AbstractLazyReachSet; kwargs...)
    return convert(TM, set(R); Δt=tspan(R), kwargs...)
end

function overapproximate(R::AbstractLazyReachSet, TM::Type{<:TaylorModelReachSet}; kwargs...)
    return overapproximate(set(R), TM; Δt=tspan(R), kwargs...)
end

# zonotope of order 2 with generators matrix G = [M; D] where M is n x n and D is n x n and diagonal
function _overapproximate_structured(Z::AbstractZonotope{N}, ::Type{<:TaylorModelReachSet};
                                     orderQ::Integer=2, orderT::Integer=8, Δt::TimeInterval=zeroI) where {N}
    n = dim(Z)
    x = set_variables("x", numvars=n, order=orderQ)

    # check structure
    order(Z) == 2 || throw(ArgumentError("this function requires that the order of the zonotope is 2, got $(order(Z))"))

    c = LazySets.center(Z)
    G = genmat(Z)

    M = view(G, :, 1:n)
    D = view(G, :, n+1:2n)
    isdiag(D) || throw(ArgumentError("the columns $(n+1) to $(2n) of the generators matrix do not form a diagonal matrix"))

    # preallocations
    vTM = Vector{TaylorModel1{TaylorN{N}, N}}(undef, n)

    # normalized time domain
    Δtn = TimeInterval(zero(N), diam(Δt))

    # for each variable i = 1, .., n, compute the linear polynomial that covers
    # the line segment corresponding to the i-th edge of Z
    @inbounds for i in 1:n
        pi = c[i] + sum(view(M, i, :) .* x)
        di = abs(D[i, i])
        rem = interval(-di, di)
        vTM[i] = TaylorModel1(Taylor1(pi, orderT), rem, zeroI, Δtn)
    end

    return TaylorModelReachSet(vTM, Δt)
end

function _overapproximate_structured(Zcp::CartesianProduct{N, <:Zonotope, <:Interval}, ::Type{<:TaylorModelReachSet};
                                     orderQ::Integer=2, orderT::Integer=8, Δt::TimeInterval=zeroI) where {N}
    n = dim(Zcp)
    x = set_variables("x", numvars=n, order=orderQ)

    # check structure
    Z = Zcp.X
    Y = Zcp.Y
    if order(Z) > 1
        Z = remove_redundant_generators(Z)
    end
    order(Z) == 1 || throw(ArgumentError("this function requires that the order of the (lower) zonotope is 1, got $(order(Z))"))

    c = LazySets.center(Z)
    G = genmat(Z)

    # preallocations
    vTM = Vector{TaylorModel1{TaylorN{N}, N}}(undef, n)

    # normalized time domain
    Δtn = TimeInterval(zero(N), diam(Δt))

    xstate = x[1:n-1]
    @inbounds for i in 1:n-1
        pi = c[i] + sum(view(G, i, :) .* xstate) + zero(TaylorN(n, order=orderQ))
        rem = interval(0)
        vTM[i] = TaylorModel1(Taylor1(pi, orderT), rem, zeroI, Δtn)
    end
    # diagonal part
    @inbounds begin
        pi = mid(Y.dat) + zero(TaylorN(1, order=orderQ))
        d = diam(Y.dat) / 2
        rem = interval(-d, d)
        vTM[n] = TaylorModel1(Taylor1(pi, orderT), rem, zeroI, Δtn)
    end
    return TaylorModelReachSet(vTM, Δt)
end

# we assume that if n = size(G, 1), then:
# (size(G) == (n, 2n + 1)) && isdiag(view(G, :, (n+2):(2n+1)))
function _overapproximate_structured_full(Zcp::CartesianProduct{N, <:Zonotope, <:Interval}, ::Type{<:TaylorModelReachSet};
                                          orderQ::Integer=2, orderT::Integer=8, Δt::TimeInterval=zeroI) where {N}
    n = dim(Zcp) - 1
    x = set_variables("x", numvars=n+1, order=orderQ)

    # check structure
    # not checking structure
    Z = Zcp.X
    c = LazySets.center(Z)
    G = genmat(Z)

    # preallocations
    vTM = Vector{TaylorModel1{TaylorN{N}, N}}(undef, n+1)

    # normalized time domain
    Δtn = TimeInterval(zero(N), diam(Δt))

    # fill rows corresponding to the "zonotope" variables: 1 to nth-variables
    @inbounds for i in 1:n
        pi = c[i] + sum(view(G, i, 1:(n+1)) .* x) + zero(TaylorN(n+1, order=orderQ))
        d = abs(G[i, n + 1 + i])
        rem = interval(-d, d)
        vTM[i] = TaylorModel1(Taylor1(pi, orderT), rem, zeroI, Δtn)
    end

    # fill the final row, which correspponds to the "interval" variable: (n+1)-th
    I = Zcp.Y.dat
    pi = mid(I) + zero(TaylorN(n+1, order=orderQ))
    d = diam(I) / 2
    rem = interval(-d, d)
    @inbounds vTM[n+1] = TaylorModel1(Taylor1(pi, orderT), rem, zeroI, Δtn)

    return TaylorModelReachSet(vTM, Δt)
end

# ======================
# Remainder handling
# ======================

using TaylorModels: shrink_wrapping!

function _shrink_wrapping(R::TaylorModelReachSet)
    rem = remainder(R)
    dt = domain(R)
    n = dim(R)

    # if all remainders are zero => nothing to do
    all(iszero, rem) && return R

    # transform into a TaylorN whose coeffs are floats
    X = set(R)
    Xev = evaluate.(X, dt)
    W = [TaylorModelN(Xev[i], zeroI, zeroBox(n), symBox(n)) for i in 1:n]
    Wfp = fp_rpa.(W)

    # absorb remainder in the polynomial part
    shrink_wrapping!(Wfp)

    # transform back to a TaylorModel1 of TaylorN
    orderT = get_order(set(R)[1])
    p = [Taylor1(TaylorN(polynomial(Wfp[i])), orderT) for i in 1:n]
    Y = [TaylorModel1(p[i], zeroI, zeroI, dt) for i in 1:n]

    return TaylorModelReachSet(Y, tspan(R))
end


text:
# ------------------------------------------------------------------------------------------
# # Iterators
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# `for` loops "lower" to `while` loops plus calls to the `iterate` function:
#
# ```jl
# for i in iter   # or  "for i = iter" or "for i ∈ iter"
#     # body
# end
# ```
#
# internally works the same as:
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ```jl
# next = iterate(iter)
# while next !== nothing
#     (i, state) = next
#     # body
#     next = iterate(iter, state)
# end
# ```
#
# The same applies to comprehensions and generators.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Note `nothing` is a singleton value (the only value of its type `Nothing`) used by
# convention when there is no value to return (a bit like `void` in C). For example
# ------------------------------------------------------------------------------------------

typeof(print("hello"))

A = ['a','b','c'];

iterate(A)

iterate(A, 2)

iterate(A, 3)

iterate(A, 4)

# ------------------------------------------------------------------------------------------
# Iteration is also used by "destructuring" assignment:
# ------------------------------------------------------------------------------------------

x, y = A

x

y

# ------------------------------------------------------------------------------------------
# Yet another user of this "iteration protocol" is so-called argument "splatting":
# ------------------------------------------------------------------------------------------

string(A)

string('a','b','c')

string(A...)

# ------------------------------------------------------------------------------------------
# ## Iteration utilities
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# `collect` gives you all elements of an iterator as an array.
# Comprehensions are actually equivalent to calling `collect` on a generator.
# ------------------------------------------------------------------------------------------

collect(pairs(A))

collect(zip(100:102,A))

# ------------------------------------------------------------------------------------------
# Some other favorites to experiment with. These are in the built-in `Iterators` module:
# - `enumerate`
# - `rest`
# - `take`
# - `drop`
# - `product`
# - `flatten`
# - `partition`
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Some iterators are infinite!
# - `countfrom`
# - `repeated`
# - `cycle`
# ------------------------------------------------------------------------------------------

I = zip(Iterators.cycle(0:1), Iterators.flatten([[2,3],[4,5]]))

collect(I)

collect(Iterators.product(I,A))

string(I...)

# ------------------------------------------------------------------------------------------
# ## Defining iterators
# ------------------------------------------------------------------------------------------

struct SimpleRange
    lo::Int
    hi::Int
end

Base.iterate(r::SimpleRange, state = r.lo) = state > r.hi ? nothing : (state, state+1)

Base.length(r::SimpleRange) = r.hi-r.lo+1

collect(SimpleRange(2,8))

# ------------------------------------------------------------------------------------------
# ## Iterator traits
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# For many algorithms, it's useful to know certain properties of an iterator up front.
#
# The most useful is whether an iterator has a fixed, known length.
# ------------------------------------------------------------------------------------------

Base.IteratorSize([1])

Base.IteratorSize(Iterators.repeated(1))

Base.IteratorSize(eachline(open("/dev/null")))

# ------------------------------------------------------------------------------------------
# ## Exercise
#
# Define an iterator giving the first N fibonacci numbers.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Index iterators
# ------------------------------------------------------------------------------------------

A = rand(3,5)

eachindex(A)

keys(A)

Av = view(A, [1,3], [1,2,5])

A[[1,3],[1,2,5]]

eachindex(Av)

# ------------------------------------------------------------------------------------------
# ### Example: $3\times 3\times \dots \times3$ boxcar filter (from a blog post by Tim Holy)
# ------------------------------------------------------------------------------------------

function boxcar3(A::AbstractArray)
    out = similar(A)
    R = CartesianIndices(size(A))
    I1, Iend = first(R), last(R)
    for I in R
        n, s = 0, zero(eltype(out))
        for J in CartesianIndices(map(:, max(I1, I-I1).I, min(Iend, I+I1).I))
            s += A[J]
            n += 1
        end
        out[I] = s/n
    end
    out
end

using Images

A = rand(256,256);

Gray.(A)

Gray.(boxcar3(A))

function sumalongdims!(B, A)
    # It's assumed that B has size 1 along any dimension that we're summing
    fill!(B, 0)
    Bmax = CartesianIndex(size(B))
    for I in CartesianIndices(size(A))
        B[min(Bmax,I)] += A[I]
    end
    B
end

B = zeros(1, 256)

sumalongdims!(B, A)

reduce(+,A,dims=(1,))

# ------------------------------------------------------------------------------------------
# `CartesianIndices` and other "N-d" iterators have a shape that propagates through
# generators.
# ------------------------------------------------------------------------------------------

[1 for i in CartesianIndices((2,3))]

B = rand(5,5)

view(B,CartesianIndices((2,3)))

# ------------------------------------------------------------------------------------------
# ## Exercise: CartesianIndex life!
#
# - Write a function `neighborhood(A::Array, I::CartesianIndex)` that returns a view of the
# 3x3 neighborhood around a location
# - Write a function `liferule(A, I)` that implements the evolution rule of Conway's life
# cellular automaton:
#   - 2 live neighbors $\rightarrow$ stay the same
#   - 3 live neighbors $\rightarrow$ 1
#   - otherwise $\rightarrow$ 0
# - Write a function `life(A)` that maps A to the next life step using these
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Some famous initial conditions:
# ------------------------------------------------------------------------------------------

A = fill(0, 128,128);

A[61:63,61:63] = [1 1 0
                  0 1 1
                  0 1 0]

A = life(A)
# `repeat` can be used to get chunky pixels to make the output easier to see
Gray.(repeat(A,inner=(4,4)))

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
function +(ps::MPoly...)
    pu = PolyUnion(ps)
    r = zero(pu)
    for pindex = 1:length(ps)
        for (m, c) in ps[pindex]
            r[newexps(pu, m, pindex)] += c
        end
    end
    r
end

function +(p1::MPoly, p2::MPoly)
    pu = PolyUnion([p1, p2])
    r = zero(pu)
    for (m, c) in p1
        r[newexps(pu, m, 1)] += c
    end
    for (m, c) in p2
        r[newexps(pu, m, 2)] += c
    end
    r
end

function +{T,U}(p::MPoly{T}, s::U)
    r = copy(convert(MPoly{promote_type(T,U)}, p))
    r[zeros(Int, nvars(p))] += s
    r
end

+(s, p::MPoly) = p + s

function -(p::MPoly, q::MPoly)
    pu = PolyUnion([p, q])
    r = zero(pu)
    for (m, c) in p
        r[newexps(pu, m, 1)] += c
    end
    for (m, c) in q
        r[newexps(pu, m, 2)] -= c
    end
    r
end

function -{T,U}(p::MPoly{T}, s::U)
    r = copy(convert(MPoly{promote_type(T, U)}, p))
    r[zeros(Int, nvars(p))] -= s
    r
end

function -(p::MPoly)
    r = zero(p)
    for (m, c) in p
        r[m] = -c
    end
    r
end

-(s, p::MPoly) = -p + s

function *(p1::MPoly, p2::MPoly)
    pu = PolyUnion([p1, p2])
    r = zero(pu)
    for (m1, c1) in p1
        for (m2, c2) in p2
            r[newexps(pu, m1, 1) + newexps(pu, m2, 2)] += c1 * c2
        end
    end
    r
end

function *{T,U}(s::T, p::MPoly{U})
    r = zero(MPoly{promote_type(T,U)}, vars=vars(p))
    for (m, c) in p
        r[m] = s * c
    end
    r
end

*(p::MPoly, s) = s * p

function Base.inv(p)
    q = copy(p)
    if length(q.terms) != 1
        throw(ArgumentError("Can only take inverse of a monomial"))
    end
    term = collect(keys(q.terms))[1]
    val = collect(values(q.terms))[1]
    delete!(q.terms, term)
    q[-term] = val
    q
end

function ^(p::MPoly, power::Integer)
    if power < 0
        return inv(p)^abs(power)
    elseif power == 0
        return one(p)
    elseif power == 1
        return p
    else
        f, r = divrem(power, 2)
        return p^(f+r) * p^f
    end
end

function /{T,U}(p::MPoly{T}, s::U)
    r = zero(MPoly{promote_type(T,U)}, vars=vars(p))
    for (m, c) in p
        r[m] = c/s
    end
    r
end

function conj(p::MPoly)
    r = zero(p)
    for (m, c) in p
        r[m] = conj(c)
    end
    r
end

function real{T<:Real}(p::MPoly{Complex{T}})
    r = zero(MPoly{T}, vars=vars(p))
    for (m, c) in p
        r[m] = real(c)
    end
    r
end

real{T<:Real}(p::MPoly{T}) = p

function imag{T<:Real}(p::MPoly{Complex{T}})
    r = zero(MPoly{T}, vars=vars(p))
    for (m, c) in p
        r[m] = imag(c)
    end
    r
end

imag{T<:Real}(p::MPoly{T}) = zero(p)


text:
using AnyBoost
using Plots
using DataFrames
include("xgb.jl")

function check(args, N::Int, types...)
    @assert length(args) == N
    @assert length(types) == N
    for i in 1:N
        @assert isa(args[i], types[i])
    end
end

function onerun(trdata, tedata, l, B, O; algo=:fw, args=([1., 2.]), kwargs=Dict())
    if algo === :xg
        check(args, 0)
        return xgb(trdata, tedata, l, args...; kwargs...)
    end

    prob = Problem(trdata.X, trdata.y, l, B, O)

    time = @elapsed if algo === :fw
        check(args, 1, Vector)
        ens = ADFWPath(prob, args...; kwargs...)
    elseif algo === :gb
        check(args, 1, Real)
        ens = GradientBoost(prob, args...; kwargs...)
    elseif algo === :sw
        check(args, 1, Real)
        ens = Stagewise(prob, args...; kwargs...)
    elseif algo === :nn
        check(args, 2, Vector, Int)
        ens = NNPath(prob, args...; kwargs...)
    else
        error("Unknown Algorithm: ", algo)
    end

    norms, cards = normcardseq(ens)

    trerrs = trainerrseq(ens)/trdata.n
    teerrs = testerrseq(ens, tedata.X, tedata.y)/tedata.n
    bayes = l(tedata.y, tedata.signal)/tedata.n

    al = altloss(l)
    altrerrs = trainerrseq(ens, al)/trdata.n
    alteerrs = testerrseq(ens, tedata.X, tedata.y, al)/tedata.n
    albayes = al(tedata.y, tedata.signal)/tedata.n

    cards, norms, trerrs, teerrs, bayes, altrerrs, alteerrs, albayes, time
end


text:
struct PosteriorPsi
   beta::Distributions.Beta
   emperical::Vector{Float64}
   psi::Float64
end

PosteriorPsi() = PosteriorPsi( Beta(1,1), Float64[], 0.0 )

# create one PosteriorPsi object given psi, N, and priors α and β
# Posterior                      Likelihood         Prior
# Beta( Ψ*N + α, (1-Ψ)*N + β ) ∝ Binomial( Ψ, N ) * Beta( α, β )
function PosteriorPsi( Ψ::Float64, N::Float64, α=1.0, β=1.0; size=1000 )
   # Posterior Beta is k + α, n-k + β
   beta = Beta( Ψ*N + α, (1-Ψ)*N + β )
   emperical = rand(beta, size)
   psi = mean(beta)
   PosteriorPsi( beta, emperical, Ψ )
end

# combine multiple beta posteriors for replicates into one
# if unpaired, we don't need to maintain sample1_rep1 to sample2_rep1
# so we can re-sample from the new joint beta
function PosteriorPsi( set::Vector{PosteriorPsi}; paired=false, size=1000 )
   emperical = vcat( map( x->x.emperical, set )... )
   beta = fit(Beta, emperical)
   psi  = mean(beta)
   if !paired
      emperical = rand(beta, size)
   end
   PosteriorPsi( beta, emperical, psi )
end

# Calculates the probability, P( a-b > amt)
function probability( a::PosteriorPsi, b::PosteriorPsi; amt=0.0 )
   success = 0
   total = min(length(a.emperical),length(b.emperical))
   for i in 1:total
      success += a.emperical[i] - b.emperical[i] > amt ? 1 : 0
   end
   @fastmath success / total
end

function open_stream( filename )
   fopen = open( filename, "r" )
   if isgzipped( filename )
      stream = ZlibInflateInputStream( fopen, reset_on_end=true )
   else
      stream = BufferedStreams.BufferedInputStream( fopen )
   end
   stream
end

function open_streams( files::Vector{String} )
   buf = Vector{BufferedStreams.BufferedInputStream}(length(files))
   for i in 1:length(files)
      buf[i] = open_stream( files[i] )
   end
   buf
end

parse_float_omit_text( str::S, header::String ) where {S <: AbstractString} = str != "NA" && str != header ? parse(Float64, str) : 0.0

function parse_psi_line( line::String; min_num=5, size=1000 )
   res  = split( line, '\t' )
   psi  = parse_float_omit_text( res[6], "Psi" )
   num  = parse_float_omit_text( res[9], "Total_Reads" )
   if psi < 0 || num < min_num
      post = PosteriorPsi()
      bool = false
   else
      post = PosteriorPsi( psi, num, size=size )
      bool = true
   end
   res,post,bool
end

mutable struct PosteriorEvent
   event::Vector{SubString{String}}
   complexity::String
   a::PosteriorPsi
   b::PosteriorPsi
end

parse_complexity( c::S ) where {S <: AbstractString} = split( c, COMPLEX_CHAR, keep=false )[1] |> x->parse(Int,x)

function process_psi_line( streams::Vector{BufferedStreams.BufferedInputStream}; min_reads=5, size=1000 )
   postvec = Vector{PosteriorPsi}()
   event   = split( "", "" )
   complex = 0
   entropy = 0.0
   i = 1
   while i <= length(streams)
      line = readline( streams[i] )
      i += 1
      if line != ""
         par,post,isok = parse_psi_line( line, min_num=min_reads, size=size )
         event = par[1:5]
         event[5] == "BS" && (i -= 1; continue)
         !isok && continue
         push!( postvec, post )
         parcomplex = parse_complexity( par[10] )
         parentropy = parse_float_omit_text( par[11], "Entropy" )
         complex = parcomplex > complex ? parcomplex : complex
         entropy = parentropy > entropy ? parentropy : entropy
      end
   end
   event,complex,entropy,postvec
end

function process_psi_files( outfile, a::Vector{BufferedStreams.BufferedInputStream},
                                     b::Vector{BufferedStreams.BufferedInputStream};
                                     min_samp=1, min_reads=5, amt=0.0, size=1000 )
   io = open( outfile, "w" )
   stream = ZlibDeflateOutputStream( io )
   output_diff_header( stream )
   while true # go through all lines until we hit eof
      a_event,a_complex,a_entropy,a_post = process_psi_line( a, min_reads=min_reads, size=size )
      b_event,b_complex,b_entropy,b_post = process_psi_line( b, min_reads=min_reads, size=size )
      if a_event == split( "", "" ) || b_event == split( "", "" )
         break # eof
      end
      complex = a_complex > b_complex ? a_complex : b_complex
      entropy = a_entropy > b_entropy ? a_entropy : b_entropy
      @assert( a_event == b_event, "Incorrect events matched!!" )
      if length(a_post) >= min_samp && length(b_post) >= min_samp
         a_post   = PosteriorPsi( a_post ) # fit new posterior
         b_post   = PosteriorPsi( b_post )
         fwdprob  = probability( a_post, b_post, amt=amt )
         revprob  = probability( b_post, a_post, amt=amt )
         prob     = fwdprob > revprob ? fwdprob : revprob
         psi_a    = mean(a_post.beta)
         psi_b    = mean(b_post.beta)
         deltapsi = psi_a - psi_b
         output_diff( stream, a_event, complex, entropy, psi_a, psi_b, deltapsi, prob )
      end
   end
   close(stream)
   close(io)
end


text:
using Documenter
using Statistics
using TimeDag

DocMeta.setdocmeta!(TimeDag, :DocTestSetup, :(using TimeDag); recursive=true)

# This is inspired by the following, to prevent a `gksqt` process blocking (and needing
# manual termination) for every plot when generating documentation on MacOS.
# https://discourse.julialang.org/t/deactivate-plot-display-to-avoid-need-for-x-server/19359/2
# Another choice could be "100":
# https://discourse.julialang.org/t/generation-of-documentation-fails-qt-qpa-xcb-could-not-connect-to-display/60988
# Empirically, both choices seem to work.
withenv("GKSwstype" => "nul") do
    makedocs(;
        modules=[TimeDag],
        authors="Invenia Technical Computing Corporation",
        repo="https://github.com/invenia/TimeDag.jl/blob/{commit}{path}#{line}",
        sitename="TimeDag.jl",
        format=Documenter.HTML(;
            prettyurls=get(ENV, "CI", "false") == "true",
            canonical="https://invenia.github.io/TimeDag.jl",
            assets=String[],
        ),
        pages=[
            "Home" => "index.md",
            "concepts.md",
            "examples.md",
            "Reference" => [
                "reference/fundamentals.md",
                "Node ops" => [
                    "reference/align.md",
                    "reference/arithmetic.md",
                    "reference/online_windowed.md",
                    "reference/misc_ops.md",
                    "reference/sources.md",
                ],
                "reference/creating_ops.md",
                "reference/internals.md",
            ],
        ],
        checkdocs=:exports,
        strict=true,
    )
end

deploydocs(; repo="github.com/invenia/TimeDag.jl", devbranch="main")


text:
using GeoData, RasterDataSources, Test, Dates, NCDatasets, ArchGDAL

# Too big to test on CI
# if !haskey(ENV, "CI")
#     @testset "load WorldClim Weather" begin
#         # Weather time-series
#         dates = (Date(2001), Date(2002))
#         ser = series(WorldClim{Weather}, (:prec,); date=dates)
#         ser[Date(2001, 1)][:prec]
#         A = geoarray(WorldClim{Weather}, :prec; date=DateTime(2001, 05), mappedcrs=EPSG(4326))
#     end
# end

@testset "load WorldClim Climate" begin
    # Weather time-series
    ser = GeoSeries(WorldClim{Climate}, :prec; res="10m", month=Jan:March, mappedcrs=EPSG(4326))
    # Select Australia, using regular lat/lon selectors
    A = ser[month=Jan]
    @test A isa GeoArray
    A[Y(Between(-10, -45)), X(Between(110, 160))]
    st = GeoStack(WorldClim{Climate}, (:prec, :tmax); month=1)
    st[:prec]
    @test st isa GeoStack
end

@testset "load WorldClim BioClim" begin
    A = GeoArray(WorldClim{BioClim}, :Bio_1; mappedcrs=EPSG(4326))
    A[Y(Between(-10, -45)), X(Between(110, 160))]
    @test A isa GeoArray
    st = stack(WorldClim{BioClim}, (1, 2))
    st[:bio1]
    @test st isa GeoStack
    @test A isa GeoArray
end

@testset "load CHELSA BioClim" begin
    A = GeoArray(CHELSA{BioClim}, 1; mappedcrs=EPSG(4326))
    @test GeoData.name(A) == :bio1
    st = GeoStack(CHELSA{BioClim}, (:bio1, :BIO2))
    @test keys(st) == (:bio1, :bio2)
    @test A isa GeoArray
    @test st isa GeoStack
    @test st[:bio2] isa GeoArray
end

@testset "load EarthEnv HabitatHeterogeneity" begin
    A = GeoArray(EarthEnv{HabitatHeterogeneity}, :cv; mappedcrs=EPSG(4326))
    A[Y(Between(-10, -45)), X(Between(110, 160))]
    st = GeoStack(EarthEnv{HabitatHeterogeneity}, (:cv, :evenness))
    @test A isa GeoArray
    @test st isa GeoStack
    @test st[:evenness] isa GeoArray
end

@testset "load EarthEnv LandCover" begin
    A = GeoArray(EarthEnv{LandCover}, 2; mappedcrs=EPSG(4326))
    @test GeoData.name(A) == :evergreen_broadleaf_trees
    A[Y(Between(-10, -45)), X(Between(110, 160))]
    @test A isa GeoArray
    st = GeoStack(EarthEnv{LandCover}, (:evergreen_broadleaf_trees, :deciduous_broadleaf_trees); mappedcrs=EPSG(4326))
    keys(st) = (:evergreen_broadleaf_trees, :deciduous_broadleaf_trees)
    @test st isa GeoStack
end

@testset "load ALWB" begin
    A = GeoArray(ALWB{Deciles,Day}, :rain_day; date=DateTime(2019, 10, 19))
    @test crs(A) == EPSG(4326)
    A = GeoArray(ALWB{Values,Day}, :ss_pct; date=DateTime(2019, 10, 19))
    @test crs(A) == EPSG(4326)
    st = GeoStack(ALWB{Values,Day}, (:s0_pct, :ss_pct); date=DateTime(2019, 10, 19))
    @test crs(st) == EPSG(4326)
    @test crs(st[:s0_pct]) == EPSG(4326)
    dates = DateTime(2019, 10, 19), DateTime(2021, 11, 20)
    s = GeoSeries(ALWB{Values,Day}, (:s0_pct, :ss_pct); date=dates)
    s[1]
    @test A isa GeoArray
    @test st isa GeoStack
    @test s isa GeoSeries
end

# Obscure .Z format may not work on windows
if Sys.islinux()
    @testset "load AWAP" begin
        A = GeoArray(AWAP, :rainfall; date=DateTime(2019, 10, 19))
        @test crs(A) == EPSG(4326)
        st = GeoStack(AWAP; date=DateTime(2019, 10, 19), resize=crop)
        @test crs(st) == EPSG(4326)
        dates = DateTime(2019, 09, 19), DateTime(2019, 11, 19)
        s = GeoSeries(AWAP; date=dates, resize=crop)
        # s = GeoSeries(AWAP; date=dates, resize=resample, crs=EPSG(4326)) TODO: all the same
        # s = GeoSeries(AWAP; date=dates, resize=extend) TODO: this is slow !!!
        @test crs(s[1][:solar]) == EPSG(4326)
        @test A isa GeoArray
        @test st isa GeoStack
        @test s isa GeoSeries
    end
end


text:
using BrainWave, Statistics, Test

sampRate = 1000
bsDur = 0.1
sigDur = 1-0.1

bsAmp = 0.5
sigAmp = 1

nSampBs = round(Int, bsDur*sampRate)
nSampSig = round(Int, sigDur*sampRate)
nSampTot = nSampBs+nSampSig
bs = zeros(nSampBs) .+ bsAmp
sig = zeros(nSampSig) .+ sigAmp

recVals = vcat(bs, sig)'

epochDur=sigDur; preDur=bsDur; events=[1,2];
rec, evtTab = simulateRecording(nChans=1, dur=10, epochDur=epochDur, preDur=preDur, events=events, sampRate=sampRate)
segs, nSegs = segment(rec, evtTab, -preDur, epochDur, sampRate, eventsList=[1, 2], eventsLabelsList=["cnd1", "cnd2"])

nSegsCnd1 = size(segs["cnd1"])[3]
nSegsCnd2 = size(segs["cnd2"])[3]

for i=1:nSegsCnd1
    segs["cnd1"][:, 1:nSampBs, i] .= bsAmp
    segs["cnd1"][:, nSampBs+1:nSampTot, i] .= sigAmp
end

for i=1:nSegsCnd2
    segs["cnd2"][:, 1:nSampBs, i] .= bsAmp
    segs["cnd2"][:, nSampBs+1:nSampTot, i] .= sigAmp
end

baselineCorrect!(segs, -bsDur, bsDur, sampRate)

@test segs["cnd1"][:,nSampBs+1,1][1] .== 0.5

for i=1:nSegsCnd1
    @test mean(segs["cnd1"][:,nSampBs+1:nSampTot,i]) == 0.5
end

for i=1:nSegsCnd2
    @test mean(segs["cnd2"][:,nSampBs+1:nSampTot,i]) == 0.5
end

#test case with `bsStart` and `bsEnd`

for i=1:nSegsCnd1
    segs["cnd1"][:, 1:50, i] .= bsAmp-0.5
    segs["cnd1"][:, 51:101, i] .= bsAmp
    segs["cnd1"][:, 102:nSampTot, i] .= sigAmp
end

for i=1:nSegsCnd2
    segs["cnd2"][:, 1:50, i] .= bsAmp-0.5
    segs["cnd2"][:, 51:101, i] .= bsAmp
    segs["cnd2"][:, 102:nSampTot, i] .= sigAmp
end

baselineCorrect!(segs, -bsDur/2, 0, preDur, sampRate)

@test segs["cnd1"][:,102,1][1] .== 0.5

for i=1:nSegsCnd1
    @test mean(segs["cnd1"][:,102:nSampTot,i]) == 0.5
end

for i=1:nSegsCnd2
    @test mean(segs["cnd2"][:,102:nSampTot,i]) == 0.5
end

@test_throws(ErrorException, baselineCorrect!(segs, 0, 0, preDur, sampRate))
@test_throws(ErrorException, baselineCorrect!(segs, -0.05, -0.06, preDur, sampRate))


text:
"""
plotnsold(half="no",c=.5,maxit=10,level=5; bigtitle="", table=false)

Makes the plots for the SIAM Review paper.

Makes convergence rate tables for Float16 computations.

Newton's Method in Mixed-Precision.

"""
function plotnsold(
    half = "no",
    c = 0.5,
    maxit = 10,
    level = 5;
    bigtitle = "",
    table = false,
)
    fmtplot = ("k-", "k--", "k-.", "k-.", "k>:")
    Fdata = zeros(4, level, maxit + 1)
    dout = zeros(maxit + 1, level)
    HeqAnnotations = fillTitles(half, level)
    titles4heq = HeqAnnotations.titles4heq
    legenddata = HeqAnnotations.legenddata
    gxlabel = HeqAnnotations.gxlabel
    gylabel = HeqAnnotations.gylabel
    fillFdata!(Fdata, level, half)
    #
    # Draw the windows for analytic/difference Jacobians.
    # The 64-32 bit comparisons are on a 2x2 plot.
    # The 16 bit results are on a 1x2 plot.
    #
    ipmax = 4
    if half == "yes"
        ipmax = 2
    end
    if table == false
        for ip = 1:ipmax
            if half == "no"
                subplot(2, 2, ip)
            else
                subplot(1, 2, ip)
            end
            for ir = 1:level
                dout[:, ir] = Fdata[ip, ir, :] ./ Fdata[ip, ir, 1]
            end
            for ir = 1:level
                semilogy(dout[:, ir], fmtplot[ir])
            end
            if ip == 1
                legend(legenddata, fontsize = "6")
            end
            xlabel(gxlabel)
            ylabel(gylabel)
            aymin = ymin(c, half)
            axis([0.0, maxit, aymin, 1.0])
            title(titles4heq[ip])
        end
        PyPlot.tight_layout()
        PyPlot.suptitle(bigtitle)
    else
        if half == "yes"
            maketable(Fdata, level, maxit)
        end
    end
end

function maketable(Fdata, level, maxit)
    headers = ["n", "1024", "2048", "4096", "8192", "16384"]
    formats = "%d & %9.5e & %9.5e & %9.5e & %9.5e & %9.5e"
    Rates = zeros(maxit, 6)
    Rates[:, 2:6] = transpose(Fdata[1, 1:level, 2:maxit+1] ./ Fdata[1, 1:level, 1:maxit])
    Rates[:, 1] = 1:maxit
    fprintTeX(headers, formats, Rates)
end

function fillTitles(half, level)
    if half == "no"
        #   titles4heq=("Double precision, analytic Jacobian",
        #     "Double precision, finite difference Jacobian",
        #     "Single precision, analytic Jacobian",
        #     "Single precision, finite difference Jacobian")
        titles4heq = ("Double, analytic ", "Double, FD", "Single, analytic", "Single, FD")
    else
        #   titles4heq=("Half, analytic",
        #     "Half, FD")
        titles4heq = (
            "Half precision, analytic Jacobian",
            "Half precision, finite difference Jacobian",
        )
    end
    legenddata = Array{String,1}(undef, level)
    for il = 1:level
        legenddata[il] = string("N=", string(512 * 2^il))
    end
    #legenddata=("N=1024","N=2048","N=4096","N=8192","N=16384")
    gxlabel = "Nonlinear Iterations"
    gylabel = L"$|| F ||/||F_0||$"
    return (
        titles4heq = titles4heq,
        legenddata = legenddata,
        gxlabel = gxlabel,
        gylabel = gylabel,
    )
end

function fillFdata!(Fdata, level, half)
    #
    # Read the data files and get organized.
    #
    loadme(fmt::String) = @eval @load($fmt)
    for ir = 1:level
        gridlev = 512 * (2^ir)
        #mixedpfile = string("paper", string(gridlev), ".jld")
        mixedpfile=string("paper",string(gridlev),".jld2")
        loadme(mixedpfile)
        if half == "no"
            Fdata[1, ir, :] .= fout64.exactout.ithist
            Fdata[2, ir, :] .= fout64.fdout.ithist
            Fdata[3, ir, :] .= fout32.exactout.ithist
            Fdata[4, ir, :] .= fout32.fdout.ithist
        else
            Fdata[1, ir, :] .= fout16.exactout.ithist
            Fdata[2, ir, :] .= fout16.fdout.ithist
        end
    end
end

function ymin(c, half)
    ymin = 1.e-15
    if half == "yes"
        if c == 1
            ymin = 1.e-6
        elseif c == 99
            ymin = 1.e-10
        end
    end
    return ymin
end

max_stars_repo_name: ctkelley/MPResults


text:
using SortingAlgorithms

struct PopulationGrouping{PersonIdx <: Integer}
  person_ids::Vector{PersonIdx}
  group_ptrs::Vector{PersonIdx}
end

function show(io::IO, g::PopulationGrouping)
  num_groups = numgroups(g)
  print(io, "PopulationGrouping with ", numgroups, " groups")
  for group_id in 1:num_groups
    group = getgroup(g, group_id)
    println("Group ", group_id, " of size ", length(group))
  end
end

function make_group_ptrs(group_ids::AbstractArray{T} where T<:Integer, num_groups::Integer, PersonIdx::DataType=UInt32)
  N = length(group_ids)
  # sorting forces individuals of the same group to be together on the list
  sorted_ids = Vector{PersonIdx}(undef, N) .= sortperm(group_ids; alg=RadixSort) #alloc x3

  group_ptrs = zeros(PersonIdx, num_groups + 1) #alloc
  for group_id in group_ids
    group_ptrs[group_id+1] += 1 # size of a group in each entry
  end

  # group i starts at index group_ptrs[i], ends at index group_ptrs[i+1]-1
  group_ptrs[1] = 1
  for i in 1:num_groups
    group_ptrs[i+1] = group_ptrs[i] + group_ptrs[i+1]
  end
  sorted_ids, group_ptrs
end

PopulationGrouping(group_ids, num_groups::Integer, PersonIdx::DataType=UInt32) =
  PopulationGrouping(make_group_ptrs(group_ids, num_groups, PersonIdx)...)

numgroups(g::PopulationGrouping) = length(g.group_ptrs) - 1
getgroup(g::PopulationGrouping, group_id::Integer) = @view g.person_ids[g.group_ptrs[group_id]:(g.group_ptrs[group_id+1]-1)]
groupsizes!(sizes::AbstractVector{T} where T<:Integer, g::PopulationGrouping) = (sizes .= (@view g.group_ptrs[2:end]) .- (@view g.group_ptrs[1:end-1]))
groupsizes(g::PopulationGrouping) = groupsizes!(Vector{PersonIdx}(undef, numgroups(g)), g)


text:
import DrWatson: quickactivate
quickactivate(@__DIR__, "Chemostat_Heerden2013")

# ------------------------------------------------------------------
@time begin
    using CSV
    using DataFrames

    # Run add https://github.com/josePereiro/Chemostat_Heerden2013.jl in the Julia Pkg REPL to install the
    # package, then you must activate the package enviroment (see README)
    import Chemostat_Heerden2013
    const ChH = Chemostat_Heerden2013
    const Hd = ChH.HeerdenData
end

## ------------------------------------------------------------------
# Load data
# Data taken from Heerden2013 https://doi.org/10.1186/1475-2859-12-80.

# Table 2 Steady state results from continuous fermentations
orig_data = CSV.read(Hd.HEERDEN_CONT_CUL_DATA_ORIG_FILE,  DataFrame;
    delim = "\t", comment = "#");
for n in names(orig_data)
    orig_data[!, n] = Vector{Union{String, Float64}}(orig_data[!, n])
    orig_data[3:end, n] .= parse.(Float64, orig_data[3:end, n])
end
first(orig_data, 6)

## ------------------------------------------------------------------
# Convert data

# Converting data c(g/L) / MM (g/mol) -> c(mM)
conv_data = DataFrame(orig_data)
conv_data[3:end, :cGLC] .= orig_data[3:end, :cGLC] .* 1e3/ 180;
conv_data[3:end, :sGLC] .= orig_data[3:end, :sGLC] .* 1e3/ 180;
conv_data[3:end, :sSUCC] .= orig_data[3:end, :sSUCC] .* 1e3/ 118;
conv_data[3:end, :sAC] .= orig_data[3:end, :sAC] .* 1e3/ 60.02;
conv_data[3:end, :sFORM] .= orig_data[3:end, :sFORM] .* 1e3/ 46.03;
conv_data[3:end, :sMAL] .= orig_data[3:end, :sMAL] .* 1e3/ 134.09;
conv_data[2, [:cGLC, :sGLC,:sSUCC,:sAC,:sFORM,:sMAL]] .= "mM"
first(conv_data, 6)

## ------------------------------------------------------------------
# Saving
CSV.write(Hd.HEERDEN_CONT_CUL_DATA_CONV_FILE, conv_data, delim = "\t");
println(relpath(Hd.HEERDEN_CONT_CUL_DATA_CONV_FILE), " created!!!")


text:
#script for solving a 2D BVP of the form
# -a0(x)*\Delta u(x,y)+a1(x)*u(x,y)=f(x,y) for (x,y) ∈ Ω:= (0,1)^2
# Here a0(x,y)=1, a1(x,y)=0, f(x,y)=-8*π^2*sin(2*π*x)*sin(2*π*y)
# with Dirichlet boundary conditions u(x,y)=0 for (x,y) ∈ ∂Ω
include("../src/Solver2D.jl")
include("../src/SplinePlotting.jl")

a0(x,y) = 1
a1(x,y) = 0
f(x,y) = 8*π^2*sin(2*π*x)*sin(2*π*y)
exact_sol(x,y) = sin(2*π*x)*sin(2*π*y)
deriv_exact_sol(x,y) = [2*π*cos(2*π*x)*sin(2*π*y), 2*π*sin(2*π*x)*cos(2*π*y)]

#Define the boundary conditions
u_bound(x,y) = 0
bound_left = Boundary2D("Dirichlet", "Left", u_bound)
bound_right = Boundary2D("Dirichlet", "Right", u_bound)
bound_down = Boundary2D("Dirichlet", "Down", u_bound)
bound_up = Boundary2D("Dirichlet", "Up", u_bound)
bound_all = [bound_left, bound_right, bound_down, bound_up]

#Define the domain geometry
cornerLowerLeft = [0., 0.]
lengthx = 1.
lengthy = 1.
uPtRight = 1.
degP = [3, 3]
numSubdiv = 80
nrb = nrbsquare(cornerLowerLeft, lengthx, lengthy, degP, numSubdiv)
#Plots.surface(reuse = false)
#plt = nrbctrlplot(nrb)
#plt = nrbkntplot(nrb)
#display(plt)
IEN, elemVertex = makeIEN(nrb.knots, numSubdiv^2, degP)
IGAmesh = genMesh(nrb)
#plotBasisParam(IGAmesh)

gauss_rule = [genGaussLegendre(degP[1]+1), genGaussLegendre(degP[2]+1)]
stiff = assemble_stiff2D(IGAmesh, a0, gauss_rule)
rhs = assemble_rhs2D(IGAmesh, f, gauss_rule)
bcdof_all, elem_all = classifyBoundary2D(IGAmesh)
lhs, rhs = applyBCnurbs(IGAmesh, bound_all, stiff, rhs, bcdof_all, elem_all, gauss_rule)
@time sol0 = lhs\rhs
plotSol(IGAmesh, sol0, "Poisson2DSquare")
plotSolError(IGAmesh, sol0, exact_sol, "Poisson2DSquareError")
relL2Err, relH1Err = compErrorNorm(IGAmesh, sol0, exact_sol, deriv_exact_sol, a0, gauss_rule)
println("Relative L2-norm error is $relL2Err")
println("Relative energy-norm error is $relH1Err")
print("Done!")