text:
#isdefined(Main, :CartesianBoxes) ||
#    include(joinpath("..", "src", "CartesianBoxes.jl"))
#push!(LOAD_PATH, normpath(joinpath(dirname(@__FILE__), "../src")))

module CartesianBoxesBenchmarks

using CartesianBoxes
using BenchmarkTools

# Deal with compatibility issues.
using Compat
@static if isdefined(Base, :CartesianIndices)
    import Base: CartesianIndices
    const CARTESIAN_REGIONS = (CartesianIndices, CartesianBox)
else
    import Base: CartesianRange
    import Compat: CartesianIndices
    const CARTESIAN_REGIONS = (CartesianIndices, CartesianRange, CartesianBox)
end
const CartesianRegions = Union{CARTESIAN_REGIONS...}
@static if isdefined(Base, :axes)
    import Base: axes
else
    import Base: indices
    const axes = indices
end

function add!(dst::AbstractArray{T,N},
              A::AbstractArray{T,N},
              B::AbstractArray{T,N}) where {T,N}
    inds = axes(dst)
    @assert axes(A) == inds
    @assert axes(B) == inds
    @inbounds @simd for i in eachindex(dst, A, B)
        dst[i] = A[i] + B[i]
    end
    return dst
end

function add!(::typeof(eachindex),
              dst::AbstractArray{T,N},
              A::AbstractArray{T,N},
              B::AbstractArray{T,N}) where {T,N}
    inds = axes(dst)
    @assert axes(A) == inds
    @assert axes(B) == inds
    @inbounds @simd for i in eachindex(dst, A, B)
        dst[i] = A[i] + B[i]
    end
    return dst
end

function add!(::Type{R},
              dst::AbstractArray{T,N},
              A::AbstractArray{T,N},
              B::AbstractArray{T,N}) where {R<:CartesianRegions,T,N}
    inds = axes(dst)
    @assert axes(A) == inds
    @assert axes(B) == inds
    @inbounds @simd for i in R(inds)
        dst[i] = A[i] + B[i]
    end
    return dst
end

function inner(A::AbstractArray{T,N},
               B::AbstractArray{T,N}) where {T<:AbstractFloat,N}
    @assert axes(A) == axes(B)
    s = zero(T)
    @inbounds @simd for i in eachindex(A, B)
        s += A[i]*B[i]
    end
    return s
end

function inner(::typeof(eachindex),
               A::AbstractArray{T,N},
               B::AbstractArray{T,N}) where {T<:AbstractFloat,N}
    @assert axes(A) == axes(B)
    s = zero(T)
    @inbounds @simd for i in eachindex(A, B)
        s += A[i]*B[i]
    end
    return s
end

function inner(::Type{R},
               A::AbstractArray{T,N},
               B::AbstractArray{T,N}) where {R<:CartesianRegions,
                                             T<:AbstractFloat,N}
    inds = axes(A)
    @assert axes(B) == inds
    s = zero(T)
    @inbounds @simd for i in R(inds)
        s += A[i]*B[i]
    end
    return s
end

T = Float32
dims = (20,30,40)
A = rand(T,dims...)
B = rand(T,dims...)
C1 = Array{T}(undef, dims)
C2 = Array{T}(undef, dims)
C3 = Array{T}(undef, dims)

brief(::typeof(eachindex)) = "eachindex";
brief(::Type{CartesianBox}) = "CartesianBox";
@static if isdefined(Base, :CartesianRange)
    brief(::Type{CartesianRange}) = "CartesianRange";
end
@static if isdefined(Base, :CartesianIndices)
    brief(::Type{CartesianIndices}) = "CartesianIndices";
end

const VERBOSE = false
add!(C1, A, B)

println("Addition of arrays of size $dims ($(prod(dims)) operations)")
for R in (eachindex, CARTESIAN_REGIONS...)
    flg = (R == last(CARTESIAN_REGIONS))
    pfx1 = (flg ? " └─ " : " ├─ ")
    pfx2 = (flg ? "    " : " │  ")
    println(pfx1, "Computations with `$(brief(R))`:")

    fill!(C2, 0)
    add!(R, C2, A, B)
    println(pfx2, " ├─ diff. ∈ ", extrema(C1 - C2))

    fill!(C3, 0)
    if VERBOSE
        t = @benchmark add!($R, $C3, $A, $B)
        show(stdout, MIME"text/plain"(), t)
    else
        print(pfx2, " └─ @btime:")
        @btime add!($R, $C3, $A, $B)
    end
end

w = inner(A, B)
println("Inner of arrays of size $dims ($(2*prod(dims)) operations)")
for R in (eachindex, CARTESIAN_REGIONS...)
    flg = (R == last(CARTESIAN_REGIONS))
    pfx1 = (flg ? " └─ " : " ├─ ")
    pfx2 = (flg ? "    " : " │  ")
    println(pfx1, "Computations with  `$(brief(R))`:")

    wp = inner(R, A, B)
    println(pfx2, " ├─ result = $wp ($w)")

    if VERBOSE
        t = @benchmark inner($R, $A, $B)
        show(stdout, MIME"text/plain"(), t)
    else
        print(pfx2, " └─ @btime:")
        @btime inner($R, $A, $B)
    end
end

#show(STDOUT, MIME"text/plain"(), @benchmark $p(dat, a, img, b))

end # module


text:
__precompile__()

# HamNNNNclassifier reporting functions

"""
Print a report detailing the data file.
"""
function printdatafiledescription(dataFile)
	printheader()
	print_with_color(:blue, "Data File:\n"; bold=true)
	# print file path
	println(dataFile)
	# print file type
	format = getorangefileformat(dataFile)
	println("Orange file format: ", format, "\n")
	# get the data file contents
	codes, names, dt = readdatafile(dataFile)
	# @show codes names
	#for the class attribute, list each class and the number of cases, % of total, as well as the total
	printclassdescription(codes, names, dt)
	printattributedescriptions(codes, names, dt)
end

function printheader()
	# print date, time, software version, platform, file
	print("\n\n")
	f1 = @sprintf "%18s ; HamNNNNclassifier version %5s" Dates.format(now(), "yyyy-m-d HH:MM:SS") version
	print_with_color(:blue, f1, "\n"; bold=true)
	println("Julia version: ", VERSION)
	println("On a Mac ", Sys.cpu_info()[1].model, " '", ccall(:jl_get_cpu_name, Ref{String}, ()), "', \nRunning ", Sys.MACHINE)
	run(`sw_vers`) # shows os version
	println()
end

function printclassdescription(codes, names, table)
	#for the class attribute, list each class and the number of cases, % of total, as well as the total
	i, classesVector = extractclassattribute(table, codes)
	classAttributeName = names[i]
	println("Class Attribute: ", classAttributeName)
	casesNo = length(classesVector)
	# extract unique values
	# as it is easy to print out the strings for the levels of a categorical array, do the conversion.
	cv = categorical(classesVector)
	# levels is simply the array of uniques, for a categorical array
	levels = CategoricalArrays.levels(cv)
	# uniques works well as keys for getting values out of a countmap dict
	uniques = unique(cv)
	# create a dict of uniques and their counts
	cm = countmap(classesVector)
	# @show cm levels uniques
	print_with_color(:blue, "\nClass        Number of cases      Percent of total\n"; bold=true)
	foreach((x,y) -> print(rpad(y,20,),rpad(get(cm,x,""),20,),signif(get(cm,x,"")*100/casesNo,4),"\n"), uniques, levels)
	print_with_color(:bold, "total: ", casesNo, " cases\n")
end

function printattributedescriptions(codes, names, table)
	# for each attribute, print its name, type, number of uniques, number and % of missing values;
	# for continuous attributes, print the minimum, maximum, mean, and median
	# start with a header
	println()
	print_with_color(:blue, " Index   Label             Type   Missing      %   Uniques   Min   Max   Mean   Median    Unique Values\n"; bold=true)
	# initialize a counter to hold the total number of missing values
	missingTotal = 0
	# initialize an array for new attribute codes
	newCodes = Array{String}(length(names))
	for i in 1:length(names)
		# print the index and the attribute name
		name = rpad(names[i], 18, )
		print(lpad(i,6,),"   ", name)
		# take out null values to get type purity for type dispatching
		d = dropnull(table[i])
		# print the orange format code; if missing, print "C" for continuous (ie numeric)
		# and "D" for discrete (ie categorical) attributes
		typeCode = generateattributecode(d, codes[i])
		# convert "C" to "D" for integer attributes in range 0 to 9
		if typeCode == "C"
			typeCode = changeattributetype(d)
		end
		print(rpad(typeCode,7,))
		# add the new codes to the existing codes, in a newCodes array
		newCodes[i] = typeCode
		# print number and percent of missing values (nulls)
		missing = length(table[i]) - length(d)
		missingTotal += missing
		missingPerCent = roundpercentage(missing * 100 / length(table[i]),2)
		print(lpad(missing,7,), missingPerCent > 0.0 ? lpad(missingPerCent,7,) : "       ")
		# print number of unique values
		print(lpad(length(unique(d)),10,))
		# for continuous attributes, print the minimum, maximum, mean, and median; for discrete
		# attributes, print the number of unique values.
		if numericattribute(d) && typeCode == "C"
			((min, max),mean,median) = describecontinuousattribute(d)
			print(lpad(min,6,), lpad(max,6,), lpad(roundpercentage(mean,1),7,), lpad(roundpercentage(median,1),9,))
		else
			printuniquestring(d)
		end
		print("\n")
	end
	# print the total number of missing values and its percentage
	valuesTotal = prod(size(table))
	print_with_color(:bold, "Total missing values: ", missingTotal, " out of ", valuesTotal)
	print_with_color(:bold, " (", roundpercentage(missingTotal * 100 / valuesTotal,2), "%)\n\n")
	# print the counts for each attribute type and the total
	print_with_color(:blue, "Attribute Type   Count\n"; bold=true)
	# print(countmap(newCodes), "\n")
	for (key, value) in countmap(newCodes)
		println(rpad(key,14,), lpad(value,8,))
	end
	print_with_color(:bold, "Total count:        ", length(newCodes), "\n")
end

"""
Print a list of unique values for an attribute. If the list is too long, print just
its head and tail.
"""
function printuniquestring(d)
	uniquesNo, uniques = describediscreteattribute(d)
	print("                                ")
	# if uniquesNo is > 10, print just the first and last 5 values
	if uniquesNo > 10
		for i in 1:4
			print(uniques[i]," ")
		end
		print("... ")
		for i in uniquesNo - 3:uniquesNo
			print(uniques[i]," ")
		end
	else
		for i in 1:uniquesNo
			print(uniques[i]," ")
		end
	end
end

function printattributeranking(dataFile)
	labels, codes, table, D = generateattributeranking(dataFile)
	@show codes
	print_with_color(:blue, "Attributes in rank order of ability to separate classes:\n"; bold=true)
	print_with_color(:blue, "  Rank   Index   Label                   Type   RankValue   Optimized Slices\n"; bold=true)
	for i in indices(D,1)
		index = lpad(D[i,2], 7)
		label = rpad(labels[D[i,2]], 24)
		code = rpad(codes[D[i,2]], 4)
		rank = lpad(D[i,1], 12)
		println(lpad(i,7), index,"   ", label, code, rank)
	end
end

ext: jl


text:
module BoardGamesGUI

using Gtk
using Cairo
using BoardGames

# Write your package code here.

export
    draw,
    correspondingmove,
    mainwindow

include("functions.jl")
include("main.jl")
include("userstrategy.jl")

end


text:
cd(@__DIR__); include("setups/simpletree.jl")
gr(dpi = 200)

##
ind_eig = [5, 6, 11, 18, 21, 25, 38]
plot(layout = Plots.grid(1, 7), size = (1400, 500))
for i = 1:length(ind_eig)
    l = ind_eig[i]
    w = 𝚽[:, l]
    scatter!(X[:, 1], X[:, 2], marker_z = w, ms = 5, c = :viridis, subplot = i, mswidth = 0)
    plot!(frame = :none, cbar = false, grid = false, legend = false, subplot = i)
    plot!(title = latexstring("\\mathbf{\\phi}_{", l - 1, "}"), titlefont = 18, subplot = i)
end
plot!(subplot = 10, grid = false, frame = :none)
plt = current()
savefig(plt, "../figs/simpletree_f_rngwf_top21_ref_eigenvectors.png")


text:
struct BufferedFunction{TF,TB} <: Any
    f::TF
    buff::TB

    function BufferedFunction(f, buff)
        new{typeof(f), typeof(buff)}(f, buff)
    end
end

function (f::BufferedFunction)(x)
    return f.f(f.buff, x)
end


text:
mutable struct ConvergenceSimulation{SolType}
  solutions::Array{SolType}
  errors
  N
  auxdata
  𝒪est
  convergence_axis
end

function ConvergenceSimulation(solutions,convergence_axis;
                               auxdata=nothing,additional_errors=nothing)
  N = size(solutions,1)
  uEltype = eltype(solutions[1].u[1])
  errors = Dict() #Should add type information
  if isempty(solutions[1].errors)
    error("Errors dictionary is empty. No analytical solution set.")
  end
  for k in keys(solutions[1].errors)
    errors[k] = [mean(sol.errors[k]) for sol in solutions]
  end
  if additional_errors != nothing
    for k in keys(additional_errors)
      errors[k] = additional_errors[k]
    end
  end
  𝒪est = Dict((calc𝒪estimates(p) for p = pairs(errors)))
  #𝒪est = Dict(map(calc𝒪estimates,errors))
  𝒪esttmp = Dict() #Makes Dict of Any to be more compatible
  for (k,v) in 𝒪est
    if length(v)==1 push!(𝒪esttmp,Pair(k,v[1]))
    else push!(𝒪esttmp,Pair(k,v))
    end
  end
  𝒪est = 𝒪esttmp
  return(ConvergenceSimulation(solutions,errors,N,auxdata,𝒪est,convergence_axis))
end

function test_convergence(dts::AbstractArray,prob::Union{AbstractRODEProblem,AbstractSDEProblem},
                          alg;trajectories,save_everystep=true,timeseries_steps=1,
                          timeseries_errors=save_everystep,adaptive=false,
                          weak_timeseries_errors=false,weak_dense_errors=false,kwargs...)
  N = length(dts)
  ensemble_prob = EnsembleProblem(prob)
  _solutions = [solve(ensemble_prob,alg;dt=dts[i],save_everystep=save_everystep,
                      timeseries_steps=timeseries_steps,adaptive=adaptive,
                      timeseries_errors=timeseries_errors,trajectories=trajectories,
                      kwargs...) for i in 1:N]
  solutions = [DiffEqBase.calculate_ensemble_errors(sim;weak_timeseries_errors=weak_timeseries_errors,weak_dense_errors=weak_dense_errors) for sim in _solutions]
  auxdata = Dict("dts" =>  dts)
  # Now Calculate Weak Errors
  additional_errors = Dict()
  for k in keys(solutions[1].weak_errors)
    additional_errors[k] = [sol.weak_errors[k] for sol in solutions]
  end
  ConvergenceSimulation(solutions,dts,auxdata=auxdata,additional_errors=additional_errors)
end

function analyticless_test_convergence(dts::AbstractArray,
                          prob::Union{AbstractRODEProblem,AbstractSDEProblem},
                          alg,test_dt;trajectories=100,
                          save_everystep=true,timeseries_steps=1,
                          timeseries_errors=save_everystep,adaptive=false,
                          weak_timeseries_errors=false,weak_dense_errors=false,kwargs...)
  _solutions = []
  tmp_solutions = Array{Any}(undef,trajectories,length(dts))
  for j in 1:trajectories
    @info "Monte Carlo iteration: $j/$trajectories"
    t = prob.tspan[1]:test_dt:prob.tspan[2]
    if prob.noise_rate_prototype === nothing
      brownian_values = cumsum([[zeros(size(prob.u0))];[sqrt(test_dt)*randn(size(prob.u0)) for i in 1:length(t)-1]])
      brownian_values2 = cumsum([[zeros(size(prob.u0))];[sqrt(test_dt)*randn(size(prob.u0)) for i in 1:length(t)-1]])
    else
      brownian_values = cumsum([[zeros(size(prob.noise_rate_prototype,2))];[sqrt(test_dt)*randn(size(prob.noise_rate_prototype,2)) for i in 1:length(t)-1]])
      brownian_values2 = cumsum([[zeros(size(prob.noise_rate_prototype,2))];[sqrt(test_dt)*randn(size(prob.noise_rate_prototype,2)) for i in 1:length(t)-1]])
    end
    np = NoiseGrid(t,brownian_values,brownian_values2)
    _prob = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,
                       noise=np,
                       noise_rate_prototype=prob.noise_rate_prototype);
    true_sol = solve(_prob,alg;adaptive=adaptive,dt=test_dt);
    for i in 1:length(dts)
      sol = solve(_prob,alg;dt=dts[i],adaptive=adaptive);
      err_sol = appxtrue(sol,true_sol)
      tmp_solutions[j,i] = err_sol
    end
  end
  _solutions = [EnsembleSolution(tmp_solutions[:,i],0.0,true) for i in 1:length(dts)]
  solutions = [DiffEqBase.calculate_ensemble_errors(sim;weak_timeseries_errors=weak_timeseries_errors,weak_dense_errors=weak_dense_errors) for sim in _solutions]
  auxdata = Dict("dts" =>  dts)
  # Now Calculate Weak Errors
  additional_errors = Dict()
  for k in keys(solutions[1].weak_errors)
    additional_errors[k] = [sol.weak_errors[k] for sol in solutions]
  end
  ConvergenceSimulation(solutions,dts,auxdata=auxdata,additional_errors=additional_errors)
end

function test_convergence(dts::AbstractArray,prob::AbstractODEProblem,alg;
                          save_everystep=true,adaptive=false,kwargs...)
  N = length(dts)
  solutions = [solve(prob,alg;dt=dts[i],save_everystep=save_everystep,adaptive=adaptive,kwargs...) for i=1:N]
  auxdata = Dict(:dts =>  dts)
  ConvergenceSimulation(solutions,dts,auxdata=auxdata)
end

function analyticless_test_convergence(dts::AbstractArray,prob::AbstractODEProblem,
                                       alg,appxsol_setup;
                                       save_everystep=true,adaptive=false,kwargs...)
  true_sol = solve(prob,appxsol_setup[:alg];appxsol_setup...);
  N = length(dts)
  _solutions = [solve(prob,alg;dt=dts[i],save_everystep=save_everystep,adaptive=adaptive,kwargs...) for i=1:N]
  solutions = [appxtrue(sol,true_sol) for sol in _solutions]
  auxdata = Dict(:dts =>  dts)
  ConvergenceSimulation(solutions,dts,auxdata=auxdata)
end

function test_convergence(probs,convergence_axis,alg;kwargs...)
  ConvergenceSimulation([solve(prob,alg;kwargs...) for prob in probs],convergence_axis)
end

function test_convergence(c::ConvergenceSetup,alg::DEAlgorithm;kwargs...)
  test_convergence(c.probs,c.convergence_axis,alg;kwargs...)
end

function calc𝒪estimates(error::Pair)
  key = error.first
  error =error.second
  if ndims(error)>1 error=mean(error,1) end
  S = Vector{eltype(error)}(undef, length(error)-1)
  for i=1:length(error)-1
    S[i] = log2(error[i+1]/error[i])
  end
  return(Pair(key,abs.(mean(S,dims=1))))
end

"""
length(simres::ConvergenceSimulation)

Returns the number of simultations in the Convergence Simulation
"""
Base.length(sim::ConvergenceSimulation) = sim.N
Base.getindex(sim::ConvergenceSimulation,i::Int) = sim.solutions[i]
Base.getindex(sim::ConvergenceSimulation,i::Int,I::Int...) = sim.solutions[i][I]
Base.lastindex(sim::ConvergenceSimulation) = lastindex(sim.solutions)
Base.firstindex(sim::ConvergenceSimulation) = firstindex(sim.solutions)
Base.length(sim::ConvergenceSetup) = sim.probs


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 78e85936-bd55-11eb-23e1-67bc4a4b7901
using Pluto, PlutoUI, ImageFiltering, CairoMakie

# ╔═╡ 96d2bcc8-70c4-436d-888e-3e63a93cf8be
begin
	function generate_map(T::DataType, map_size, flower_width)
		n, m = map_size
		nan = convert(T, NaN)
		s1 = n ÷ 2 - flower_width ÷ 2
		s2 = s1 - flower_width % 2
		crop = [ones(T, s1, m); fill(nan, flower_width, m); ones(T, s2, m)]
		return crop
	end
	function generate_map(T::DataType, msize, fwidth, splits)
		h, w = msize
		nan = convert(T, NaN)
		s1 = (fwidth ÷ splits)^2
		s2 = (h ÷ splits) * (fwidth ÷ splits)
		cases = Iterators.product(1:h, 1:w) |> collect
		p = sortperm(abs.([(m * s1 + n * s2 - fwidth * h) * abs(m - n) for (m, n) in cases])[:])
		cases[p[1:5]]
	end
	generate_map(args...) = generate_map(Float32, args...)
end

# ╔═╡ c7dfd98d-49fa-4d73-bf99-f2cc33878e67
flower_width = @bind flower_width PlutoUI.Slider(10:20, default=15, show_value=true)

# ╔═╡ 5baa737f-005e-4733-8eaf-fdd2643b390e
generate_map((100, 100), flower_width, 3)

# ╔═╡ 6a34a0bf-b66d-4fbe-bb79-ef7b61a03668
begin
	crop = generate_map((100, 100), flower_width);
	count(isnan, crop), prod(size(crop))
end

# ╔═╡ cf84c816-b9d9-40b5-8af6-73e8e265a243
heatmap_kwargs = (
	nan_color=:yellow,
	colorrange = (0, 1),
	colormap = RGBAf0.(0, 1, 0, 0:0.01:1)
	)

# ╔═╡ 8fcf971f-08ff-46c4-9e33-308d6233c0df
heatmap(crop; heatmap_kwargs...)

# ╔═╡ 4e92d1ae-8e86-477d-8de6-fe358cba0721
pr_killed = imfilter(isnan.(crop) .* 0.2, Kernel.gaussian(3));

# ╔═╡ 4b8a5719-2832-48f2-9dbd-e1009bccdd61
heatmap(pr_killed)

# ╔═╡ Cell order:
# ╠═78e85936-bd55-11eb-23e1-67bc4a4b7901
# ╠═96d2bcc8-70c4-436d-888e-3e63a93cf8be
# ╠═c7dfd98d-49fa-4d73-bf99-f2cc33878e67
# ╠═5baa737f-005e-4733-8eaf-fdd2643b390e
# ╠═6a34a0bf-b66d-4fbe-bb79-ef7b61a03668
# ╠═cf84c816-b9d9-40b5-8af6-73e8e265a243
# ╠═8fcf971f-08ff-46c4-9e33-308d6233c0df
# ╠═4e92d1ae-8e86-477d-8de6-fe358cba0721
# ╠═4b8a5719-2832-48f2-9dbd-e1009bccdd61

max_stars_repo_name: ndgnuh/RiceBPH

max_issues_repo_name: ndgnuh/RiceBPH

max_forks_repo_name: ndgnuh/RiceBPH


text:
struct GridIndex{T, Index <: AbstractVector{Int}} <: AbstractIndex
    x0::T
    y0::T
    spacing::T
    n_x::Int
    n_y::Int
    start_indices::Vector{Int}
    index::Index
end

Base.summary(i::GridIndex) = "GridIndex ($(i.n_x)×$(i.n_y) cells of spacing $(i.spacing))"

Base.@propagate_inbounds Base.getindex(g::GridIndex, i, j) = g[i+g.n_x*(j-1)]
Base.@propagate_inbounds Base.getindex(g::GridIndex, i) = (g.start_indices[i]):(g.start_indices[i+1]-1)

function AcceleratedArrays.accelerate(points::AbstractVector{<:StaticVector{3, T}}, ::Type{GridIndex}; spacing = one(T)) where {T}
    grid = GridIndex(points; spacing = spacing)
    return AcceleratedArray(points, grid)
end

function AcceleratedArrays.accelerate!(points::AbstractVector{<:StaticVector{3, T}}, ::Type{GridIndex}; spacing = one(T), perm = Ref{Vector{Int}}()) where {T}
    grid = GridIndex(points; spacing = spacing)
    perm[] = grid.index
    permute!(points, grid.index)
    index = keys(points)
    return AcceleratedArray(points, GridIndex(grid.x0, grid.y0, grid.spacing, grid.n_x, grid.n_y, grid.start_indices, index))
end

"""
    GridIndex(points; spacing = 1)

An acceleration index for 3D points using a 2D grid. Inside each grid cell, the points are
ordered in the third dimension. Suitable for "two-and-a-half dimensional data" such as point
clouds on the surface of the Earth, and where the query radius is close to the grid
`spacing`.

Note that `GridIndex`s are typically created by the `accelerate` and `accelerate!`
functions (which also accept the `spacing` keyword argument. In the latter case, the points
are re-ordered to increase query speed.
"""
function GridIndex(points::AbstractVector{<:StaticVector{3, T}}; spacing = one(T)) where {T}
    grid_spacing = convert(T, spacing)
    n_points = length(points)
    if n_points < 1
        return GridIndex(zero(T), zero(T), grid_spacing, 1, 1, [1, 1], Int[])
    end
    if length(points[1]) != 3
        error("Points must be of length 3")
    end

    # Initialize some grid parameters
    (xmin, xmax, ymin, ymax) = bounds(points)
    xmin = prevfloat(xmin)
    ymin = prevfloat(ymin)

    n_x = Int(cld(xmax - xmin, grid_spacing))
    n_y = Int(cld(ymax - ymin, grid_spacing))

    # Pass over the points to construct the grid
    cells = Matrix{Vector{Int}}(undef, n_x, n_y)
    for i = 1:(n_x*n_y)
        @inbounds cells[i] = Vector{Int}()
    end

    @inbounds for i ∈ 1:n_points
        p = points[i]
        cell_x = Int(cld(p[1] - xmin, grid_spacing))
        cell_y = Int(cld(p[2] - ymin, grid_spacing))
        push!(cells[cell_x, cell_y], i)
    end

    # Create the permutation and sort the cells
    permutation = Vector{Int}(undef, n_points)
    start_indices = Vector{Int}(undef, n_x::Int*n_y::Int+1)
    k = 1
    for iy ∈ 1:n_y
        for ix ∈ 1:n_x
            @inbounds start_indices[ix + n_x*(iy-1)] = k
            @inbounds cell = cells[ix,iy]
            sort!(cell, Base.Sort.QuickSortAlg(), Base.Order.By(i -> (Base.@_inline_meta; @inbounds return points[i][3])))
            k_end = k + length(cell) - 1
            @inbounds permutation[k:k_end] = cell
            k = k_end + 1
        end
    end
    @inbounds start_indices[end] = n_points+1

    grid = GridIndex(xmin, ymin, grid_spacing, n_x, n_y, start_indices, permutation)
end

function bounds(points::AbstractVector)
    PointType = eltype(points)
    T = eltype(PointType)
    xmin = typemax(T)
    xmax = typemin(T)
    ymin = typemax(T)
    ymax = typemin(T)

    @inbounds for p ∈ points
        @inbounds x = p[1]
        if x < xmin
            xmin = x
        end
        if x > xmax
            xmax = x
        end

        @inbounds y = p[2]
        if y < ymin
            ymin = y
        end
        if y > ymax
            ymax = y
        end
    end

    return (xmin, xmax, ymin, ymax)
end

function cells_bbox(grid::GridIndex, region::BoundingBox)
    # Find the rough rectangle of intersecting cells
    # Axis-aligned bounding box:
    xmin = region.xmin
    xmax = region.xmax
    ymin = region.ymin
    ymax = region.ymax

    # Make a square grid
    gx_min = Int(cld(xmin - grid.x0, grid.spacing))
    if gx_min < 1
        gx_min = 1
    end
    gx_max = Int(cld(xmax - grid.x0, grid.spacing))
    if gx_max > grid.n_x
        gx_max = grid.n_x
    end

    gy_min = Int(cld(ymin - grid.y0, grid.spacing))
    if gy_min < 1
        gy_min = 1
    end
    gy_max = Int(cld(ymax - grid.y0, grid.spacing))
    if gy_max > grid.n_y
        gy_max = grid.n_y
    end

    return (gx_min, gx_max, gy_min, gy_max)
end

# Accelerations

function Base.findall(pred::Base.Fix2{typeof(in), <:AbstractRegion}, points::AcceleratedArray{<:Any, <:Any, <:Any, <:GridIndex})
    indices = Vector{Int}()
    findall!(indices, pred, points)
    return indices
end

function findall!(indices::Vector{<:Integer}, pred::Base.Fix2{typeof(in), <:AbstractRegion}, points::AcceleratedArray{<:Any, <:Any, <:Any, <:GridIndex})
    # First find the relevant grid cells
    grid = points.index
    bbox = boundingbox(pred.x)
    (gx_min, gx_max, gy_min, gy_max) = cells_bbox(grid, bbox)

    # Now we'll iterate over these cells and fill the indices
    @inbounds for gy ∈ gy_min:gy_max
        for gx ∈ gx_min:gx_max
            for i ∈ grid[gx, gy]
                # Could use min/max height information here, as in old code
                j = grid.index[i]
                p = points[j]

                if p[3] > bbox.zmax
                    break
                end

                if pred(p)
                    push!(indices, j)
                end
            end
        end
    end

    return indices
end

function Base.count(pred::Base.Fix2{typeof(in), <:AbstractRegion}, points::AcceleratedArray{<:Any, <:Any, <:Any, <:GridIndex})
    # First find the relevant grid cells
    grid = points.index
    (gx_min, gx_max, gy_min, gy_max) = cells_bbox(grid, boundingbox(pred.x))

    out = 0

    # Now we'll iterate over these cells and fill the indices
    @inbounds for gy ∈ gy_min:gy_max
        for gx ∈ gx_min:gx_max
            for i ∈ grid[gx, gy]
                # Could use min/max height information here, as in commented code below
                j = grid.index[i]
                p = points[j]

                out += pred(p)
            end
        end
    end

    return out
end

function Base.filter(pred::Base.Fix2{typeof(in), <:AbstractRegion}, points::AcceleratedArray{<:Any, <:Any, <:Any, <:GridIndex})
    # First find the relevant grid cells
    grid = points.index
    (gx_min, gx_max, gy_min, gy_max) = cells_bbox(grid, boundingbox(pred.x))

    out = empty(points.parent)

    # Now we'll iterate over these cells and fill the indices
    @inbounds for gy ∈ gy_min:gy_max
        for gx ∈ gx_min:gx_max
            for i ∈ grid[gx, gy]
                # Could use min/max height information here, as in commented code below
                j = grid.index[i]
                p = points[j]

                if pred(p)
                    push!(out, p)
                end
            end
        end
    end

    return out
end


text:
export groupcode, read_value
export groupcode, groupmatch
export  DXFObject, DXFGroup, RawGroup, StringGroup, IntegerGroup,
    FloatGroup, PointX, PointY, PointZ

abstract type DXFObject end
abstract type DXFGroup <: DXFObject end

DXFGroupCode = Int32

group_code_registry = SortedDict{DXFGroupCode, Type{<:DXFGroup}}()

function show_group_code_registry(io::IO)
    for (key, value) in group_code_registry
        @printf(io, "%4d\t%s\n", key,
                join(pedigree(value), " <: "))
    end
end

function groupmatch(a::DXFGroup, b::DXFGroup)
    groupcode(a) == groupcode(b) && a.value == b.value
end

function defDXFGroup_(group_supertype::Type, code::Integer, name=nothing)
    defDXFGroup_(group_supertype, DXFGroupCode(code), name)
end

function defDXFGroup_(group_supertype::Type, code::DXFGroupCode, name=nothing)
    @assert group_supertype <: DXFGroup
    if name == nothing
        name = Symbol("Group_$code")
    end
    definitions = quote
        struct $name <: $group_supertype
            line::Int
            value::$(valuetype(group_supertype))
        end
        function $name(value::$(valuetype(group_supertype)))
            $name(-1, value)
        end
        if !haskey(group_code_registry, $code)
            group_code_registry[$code] = $name
        end
        function groupcode(::$name)::DXFGroupCode
            return $code
        end
        function groupcode(::Type{$name})::DXFGroupCode
            return $code
        end
        export $name
    end
    return name, definitions
end

function defDXFGroup(group_supertype::Type, code::Integer, name=nothing)
    name, defs = defDXFGroup_(group_supertype, code, name)
    eval(defs)
    eval(name)
end

abstract type RawGroup <: DXFGroup end

function valuetype(::Type{T}) where {T <: RawGroup}
    String
end

function Base.one(t::Type{<:DXFGroup})
    t(one(valuetype(t)))
end

function read_value(group_type::Type{T}, in) where {T <: RawGroup}
    group_type(in.line, readline(in))
end

function lookup_group_code(code::DXFGroupCode)
    if haskey(group_code_registry, code)
        return group_code_registry[code]
    end
    return defDXFGroup(RawGroup, code)
end

abstract type StringGroup <: DXFGroup end

function valuetype(::Type{T}) where {T <: StringGroup}
    String
end

function read_value(group_type::Type{T}, in) where {T <: StringGroup}
    group_type(in.line, strip(readline(in)))
end

# It makes it easier to write the parser if we subtype group 0 for
# each different kind of element type.

abstract type EntityType <: DXFGroup end

export EntityType

function valuetype(::Type{T}) where {T <: EntityType}
    String
end

function ensureEntityType(value)
    etname = Symbol("EntityType_$value")
    if !any([t.name.name == etname
             for t in allsubtypes(EntityType)])
        name, defs = defDXFGroup_(EntityType, 0, etname)
        eval(defs)
    end
    eval(etname)
end

function EntityType(line::Int, value::String)
    Base.invokelatest(ensureEntityType(value), line, value)
end

function EntityType(value::String)
    Base.invokelatest(ensureEntityType(value), value)
end

group_code_registry[0] = EntityType

groupcode(::Type{<:EntityType})::DXFGroupCode = 0
groupcode(::EntityType)::DXFGroupCode = 0

function read_value(group_type::Type{T}, in) where {T <: EntityType}
    value = strip(readline(in))
    Base.invokelatest(ensureEntityType(value),
                      in.line, value)
end

# We need some of these EntityTypes at compile time to compile
# parser.jl:

ensureEntityType("SECTION")
ensureEntityType("ENDSEC")
ensureEntityType("EOF")
ensureEntityType("BLOCK")
ensureEntityType("ENDBLK")

#=
ensureEntityType should have defined EntityType_CLASS
ensureEntityType should have defined EntityType_TABLE
ensureEntityType should have defined EntityType_VPORT
ensureEntityType should have defined EntityType_ENDTAB
ensureEntityType should have defined EntityType_LTYPE
ensureEntityType should have defined EntityType_LAYER
ensureEntityType should have defined EntityType_STYLE
ensureEntityType should have defined EntityType_APPID
ensureEntityType should have defined EntityType_DIMSTYLE
ensureEntityType should have defined EntityType_BLOCK_RECORD
ensureEntityType should have defined EntityType_LINE
ensureEntityType should have defined EntityType_DICTIONARY
ensureEntityType should have defined EntityType_ACDBDICTIONARYWDFLT
ensureEntityType should have defined EntityType_XRECORD
ensureEntityType should have defined EntityType_SORTENTSTABLE
ensureEntityType should have defined EntityType_LAYOUT
ensureEntityType should have defined EntityType_MATERIAL
ensureEntityType should have defined EntityType_MLEADERSTYLE
ensureEntityType should have defined EntityType_MLINESTYLE
ensureEntityType should have defined EntityType_ACDBPLACEHOLDER
ensureEntityType should have defined EntityType_SCALE
ensureEntityType should have defined EntityType_TABLESTYLE
ensureEntityType should have defined EntityType_VISUALSTYLE
ensureEntityType should have defined EntityType_DICTIONARYVAR
=#

defDXFGroup(StringGroup, 1, :PrimaryText)
defDXFGroup(StringGroup, 2, :Name)
defDXFGroup(StringGroup, 3, :EntryName)
defDXFGroup(StringGroup, 4)
defDXFGroup(StringGroup, 5, :EntityHandle)
defDXFGroup(StringGroup, 6)
defDXFGroup(StringGroup, 7, :TextStyleName)
defDXFGroup(StringGroup, 8, :LayerName)
defDXFGroup(StringGroup, 9, :HeaderVariableName)

DXFIntegerType = Int64

abstract type IntegerGroup <: DXFGroup end

function valuetype(::Type{T}) where {T <: IntegerGroup}
    DXFIntegerType
end

function read_value(group_type::Type{T}, in) where {T <: IntegerGroup}
    group_type(in.line, parse(Int64, strip(readline(in))))
end

for code in 70:78
    defDXFGroup(IntegerGroup, code)
end

DXFFloatType = Float64

abstract type FloatGroup <: DXFGroup end

function valuetype(::Type{T}) where {T <: FloatGroup}
    DXFFloatType
end

function read_value(group_type::Type{T}, in) where {T <: FloatGroup}
    group_type(in.line, parse(DXFFloatType, strip(readline(in))))
end

abstract type PointX <: FloatGroup end
abstract type PointY <: FloatGroup end
abstract type PointZ <: FloatGroup end

export PointX, PointY, PointZ

for code in 10:18
    defDXFGroup(PointX, code,      Symbol("Point$(code)_X"))
    defDXFGroup(PointY, code + 10, Symbol("Point$(code + 10)_Y"))
    defDXFGroup(PointZ, code + 20, Symbol("Point$(code + 20)_Z"))
end

defDXFGroup(PointX, 210, :ExtrusionDirectionX)
defDXFGroup(PointY, 220, :ExtrusionDirectionY)
defDXFGroup(PointZ, 230, :ExtrusionDirectionZ)

for code in 40:47
    defDXFGroup(FloatGroup, code)
end

defDXFGroup(FloatGroup, 48, :LinetypeScale)

for code in 49:59
    defDXFGroup(FloatGroup, code)
end

defDXFGroup(StringGroup, 100, :SubclassMarker)

defDXFGroup(StringGroup, 102, :ADGroupStartEnd)

for code in 110:149
    defDXFGroup(FloatGroup, code)
end

for code in 160:179
    defDXFGroup(IntegerGroup, code)
end

for code in 210:239
    defDXFGroup(FloatGroup, code)
end

for code in 270:289
    defDXFGroup(IntegerGroup, code)
end

abstract type HexIdString <: StringGroup end

for code in 330:369
    defDXFGroup(HexIdString, code)
end

for code in 370:389
    defDXFGroup(IntegerGroup, code)
end

abstract type HexHandleString <: StringGroup end

for code in 390:399
    defDXFGroup(HexHandleString, code)
end

defDXFGroup(StringGroup, 999, :Comment)


text:
using Oceananigans.AbstractOperations: GridMetricOperation

const c = Center()
const f = Face()
const IBG = ImmersedBoundaryGrid

# Grid metrics for ImmersedBoundaryGrid
#
# All grid metrics are defined here.
#
# For non "full-cell" immersed boundaries, grid metric functions
# must be extended for the specific immersed boundary grid in question.
#
for LX in (:ᶜ, :ᶠ), LY in (:ᶜ, :ᶠ), LZ in (:ᶜ, :ᶠ)
    for dir in (:x, :y, :z), operator in (:Δ, :A)

        metric = Symbol(operator, dir, LX, LY, LZ)
        @eval begin
            import Oceananigans.Operators: $metric
            @inline $metric(i, j, k, ibg::IBG) = $metric(i, j, k, ibg.underlying_grid)
        end
    end

    volume = Symbol(:V, LX, LY, LZ)
    @eval begin
        import Oceananigans.Operators: $volume
        @inline $volume(i, j, k, ibg::IBG) = $volume(i, j, k, ibg.underlying_grid)
    end
end

@inline Δzᵃᵃᶜ(i, j, k, ibg::IBG) = Δzᵃᵃᶜ(i, j, k, ibg.underlying_grid)
@inline Δzᵃᵃᶠ(i, j, k, ibg::IBG) = Δzᵃᵃᶠ(i, j, k, ibg.underlying_grid)


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

# 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 Hd  = Chemostat_Heerden2013.HeerdenData;
const ECC = Chemostat_Heerden2013.EColiCore

# run add "https://github.com/josePereiro/Chemostat" in the
# julia Pkg REPL for installing the package
import Chemostat
import Chemostat.LP.MathProgBase
const ChU = Chemostat.Utils
const ChSS = Chemostat.SteadyState
const ChLP = Chemostat.LP
const ChEP = Chemostat.MaxEntEP
const ChSU = Chemostat.SimulationUtils

using Plots

## -------------------------------------------------------------------
# save results
DATA = ChU.load_data(ECC.MAXENT_RES_FILE)
SDATA = sort(collect(DATA); by = first);

## -------------------------------------------------------------------
# flux vs beta
let
    p = plot(title = "EColiCore", xlabel = "beta", ylabel = "biom")
    for (exp, D) in SDATA
        model = D["model"]
        objidx = ChU.rxnindex(model, ECC.BIOMASS_IDER)
        epouts = D["epouts"]
        exp_beta = D["exp_beta"]
        exp_xi = D["exp_xi"]
        scatter!(p, [exp_beta], [Hd.val("D", exp)], ms = 12, color = :white, label = "")

        betas = collect(keys(epouts)) |> sort
        bioms = [ChU.av(model, epouts[beta], objidx) for beta in betas]
        scatter!(p, betas, bioms, label = "", color = :black, alpha = 0.2)

    end
    p
end

## -------------------------------------------------------------------
# beta_exp corr
let
    p = plot(title = "EColiCore", xlabel = "exp biom", ylabel = "model biom")
    for (exp, D) in SDATA

        model = D["model"]
        objidx = ChU.rxnindex(model, ECC.BIOMASS_IDER)
        exp_beta = D["exp_beta"]
        epout = D["epouts"][exp_beta]

        scatter!(p, [ChU.av(model, epout, objidx)], [Hd.val("D", exp)],
            color = :white, label = "")

    end
    p
end

## -------------------------------------------------------------------
# total correlations
let
    ep_p = plot(title = "EColiCore (EP)",
        xlabel = "model flux", ylabel = "exp flux")
    fba_p = plot(title = "EColiCore (FBA)",
        xlabel = "model flux", ylabel = "exp flux")

    m, M = Inf, -Inf
    for (exp, D) in SDATA

        model = D["model"]
        objidx = ChU.rxnindex(model, ECC.BIOMASS_IDER)
        exp_beta = D["exp_beta"]
        epout = D["epouts"][exp_beta]
        fbaout = D["fbaout"]

        ep_xs, ep_errs = [], []
        fba_xs = []
        ys = []

        # Biomass
        push!(fba_xs, ChU.av(model, fbaout, objidx))
        push!(ep_xs, ChU.av(model, epout, objidx))
        push!(ys, Hd.val("D", exp))

        # mets
        for Hd_met in Hd.msd_mets
            try
                model_met = ECC.Hd_mets_map[Hd_met]
                model_exch = ECC.exch_met_map[model_met]

                push!(fba_xs, ChU.av(model, fbaout, model_exch))
                push!(ep_xs, ChU.av(model, epout, model_exch))
                push!(ep_errs, sqrt(ChU.va(model, epout, model_exch)))
                push!(ys, Hd.val("u$Hd_met", exp))
            catch
                @warn string(Hd_met, " fails")
            end
        end

        scatter!(ep_p, ep_xs, ys; xerr = ep_errs, color = :black, label = "")
        scatter!(fba_p, fba_xs, ys, color = :black, label = "")
        m = min(m, minimum(ep_xs), minimum(fba_xs), minimum(ys))
        M = max(M, maximum(ep_xs), maximum(fba_xs), maximum(ys))

    end

    plot!(ep_p, [m,M], [m,M]; ls = :dash, color = :black, label = "")
    plot!(fba_p, [m,M], [m,M]; ls = :dash, color = :black, label = "")
    plot([ep_p, fba_p]..., size = [700, 400])
end

## -------------------------------------------------------------------
# total correlations (conc)
let
    ep_p = plot(title = "EColiCore (EP)",
        xlabel = "model conc", ylabel = "exp conc")
    fba_p = plot(title = "EColiCore (FBA)",
        xlabel = "model conc", ylabel = "exp conc")

    m, M = Inf, -Inf
    for (exp, D) in SDATA

        model = D["model"]
        objidx = ChU.rxnindex(model, ECC.BIOMASS_IDER)
        exp_beta = D["exp_beta"]
        exp_xi = Hd.val("xi", exp)
        epout = D["epouts"][exp_beta]
        fbaout = D["fbaout"]

        ep_xs, ep_errs = [], []
        fba_xs = []
        ys = []

        # mets
        for Hd_met in Hd.msd_mets
            try
                model_met = ECC.Hd_mets_map[Hd_met]
                model_exch = ECC.exch_met_map[model_met]

                # conc (s = c + u*xi)
                c = Hd.val("c$Hd_met", exp, 0.0)

                u = ChU.av(model, fbaout, model_exch)
                push!(fba_xs, max(c + u*exp_xi, 0.0))

                u = ChU.av(model, epout, model_exch)
                push!(ep_xs, max(c + u*exp_xi, 0.0))

                uerr = sqrt(ChU.va(model, epout, model_exch)) * exp_xi
                push!(ep_errs, uerr)

                push!(ys, Hd.val("s$Hd_met", exp))

            catch
                @warn string(Hd_met, " fails")
            end
        end

        @assert length(ep_xs) == length(fba_xs) == length(ys)
        scatter!(ep_p, ep_xs, ys; xerr = ep_errs, color = :black, label = "")
        scatter!(fba_p, fba_xs, ys, color = :black, label = "")
        m = min(m, minimum(ep_xs), minimum(fba_xs), minimum(ys))
        M = max(M, maximum(ep_xs), maximum(fba_xs), maximum(ys))

    end

    # plot!(ep_p, [m,M], [m,M]; ls = :dash, color = :black, label = "")
    # plot!(fba_p, [m,M], [m,M]; ls = :dash, color = :black, label = "")
    plot([ep_p, fba_p]..., size = [700, 400])
end


text:
using DynamicGrids, StaticArrays, Test, FileIO, Colors, FixedPointNumbers
using DynamicGrids: SimData, NoDisplayImageOutput

@testset "CellRule that multiples a StaticArray" begin
    rule = Cell{:grid1}() do data, state, I
         2state
    end
    init = (grid1 = fill(SA[1.0, 2.0], 5, 5),)
    output = ArrayOutput(init; tspan=1:3)
    sim!(output, rule)
    @test output[2][:grid1] == reshape([ # Have to use reshape to construct this
        SA[2.0, 4.0], SA[2.0, 4.0], SA[2.0, 4.0], SA[2.0, 4.0], SA[2.0, 4.0],
        SA[2.0, 4.0], SA[2.0, 4.0], SA[2.0, 4.0], SA[2.0, 4.0], SA[2.0, 4.0],
        SA[2.0, 4.0], SA[2.0, 4.0], SA[2.0, 4.0], SA[2.0, 4.0], SA[2.0, 4.0],
        SA[2.0, 4.0], SA[2.0, 4.0], SA[2.0, 4.0], SA[2.0, 4.0], SA[2.0, 4.0],
        SA[2.0, 4.0], SA[2.0, 4.0], SA[2.0, 4.0], SA[2.0, 4.0], SA[2.0, 4.0]
    ], (5, 5))
    @test output[3][:grid1] == reshape([ # Have to use reshape to construct this
        SA[4.0, 8.0], SA[4.0, 8.0], SA[4.0, 8.0], SA[4.0, 8.0], SA[4.0, 8.0],
        SA[4.0, 8.0], SA[4.0, 8.0], SA[4.0, 8.0], SA[4.0, 8.0], SA[4.0, 8.0],
        SA[4.0, 8.0], SA[4.0, 8.0], SA[4.0, 8.0], SA[4.0, 8.0], SA[4.0, 8.0],
        SA[4.0, 8.0], SA[4.0, 8.0], SA[4.0, 8.0], SA[4.0, 8.0], SA[4.0, 8.0],
        SA[4.0, 8.0], SA[4.0, 8.0], SA[4.0, 8.0], SA[4.0, 8.0], SA[4.0, 8.0]
    ], (5, 5))
end

@testset "Neighborhood Rule that sums a Neighborhood of StaticArrays" begin
    rule = Neighbors{Tuple{:grid1,:grid2},:grid1}(Moore(1)) do data, hood, (s1, s2), I
        sum(hood) .+ s2
    end
    init = (
        grid1 = fill(SA[1.0, 2.0], 5, 5),
        grid2 = fill(0.5, 5, 5),
    )
    output = ArrayOutput(init; tspan=1:2)
    sim!(output, rule)
    @test output[2][:grid1] == reshape([ # Have to use reshape to construct this
        SA[3.5, 6.5],  SA[5.5, 10.5], SA[5.5, 10.5], SA[5.5, 10.5], SA[3.5, 6.5],
        SA[5.5, 10.5], SA[8.5, 16.5], SA[8.5, 16.5], SA[8.5, 16.5], SA[5.5, 10.5],
        SA[5.5, 10.5], SA[8.5, 16.5], SA[8.5, 16.5], SA[8.5, 16.5], SA[5.5, 10.5],
        SA[5.5, 10.5], SA[8.5, 16.5], SA[8.5, 16.5], SA[8.5, 16.5], SA[5.5, 10.5],
        SA[3.5, 6.5],  SA[5.5, 10.5], SA[5.5, 10.5], SA[5.5, 10.5], SA[3.5, 6.5]
    ], (5, 5))
end

@testset "SetCell randomly updates a StaticArray" begin
    rule = SetCell{:grid1}() do data, state, I
        if I == (2, 2) || I == (1, 3)
            data[:grid1][I...] = SA[99.0, 100.0]
        end
    end
    init = (grid1 = fill(SA[0.0, 0.0], 3, 3),)
    output = ArrayOutput(init; tspan=1:2)
    sim!(output, rule)
    @test output[2][:grid1] == reshape([ # Have to use reshape to construct this
        SA[0.0, 0.0], SA[0.0, 0.0], SA[0.0, 0.0],
        SA[0.0, 0.0], SA[99.0, 100.0], SA[0.0, 0.0],
        SA[99.0, 100.0], SA[0.0, 0.0], SA[0.0, 0.0]
    ], (3, 3))
end

struct TestStruct{A,B}
    a::A
    b::B
end
const TS = TestStruct

Base.:*(ts::TestStruct, x::Number) = TestStruct(ts.a * x, ts.b * x)
Base.:*(x::Number, ts::TestStruct) = TestStruct(x * ts.a, x * ts.b)
Base.:/(ts::TestStruct, x::Number) = TestStruct(ts.a / x, ts.b / x)
Base.:+(ts1::TestStruct, ts2::TestStruct) = TestStruct(ts1.a + ts2.a, ts1.b + ts2.b)
Base.:-(ts1::TestStruct, ts2::TestStruct) = TestStruct(ts1.a - ts2.a, ts1.b - ts2.b)
Base.:+(ts::TestStruct, x::Number) = TestStruct(ts.a + x, ts.b + x)
Base.:-(ts::TestStruct, x::Number) = TestStruct(ts.a - x, ts.b - x)

Base.isless(a::TestStruct, b::TestStruct) = isless(a.a, b.a)
Base.zero(::Type{<:TestStruct{T1,T2}}) where {T1,T2} = TestStruct(zero(T1), zero(T2))
Base.oneunit(::Type{<:TestStruct{T1,T2}}) where {T1,T2} = TestStruct(oneunit(T1), oneunit(T2))

DynamicGrids.to_rgb(scheme::ObjectScheme, obj::TestStruct) = ARGB32(obj.a)
DynamicGrids.to_rgb(scheme, obj::TestStruct) = get(scheme, obj.a)

@testset "CellRule that multiples a struct" begin
    rule = Cell{:grid1,:grid1}() do data, state, I
         2state
    end
    init = (grid1 = fill(TS(1.0, 2.0), 5, 5),)
    output = ArrayOutput(init; tspan=1:3)
    sim!(output, rule)
    @test output[2][:grid1] == reshape([ # Have to use reshape to construct this
        TS(2.0, 4.0), TS(2.0, 4.0), TS(2.0, 4.0), TS(2.0, 4.0), TS(2.0, 4.0),
        TS(2.0, 4.0), TS(2.0, 4.0), TS(2.0, 4.0), TS(2.0, 4.0), TS(2.0, 4.0),
        TS(2.0, 4.0), TS(2.0, 4.0), TS(2.0, 4.0), TS(2.0, 4.0), TS(2.0, 4.0),
        TS(2.0, 4.0), TS(2.0, 4.0), TS(2.0, 4.0), TS(2.0, 4.0), TS(2.0, 4.0),
        TS(2.0, 4.0), TS(2.0, 4.0), TS(2.0, 4.0), TS(2.0, 4.0), TS(2.0, 4.0)
    ], (5, 5))
    @test output[3][:grid1] == reshape([ # Have to use reshape to construct this
        TS(4.0, 8.0), TS(4.0, 8.0), TS(4.0, 8.0), TS(4.0, 8.0), TS(4.0, 8.0),
        TS(4.0, 8.0), TS(4.0, 8.0), TS(4.0, 8.0), TS(4.0, 8.0), TS(4.0, 8.0),
        TS(4.0, 8.0), TS(4.0, 8.0), TS(4.0, 8.0), TS(4.0, 8.0), TS(4.0, 8.0),
        TS(4.0, 8.0), TS(4.0, 8.0), TS(4.0, 8.0), TS(4.0, 8.0), TS(4.0, 8.0),
        TS(4.0, 8.0), TS(4.0, 8.0), TS(4.0, 8.0), TS(4.0, 8.0), TS(4.0, 8.0)
    ], (5, 5))
end

@testset "Neighborhood Rule that sums a Neighborhood of stucts" begin
    rule = Neighbors{Tuple{:grid1,:grid2},:grid1}(Moore(1)) do data, hood, (s1, s2), I
        sum(hood) * s2
    end
    init = (
        grid1 = fill(TS(1.0, 2.0), 5, 5),
        grid2 = fill(0.5, 5, 5),
    )
    output = ArrayOutput(init; tspan=1:2)
    sim!(output, rule)
    @test output[2][:grid1] == reshape([ # Have to use reshape to construct this
        TS(1.5, 3.0), TS(2.5, 5.0), TS(2.5, 5.0), TS(2.5, 5.0), TS(1.5, 3.0),
        TS(2.5, 5.0), TS(4.0, 8.0), TS(4.0, 8.0), TS(4.0, 8.0), TS(2.5, 5.0),
        TS(2.5, 5.0), TS(4.0, 8.0), TS(4.0, 8.0), TS(4.0, 8.0), TS(2.5, 5.0),
        TS(2.5, 5.0), TS(4.0, 8.0), TS(4.0, 8.0), TS(4.0, 8.0), TS(2.5, 5.0),
        TS(1.5, 3.0), TS(2.5, 5.0), TS(2.5, 5.0), TS(2.5, 5.0), TS(1.5, 3.0),
    ], (5, 5))
end

@testset "SetCell rule randomly updates a struct" begin
    rule = SetCell{:grid1,:grid1}() do data, state, I
        if I == (2, 2) || I == (1, 3)
            data[:grid1][I...] = TS(99.0, 100.0)
        end
    end
    init = (grid1 = fill(TS(0.0, 0.0), 3, 3),)
    output = ArrayOutput(init; tspan=1:2)
    sim!(output, rule)
    @test output[2][:grid1] == reshape([ # Have to use reshape to construct this
        TS(0.0, 0.0), TS(0.0, 0.0), TS(0.0, 0.0),
        TS(0.0, 0.0), TS(99.0, 100.0), TS(0.0, 0.0),
        TS(99.0, 100.0), TS(0.0, 0.0), TS(0.0, 0.0)
    ], (3, 3))
end

@testset "object grid can generate an image" begin
    @testset "normalise" begin
        @test DynamicGrids.to_rgb(ObjectScheme(), TestStruct(99.0, 1.0) / 99) == ARGB32(1.0)
        @test DynamicGrids.to_rgb(ObjectScheme(), TestStruct(00.0, 0.0) / 99) == ARGB32(0.0)
        @test DynamicGrids.to_rgb(ObjectScheme(), DynamicGrids.normalise(TestStruct(99.0, 1.0), nothing, 99)) == ARGB32(1.0)
    end
    rule = SetCell{:grid1,:grid1}() do data, state, I
        if I == (2, 2) || I == (1, 3)
            data[:grid1][I...] = TS(99.0, 100.0)
        end
    end
    init = (grid1=fill(TS(0.0, 0.0), 3, 3),)
    # These should have the same answer
    output1 = GifOutput(init;
        filename="objectgrid.gif", store=true, tspan=1:2, maxval=[99.0], text=nothing
    )
    output2 = GifOutput(init;
        filename="objectgrid_greyscale.gif", scheme=Greyscale(), store=true, tspan=1:2,
        maxval=reshape([99.0], 1, 1), text=nothing
    )
    sim!(output1, rule)
    sim!(output2, rule)
    @test output1[2][:grid1] ==
       [TS(0.0, 0.0) TS(0.0, 0.0) TS(99.0, 100.0)
        TS(0.0, 0.0) TS(99.0, 100.0) TS(0.0, 0.0)
        TS(0.0, 0.0) TS(0.0, 0.0) TS(0.0, 0.0)]
    @test RGB.(output1.gif[:, :, 2]) ==
          RGB.(output2.gif[:, :, 2]) ==
          load("objectgrid.gif")[:, :, 2] ==
          load("objectgrid_greyscale.gif")[:, :, 2] ==
          map(xs -> RGB{N0f8}(xs...),
              [(0.298,0.298,0.298) (0.298,0.298,0.298) (1.0,1.0,1.0)
               (0.298,0.298,0.298) (1.0,1.0,1.0) (0.298,0.298,0.298)
               (0.298,0.298,0.298) (0.298,0.298,0.298) (0.298,0.298,0.298)]
          )
end

@testset "static arrays grid can generate an image" begin
    rule = Cell{:grid1}() do data, state, I
         2state
    end
    init = (grid1 = fill(SA[1.0, 2.0], 5, 5),)

    # We can index into the SArray or access
    # it with a function, defined using a Pair
    output = GifOutput(init;
        filename="sa.gif",
        tspan=1:3,
        store=true,
        layout=[:grid1=>1 :grid1=>x->x[2]],
        scheme=[Greyscale() Greyscale()],
        minval=[0.0 0.0], maxval=[10.0 10.0],
        text=nothing,
    )

    sim!(output, rule)
    a02 = ARGB32(0.2)
    a04 = ARGB32(0.4)
    @test output.gif[:, :, 2] ==
        [a02 a02 a02 a02 a02 a04 a04 a04 a04 a04
         a02 a02 a02 a02 a02 a04 a04 a04 a04 a04
         a02 a02 a02 a02 a02 a04 a04 a04 a04 a04
         a02 a02 a02 a02 a02 a04 a04 a04 a04 a04
         a02 a02 a02 a02 a02 a04 a04 a04 a04 a04]
end


text:
# Copyright (c) 2014 Aubrey Barnard.  This is free software.  See
# LICENSE for details.

# List with block-based storage that supports random access.  The
# benefits are that allocation and deallocation do not require copying
# and that at most a constant amount of memory is unused (the block
# size).  The costs are an extra level of indirection for access and
# potentially more allocations (although no copies) for a given data
# size.

module jltbl

import Base: length, size, getindex, setindex!, append!, start, next, done

# Constants

const _defaultBlockSize = 1024
const _defaultTopLevelCapacity = 16

# BlockedList
type BlockedList{T} <: AbstractVector{T}
    blockSize::Int
    blocks::Vector{Vector{T}} # The top-level array is pointer-based
    blocksLength::Int # Number of valid top-level elements
    length::Int # Total number of valid elements

    function BlockedList(blockSize::Int, blocks::Vector{Vector{T}})
        # Check for valid arguments
        blockSize > 0 || error("Block size not positive: $(blockSize)")
        # Construct
        return new(blockSize, blocks, 0, 0)
    end
end

# Constructors

BlockedList{T}(eltype::Type{T}) = BlockedList(eltype, _defaultBlockSize)
BlockedList{T}(eltype::Type{T}, blockSize::Int) = BlockedList{T}(blockSize, Array(Vector{T}, _defaultTopLevelCapacity))
BlockedList{T}(blockSize::Int, blocks::Vector{Vector{T}}) = BlockedList{T}(blockSize, blocks)

# Length, size, capacity

length(list::BlockedList) = list.length

size(list::BlockedList) = (list.length,)

capacity(list::BlockedList) = list.blocksLength * list.blockSize

# Element access

function getindex(list::BlockedList, index::Int)
    checkbounds(list.length, index)
    index, subindex = divrem(index - 1, list.blockSize)
    return list.blocks[index + 1][subindex + 1]
end

function setindex!{T}(list::BlockedList{T}, value::T, index::Int)
    checkbounds(list.length, index)
    return _setindex!(list, value, index)
end

function _setindex!{T}(list::BlockedList{T}, value::T, index::Int)
    index, subindex = divrem(index - 1, list.blockSize)
    list.blocks[index + 1][subindex + 1] = value
    return list
end

# Appending

function ensureCapacity!{T}(list::BlockedList{T}, newCapacity::Int)
    # How many blocks are needed for the new capacity?
    newBlocksLength = iceil(newCapacity / list.blockSize)
    blocksCapacity = length(list.blocks)
    # Allocate top level as needed
    newBlocksCapacity = blocksCapacity
    while newBlocksCapacity < newBlocksLength
        newBlocksCapacity *= 2
    end
    if newBlocksCapacity > blocksCapacity
        resize!(list.blocks, newBlocksCapacity)
    end
    # Allocate blocks as needed
    while list.blocksLength < newBlocksLength
        list.blocksLength += 1
        list.blocks[list.blocksLength] = Array(T, list.blockSize)
    end
    return list
end

# Append like in the standard library, where the number of new elements
# is known
function append!{T}(list::BlockedList{T}, items::AbstractVector{T})
    numberItems = length(items)
    ensureCapacity!(list, list.length + numberItems)
    for index in 1:numberItems
        _setindex!(list, items[index], list.length + index)
    end
    list.length += numberItems
    return list
end

# Append like in the standard library, where the number of new elements
# is not known
function append!{T}(list::BlockedList{T}, items) # TODO items::Iterable
    baseIndex = list.length
    cap = capacity(list)
    index = 1
    for (index, item) in enumerate(items) # TODO start index at baseIndex
        # Check the type of the item
        typeof(item) == T || error("Wrong element type: $(typeof(item)) is not $(T): $(item)")
        # Add capacity for this item if needed
        if baseIndex + index > cap
            ensureCapacity!(list, baseIndex + index)
            cap = capacity(list)
        end
        _setindex!(list, item, baseIndex + index)
    end
    list.length += index
    return list
end

# Append like in Go
function append!{T}(list::BlockedList{T}, items::T...)
    return append!(list, items)
end

# Iteration

start(list::BlockedList) = Int[1, 1]

function next(list::BlockedList, state)
    # Get the next item
    item = list.blocks[state[1]][state[2]]
    # Advance the indices
    state[2] += 1
    if state[2] > list.blockSize
        state[1] += 1
        state[2] = 1
    end
    return (item, state)
end

function done(list::BlockedList, state)
    lastLength = (list.length - 1) % list.blockSize + 1
    return (state[1] > list.blocksLength ||
        (state[1] == list.blocksLength && state[2] > lastLength))
end

# TODO implement the following functions:
# deleteat! (2 variants) (?)
# splice! (2 variants) (?)
# insert! (?)

end # module jltbl

max_stars_repo_name: afbarnard/jltbl

max_issues_repo_name: afbarnard/jltbl

max_forks_repo_name: afbarnard/jltbl


text:
using PkgSRA, Plots, LinearAlgebra, Statistics
using AxisArrays, DataFrames, DataStructures
using BSON: @load
pyplot()

include("paper_plotting.jl")
include("paper_settings.jl")

#####
##### Plot AND save using external function
#####

## Lotka Volterra
this_dat_name = "dat_library_with_noise_lv_"
plot_final = plot_library_noise(this_dat_name, "Lotka Volterra";
                        plot_opt=plot_opt)

this_fig_name = FIGURE_FOLDERNAME*"fig_library_with_noise_";
fname = this_fig_name * "lotkaVolterra.png";
savefig(plot_final, fname)

## FitzHugh Nagumo
this_dat_name = "dat_library_with_noise_fhn_"
plot_final = plot_library_noise(this_dat_name, "FitzHugh Nagumo";
                        plot_opt=plot_opt)

this_fig_name = FIGURE_FOLDERNAME*"fig_library_with_noise_";
fname = this_fig_name * "fitzHughNagumo.png";
savefig(plot_final, fname)


text:
# Let's begin, shall we?

# I'm going to begin with just 3 oscillators, because
# I haven't touched differential equations with Julia
# in quite a while. Hopefully the code will generalize.

#######################################################################
# First, of course, we're going to need our tools.

# To get these, you go into the REPL, type a
#
# ]
#
# at which point it'll change from "julia>" to "pkg>".
#
# Then you can just write in
#
# add Plots
# add DifferentialEquations
#
# or whatever other packages we end up using.

using Plots
using DifferentialEquations

# I'm being super verbose here so I don't forget anything.
#
# https://en.wikibooks.org/wiki/Introducing_Julia/Modules_and_packages#Modules_and_packages
# ^ source for my info.
########################################################################
# So first I'm just going to implement a system of equations, nothing
# fancy yet.
#
# http://docs.juliadiffeq.org/latest/tutorials/ode_example.html#Example-2:-Solving-Systems-of-Equations-1
#
# to see the example I'm working from.
#
# https://en.wikipedia.org/wiki/Kuramoto_model
#
# First diff eq is what I'm implementing here; might post LaTeX code
# later if I feel like it, and figure out some literate programming
# approach.

# Let's implement the trivial case, of N = 1 single oscillator.

function kuramoto_one(dtheta, theta, p, t)
    dtheta[1] = 1 + (1/3) * (0) # Included as a formality. 😜
end

theta0 = [0.0]
tspan = (0.0, 100.0)
prob = ODEProblem(kuramoto_one,
                  theta0,
                  tspan)

sol = solve(prob)

plot(sol)

# If I've coded this right, we should literally just get a straight line
# for output. Let's see if that happens.
#
# And it does! Now, the question is can we wrap it around 2pi with a
# modulo operation?
#
# It appears not. The default way that DiffEq solves this actually only
# gives us a few data points, since it (correctly) guesses somehow that
# the output is going to be very simple. I wonder if we can adjust that
# somehow.

sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)

# That didn't work.
# Aha! I see what we're looking for now. We want a "fixed timestep"
# method. The Euler method is the canonical choice, of course.

fixed_timestep=0.1

sol = solve(prob,Euler(),dt=fixed_timestep)

# Now we need to iterate through sol.u and modulo everything by 2pi.
# Oh cool, there's actually an exact mod2pi function in Base.Math --
# which, apparently, doesn't even need to be imported? Okay. I guess
# anything Base.X is already there, and the X is just to keep things
# organized.

for x in sol.u
    print(mod2pi(x[1]), "\n")
    x[1] = mod2pi(x[1])
end

plot(sol)

# Alright, now THAT is what we're looking for. 😄

########################################################################
# Okay, now let's try to implement this cleanly.

function kuramoto_one_cleaned_up(dtheta, theta, p, t)
    dtheta[1] = 1
    theta[1] = mod2pi(theta[1]) # This seems to work fine.
end

theta0 = [0.0]         # For type reasons, the .0 is significant.
tspan = (0.0, 100.0)    # Otherwise Julia expects the exactness of Ints.
prob = ODEProblem(kuramoto_one_cleaned_up, theta0, tspan)
solution = solve(prob, Euler(), dt=0.01)
plot(solution)

# Great.
########################################################################
# Now let's try an N=2 Kuramoto model. Again, hard coding it, because
# I am not super comfortable yet in Julia.

# K is our coupling constant. Raising it will make our values converge
# faster, and it should be a floating point. Setting this to 0.0 means
# no stabilization whatsoever should occur, but so long as it's non-
# negative, it should work fine. (I guess negative would be like, a
# destabilizing force? Haha!)
K = 0.1

# N is the number of oscillators we have paired. It should stay at 2
# for now.
N = 2

# omega_1 and omega_2 are the "natural frequencies" of our two
# oscillators. They can be whatever value you want. Keep them as
# floats, though, just to be safe.
omega_1 = 1.0
omega_2 = 1.0

function kuramoto_two(dtheta, theta, p, t)
    # First, let's make sure our thetas are wrapped around.
    theta[1] = mod2pi(theta[1])
    theta[2] = mod2pi(theta[2])

    # Now, let's implement the actual Diff EQs.
    dtheta[1] = omega_1 + (K / N) * (sin(theta[2] - theta[1]))
    dtheta[2] = omega_2 + (K / N) * (sin(theta[1] - theta[2]))

    # Finally, we'll wrap around again for good measure.
    theta[1] = mod2pi(theta[1])
    theta[2] = mod2pi(theta[2])
end

theta0 = [0.0, 0.0]
tspan = (0.0, 100.0)
prob = ODEProblem(kuramoto_two, theta0, tspan)
solution = solve(prob, Euler(), dt=0.01)
plot(solution)
# And as you can see, the two are exactly the same.
# Which makes sense - they have the same natural frequencies, and
# the same starting positions.
#
# But what if we vary up the initial conditions a bit?
########################################################################
theta0 = [mod2pi(rand()*(2*pi)), mod2pi(rand()*(2*pi))]
tspan = (0.0, 10.0)
prob = ODEProblem(kuramoto_two, theta0, tspan)
solution = solve(prob, Euler(), dt=0.01)
plot(solution)
# Eheheheheh, *now* we start to see it at work. 😈
########################################################################
# Okay, so that's pretty cool. But let's say we wanted to Monte Carlo
# this up. How might we do that?

# That's where things get a bit tricky. I'll have to read up on this.
# http://docs.juliadiffeq.org/latest/features/ensemble.html#Parallel-Ensemble-Simulations-1

# I'll leave it here for tonight.
########################################################################


text:
# Pressurized hyperboloid clamped around the ends.
#
# Example introduced in
# @article{Lee2004,
#    author = {Lee, P. S. and Bathe, K. J.},
#    title = {Development of MITC isotropic triangular shell finite elements},
#    journal = {Computers & Structures},
#    volume = {82},
#    number = {11-12},
#    pages = {945-962},
#    ISSN = {0045-7949},
#    DOI = {10.1016/j.compstruc.2004.02.004},
#    year = {2004},
#    type = {Journal Article}
# }

module cos_2t_press_hyperboloid_fixed_examples

using LinearAlgebra
using FinEtools
using FinEtools.MeshModificationModule: distortblock
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FESetShellQ4Module: FESetShellQ4
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json
using FinEtools.MeshExportModule.VTKWrite: vtkwrite

# Parameters:
E = 2.0e5
nu = 1/3;
pressure = 1.0;
Length = 2.0;

# The hyperboloid axis is parallel to Y

function hyperbolic!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    n = cross(tangents[:, 1], tangents[:, 2])
    n = n/norm(n)
    # r = vec(XYZ); r[2] = 0.0
    csmatout[:, 3] .= n
    csmatout[:, 2] .= (0.0, 1.0, 0.0)
    cross3!(view(csmatout, :, 1), view(csmatout, :, 2), view(csmatout, :, 3))
    return csmatout
end

function computetrac!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    r = vec(XYZ); r[2] = 0.0
    r .= vec(r)/norm(vec(r))
    theta = atan(r[3], r[1])
    n = cross(tangents[:, 1], tangents[:, 2])
    n = n/norm(n)
    forceout[1:3] = n*pressure*cos(2*theta)
    forceout[4:6] .= 0.0
    # @show dot(n, forceout[1:3])
    return forceout
end

function _execute(formul, n = 8, thickness = Length/2/100, visualize = false, distortion = 0.0)
    tolerance = Length/n/100
    fens, fes = distortblock(T3block, 90/360*2*pi, Length/2, n, n, distortion, distortion);
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        a=fens.xyz[i, 1]; y=fens.xyz[i, 2];
        R = sqrt(1 + y^2)
        fens.xyz[i, :] .= (R*sin(a), y, R*cos(a))
    end

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    stiffness = formul.stiffness
    associategeometry! = formul.associategeometry!

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # plane of symmetry perpendicular to Z
    l1 = selectnode(fens; box = Float64[-Inf Inf -Inf Inf 0 0], inflate = tolerance)
    for i in [3,4,5]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to Y
    l1 = selectnode(fens; box = Float64[-Inf Inf 0 0 -Inf Inf], inflate = tolerance)
    for i in [2,4,6]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to X
    l1 = selectnode(fens; box = Float64[0 0 -Inf Inf -Inf Inf], inflate = tolerance)
    for i in [1,5,6]
        setebc!(dchi, l1, true, i)
    end
    # clamped edge perpendicular to Y
    l1 = selectnode(fens; box = Float64[-Inf Inf Length/2 Length/2 -Inf Inf], inflate = tolerance)
    for i in [1,2,3,4,5,6]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    # Midpoint of the fixed edge
    # nl = selectnode(fens; box = Float64[R R L/2 L/2 -Inf Inf], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(fes, TriRule(3)))

    fi = ForceIntensity(FFlt, 6, computetrac!);
    F = distribloads(lfemm, geom0, dchi, fi, 2);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    strainenergy = 1/2 * U' * K * U
    @info "Strain Energy: $(round(strainenergy, digits = 9))"

    # Generate a graphical display of resultants
    ocsys = CSys(3, 3, hyperbolic!)
    scalars = []
    for nc in 1:3
        fld = fieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
        push!(scalars, ("m$nc", fld.values))
    end
    vtkwrite("cos_2t_press_hyperboloid_fixed-$(n)-$(thickness)-$(distortion)-m.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
    scalars = []
    for nc in 1:3
        fld = fieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
        push!(scalars, ("n$nc", fld.values))
    end
    vtkwrite("cos_2t_press_hyperboloid_fixed-$(n)-$(thickness)-$(distortion)-n.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
    scalars = []
    for nc in 1:2
        fld = fieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
        push!(scalars, ("q$nc", fld.values))
    end
    vtkwrite("cos_2t_press_hyperboloid_fixed-$(n)-$(thickness)-$(distortion)-q.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])

    # Visualization
    if visualize
        scattersysvec!(dchi, (Length/8)/maximum(abs.(U)).*U)
        update_rotation_field!(Rfield0, dchi)
        plots = cat(plot_space_box([[0 0 -Length/2]; [Length/2 Length/2 Length/2]]),
            #plot_nodes(fens),
            plot_midsurface(fens, fes; x = geom0.values, facecolor = "rgb(12, 12, 123)"),
            plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
            dims = 1)
        pl = render(plots)
    end

    return strainenergy
end

function test_convergence(formul, thicknessmult = 1/100, distortion = 0.0)
    @info "Pressurized Hyperbolic shell, fixed ends, formulation=$(formul)"
    results = []
    ns = [16, 32, 64, 128, 256]
    for n in ns
        push!(results, _execute(formul, n, Length/2*thicknessmult, false, 2*distortion/n))
    end
    return ns, results
end

end # module

using .cos_2t_press_hyperboloid_fixed_examples
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using PGFPlotsX

let

    for distortion in [2.0, 0.0]
        ns, results100 = cos_2t_press_hyperboloid_fixed_examples.test_convergence(FEMMShellT3FFModule, 1/100, distortion)
        ns, results1000 = cos_2t_press_hyperboloid_fixed_examples.test_convergence(FEMMShellT3FFModule, 1/1000, distortion)
        ns, results10000 = cos_2t_press_hyperboloid_fixed_examples.test_convergence(FEMMShellT3FFModule, 1/10000, distortion)
        ns, results100000 = cos_2t_press_hyperboloid_fixed_examples.test_convergence(FEMMShellT3FFModule, 1/100000, distortion)

        errors100 = diff(results100)/results100[end]
        errors1000 = diff(results1000)/results1000[end]
        errors10000 = diff(results10000)/results10000[end]
        errors100000 = diff(results100000)/results100000[end]

        objects = []

        all_results = [("t=L/100", errors100, "*"), ("t=L/1000", errors1000, "x"), ("t=L/10000", errors10000, "triangle"), ("t=L/100000", errors100000, "diamond"), ]

        for r in  all_results
            @pgf p = PGFPlotsX.Plot(
            {
            color = "black",
            line_width  = 0.7,
            style = "solid",
            mark = "$(r[3])"
            },
            Coordinates([v for v in  zip(1 ./ ns, r[2]) if v[2] !== missing])
            )
            push!(objects, p)
            push!(objects, LegendEntry("$(r[1])"))
        end

        @pgf p = PGFPlotsX.Plot(
        {
        color = "red",
        line_width  = 0.7,
        style = "dashed"
        },
        Coordinates([(0.01, 0.0001), (0.1, 0.01)])
        )
        push!(objects, p)
        push!(objects, LegendEntry("slope=2"))

        @pgf ax = Axis(
        {
        xlabel = "Relative element size [ND]",
        ylabel = "Approximate Normalized Error [ND]",
        xmin = 0.005,
        # xmax = range[2],
        ymin = 0.0001,
        xmode = "log",
        ymode = "log",
        yminorgrids = "true",
        grid = "both",
        legend_pos = "south east"
        },
        objects...
        )

        display(ax)
        pgfsave("cos_2t_press_hyperboloid_fixed_$(distortion)-convergence.pdf", ax)

    end
end


text:
using RecursiveSparseBlocks

function run(T)
    out = open("benchmark.csv", "w")
    println(out, "method,n,p,t")

    N = 100
    for (maxn, p) in [(1000000, 1e-4)]
        for _ in 1:N
            n = rand(1:maxn)
            A = sprand(T, n, n, p)
            Arsb = SparseMatrixRSB(A)
            x = rand(T, n)

            y = A*x
            tic()
            y = A*x
            t = toc()
            println(out, "SparseMatrixCSC,", n, ",", p, ",", t)

            y = Arsb*x
            tic()
            y = Arsb*x
            t = toc()
            println(out, "SparseMatrixRSB,", n, ",", p, ",", t)
        end
    end

    close(out)
end

run(Float64)

#@time bench_csc()
#@time bench_csc()
#@profile bench_csc()
#Profile.print(C=true)

#@time bench_rsb()
#@time bench_rsb()
#@profile bench_rsb()
#Profile.print(C=true)


text:
using FlightMechanicsUtils
using Test

@testset "Ellipsoid WGS84" begin
    # Rogers, R. M. (2007). Applied mathematics in integrated navigation systems.
    # American Institute of Aeronautics and Astronautics. (Page 77, table 4.2)
    @test isapprox(WGS84.a, 6378137)
    @test isapprox(WGS84.b, 6356752.3142)
    @test isapprox(WGS84.finv, 298.257223563)
    @test isapprox(WGS84.e2, 0.0818191908426^2)
end

@testset "llh <-> ECEF" begin
    # llh ECEF (using data from
    # - [1] Bowring, B. R. (1976). Transformation from spatial to geographical
    # coordinates. Survey review, 23(181), 323-327.)

    # From [1] Example 1 (page 325)
    x = 4114496.258  # m
    y = 0.0          # m
    z = 4870157.031  # m
    # Ellipsoid does not match any known standard
    a_bowring1976 = 6378249.145  # m
    b_bowring1976 = 6356514.870  # m
    finv_bowring1976 = a_bowring1976 / (a_bowring1976 - b_bowring1976)
    Bowring1976 = Ellipsoid(a_bowring1976, finv_bowring1976, "Bowring 1976")

    exp_llh = [deg2rad(50.0), 0, 10000.0]
    llh = ecef2llh(x, y, z; ellipsoid=Bowring1976)
    @test isapprox(llh, exp_llh)

    exp_xyz = [4114496.258, 0.0, 4870157.031]  # m
    xyz = llh2ecef(deg2rad(50.0), 0, 10000.0; ellipsoid=Bowring1976)
    @test isapprox(xyz, exp_xyz)

    # http://www.sysense.com/products/ecef_lla_converter/index.html
    x, y, z = (4114496.258, 0.0, 4870157.031)  # m
    exp_llh = [deg2rad(49.996908), 0.000000, 9907.31]
    llh = ecef2llh(x, y, z)
    # Longitude should be exact, latitude max error: 0.0018" according to [1]
    # Altitude 0.17 m
    # The implemented function must be better, but I don't know which one is using
    # this web page. The ellipsoid is assumed to be WGS84
    @test isapprox(llh[:2], exp_llh[:2], atol=deg2rad(0.0018/3600.0))
    @test isapprox(llh[3], exp_llh[3], atol=0.17)

    exp_xyz = [4114291.97, 0.00, 4870449.48]  # m
    xyz = llh2ecef(deg2rad(50.0), 0, 10000.0)
    @test isapprox(xyz, exp_xyz)

    # non-zero longitude
    x, y, z = (3814496.258, 1514496.258, 4870157.031)  # m
    exp_llh = [deg2rad(50.068095), deg2rad(21.654910), 3263.97]
    llh = ecef2llh(x, y, z)
    # Longitude should be exact, latitude max error: 0.0018" according to [1]
    # Altitude 0.17 m
    # The implemented function must be better, but I don't know which one is using
    # this web page. The ellipsoid is assumed to be WGS84
    @test isapprox(llh[:2], exp_llh[:2], atol=deg2rad(0.0018/3600.0))
    @test isapprox(llh[3], exp_llh[3], atol=0.17)

    exp_xyz = [3814496.258, 1514496.258, 4870157.031]  # m
    xyz = llh2ecef(deg2rad(50.068095), deg2rad(21.654910), 3263.97)
    @test isapprox(xyz, exp_xyz)
end


text:
# Load packages
using CSV
using DataFramesMeta
using GBIF
using Pipe
using Plots
using SimpleSDMLayers
default(dpi=200)

## Preparation

# Define range for Quebec
# extent(-80.0, -56.0, 44.0, 62.0)
# spatialrange = (left=-80., right=-50., bottom=45., top=65.)
spatialrange = (left=-80., right=-50., bottom=44., top=65.)
lon = (spatialrange.left, spatialrange.right)
lat = (spatialrange.bottom, spatialrange.top)

# Define reference layer
temperature = SimpleSDMPredictor(WorldClim, BioClim, 1; spatialrange...)

# Define species
species = ["Vulpes vulpes", "Tamias striatus", "Marmota monax"]
species_en = ["Red Fox", "Eastern Chipmunk", "Groundhog"]

## Extract data

# Wrap everything in one function
function gimme_mammals_pls(sp::String)
    # Replace spaces
    _sp = replace(sp, " " => "_")

    # Define taxons
    sp_code = taxon(sp, rank = :SPECIES)

    # Query parameters
    observations = occurrences(
        sp_code,
        "limit" => 300,
        "hasCoordinate" => "true",
        "decimalLatitude" => lat,
        "decimalLongitude" => lon,
    )

    # Retrieve all
    while length(observations) < size(observations)
        occurrences!(observations)
    end

    # Check the data
    obs_df = DataFrame(observations)

    ## Layers

    # Plot to see occurrences
    p_obs = contour(temperature; c=:cork, frame=:box, fill=true, clim=(-20, 20), levels=6, title=sp)
    scatter!(
        longitudes(observations), latitudes(observations); lab="", c=:white, msc=:orange, ms=3
    )
    savefig(joinpath("figures", "obs_$(_sp).png"))

    # Layer with the number of observations per cell
    p_counts = counts = mask(temperature, observations, Float32)
    plot(log1p(counts); c=:tokyo, title=sp)
    savefig(joinpath("figures", "counts_$(_sp).png"))

    # Check the number of pixels with observations
    # presabs = replace(counts, 0.0 => nothing)
    presabs = @pipe mask(temperature, observations, Bool) |>
        convert(Float32, _) |>
        replace(_, 0.0 => nothing)

    # Plot
    p_presabs = plot(temperature, c=:lightgrey)
    plot!(presabs; c=:BuPu, cb=:none, title=sp)
    savefig(joinpath("figures", "presabs_$(_sp).png"))

    # Export elements
    results = (
        occurrences = occurrences,
        obs_df = obs_df,
        p_obs = p_obs,
        counts = counts,
        p_counts = p_counts,
        presabs = presabs,
        p_presabs = p_presabs
    )

    return results
end

# Get the data for the three species
sp1 = gimme_mammals_pls(species[1]);
sp2 = gimme_mammals_pls(species[2]);
sp3 = gimme_mammals_pls(species[3]);

# Prepare a prettier figure as an exemple
foxplot = sp1.p_obs
title!(foxplot, "Observations du Renard roux")
xaxis!("Longitude")
yaxis!("Latitude")
plot!(colorbar_title="Température moyenne")
savefig(joinpath("figures", "exemple_renard.png"))

## Assemble dataset

# Combine observations in single DataFrame
combined_df = vcat(sp1.obs_df, sp2.obs_df, sp3.obs_df)
unique(combined_df.name)

# Export
CSV.write(joinpath("data", "mammals_occurrences.csv"), combined_df)

# Assemble layers
layers = [sp1.presabs, sp2.presabs, sp3.presabs]

# Export
geotiff(joinpath("data", "mammals_layers.csv"), layers)

# Convert layers to DataFrame
layers_df_full = DataFrame(layers)
rename!(layers_df_full, "x1" => "sp1", "x2" => "sp2", "x3" => "sp3")

# Remove sites without observations
layers_df = filter(x -> any(!ismissing, [x.sp1, x.sp2, x.sp3]), layers_df_full)

# Verify species co-occurrence
layers_df.sum = map(x -> sum(skipmissing([x.sp1, x.sp2, x.sp3])), eachrow(layers_df))
layers_df
unique(layers_df.sum)
@chain begin layers_df
    groupby(:sum)
    combine(nrow)
end

# Verify the occupancy ratio & pixel count
@chain layers_df begin
    describe(:min, :max, :nmissing)
    @rtransform(:n = nrow(layers_df) - :nmissing)
    @rsubset(String(:variable) in ["sp1", "sp2", "sp3"])
    @rtransform(:ratio = :n / nrow(layers_df))
    @transform(:equal_layer = :n .== length.(layers))
end
# Seems reasonable

## Combine with environmental data

# Load glossary
glo_path = "$(homedir())/github/betadiversity-hotspots/data/proc/glossary.csv"
glossary = CSV.read(glo_path, DataFrame)
filter!(:type => !=("species"), glossary)
glossary = vcat(
    filter(:type => ==("climate"), glossary),
    filter(:type => ==("landcover"), glossary)
)

# Add new species
glossary_sp = DataFrame(
    variable = ["sp1", "sp2", "sp3"],
    type = "species",
    full_name = replace.(species, " " => "_"),
    description = species_en
)
glossary = vcat(glossary_sp, glossary)

# Load env data
nlayers = nrow(filter(:type => !=("species"), glossary))
env_path = "$(homedir())/github/betadiversity-hotspots/data/raster/env_stack.tif"
env_layers = [geotiff(SimpleSDMPredictor, env_path, i; spatialrange...) for i in 1:nlayers]

# Verify the data
size(env_layers[1]) == size(layers[1])
boundingbox(env_layers[1]) == boundingbox(layers[1])
env_layers[1] == temperature

# Arrange as DataFrame
env_df = DataFrame(env_layers)
env_names = filter(:type => !=("species"), glossary).variable
rename!(env_df, ["longitude", "latitude", env_names...])

# Combine with species data
full_df = leftjoin(select(layers_df, Not(:sum)), env_df, on=[:longitude, :latitude])

# Export
CSV.write("./data/mammals_complete.csv", full_df)
CSV.write("./data/mammals_env.csv", env_df)
CSV.write("./data/mammals_glossary.csv", glossary)


text:
using ConstrainedDynamics
using ConstrainedDynamics: skewplusdiag

function dynTvel(vars)
    ezg = [0; 0; -9.81]
    Δt = 0.01
    m = 1.

    vc1 = vars[1:3]
    vc2 = vars[4:6]

    return m * ((vc2 - vc1) / Δt + ezg)
end

function dynRvel(vars)
    Δt = 0.01

    J = diagm(ones(3))
    ωc1 = vars[1:3]
    ωc2 = vars[4:6]
    sq2 = sqrt(4 / Δt^2 - ωc2' * ωc2)
    sq1 = sqrt(4 / Δt^2 - ωc1' * ωc1)

    return skewplusdiag(ωc2, sq2) * (J * ωc2) - skewplusdiag(ωc1, sq1) * (J * ωc1)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ f64b7897-2e6e-4654-bb84-43bf3807f771
html"<button onclick='present()'>presentation mode</button>"

# ╔═╡ 4752d352-f864-11eb-3cdc-eff915ddae6e
md"""
# Some packages for Data Science
- My personal recommendations of Julia packages that can be useful to Data Scientists
    - Some packages are mature - some aren't
    - Some packages have volunteer maintainers - others are developed within companies
    - Most data science packages are still at the "volunteer maintainer" stage
- Areas I will mention
    - Data manipulation
    - Exploratory/presentation graphics
    - Model fitting (ML, Bayesian inference)
## Data manipulation
- For column-oriented tables, [DataFrames](https://github.com/JuliaData/DataFrames.jl) is the workhorse.
    - has reached a 1.0.0 release
    - has a mini-language for split-apply-combine types of operations
    - provides various types of joins of tables - a lot of work has gone into performance enhancements here
    - language, while consistent, can be a bit awkward
- I prefer the [DataFrameMacros](https://github.com/jkrumbiegel/DataFrameMacros.jl) language in combination with [Chain](https://github.com/jkrumbiegel/Chain.jl)
- The [Tables](https://github.com/JuliaData/Tables.jl) package provides a clearing house for various representations of column-oriented or row-oriented table representations.
    - provides a general formulation of `rowtable` (`Vector` of `NamedTuple`s) or `columntable` (`NamedTuple` of `Vector`s) and acts as a hub for transfers between them.
        + e.g. suppose I read a csv file using the [CSV](https://github.com/JuliaData/CSV.jl) package, how do I create a DataFrame from the result or from an [Arrow](https://github.com/JuliaData/Arrow.jl.git) table or write it out to an [SQLite](https://github.com/JuliaDatabases/SQLite.jl.git) table or ...?
- Many of the JuliaData packages are primarily written by Jacob Quinn who has an amazing [tutorial](https://www.youtube.com/watch?v=uLhXgt_gKJc) on microservices from JuliaCon 2020.
## Exploratory/presentation graphics
- There are many graphics systems in Julia.
- I prefer the ones based on [Makie](https://github.com/JuliaPlots/Makie.jl.git) - check out the gallery at [BeautifulMakie](https://lazarusa.github.io/BeautifulMakie/).
    - the calls to produce graphics elements and to lay them out are common to all the back ends (that's what Makie provides)
    - the actual display requires a back end package like [CairoMakie](https://github.com/JuliaPlots/Makie.jl.git), which we have used, or [GLMakie](https://github.com/JuliaPlots/Makie.jl.git), based on OpenGL
    - a higher-level (but still experimental) approach is available in [AlgebraOfGraphics](https://github.com/JuliaPlots/AlgebraOfGraphics.jl.git)
## Model fitting
- For machine learning many people like [Flux](https://github.com/FluxML/Flux.jl.git) (I haven't used it myself).  See also [MLJ](https://github.com/alan-turing-institute/MLJ.jl.git)
- [Turing](https://github.com/TuringLang/Turing.jl.git) is a popular framework for Bayesian inference. [Soss](https://github.com/cscherrer/Soss.jl.git) and associated packages like [MeasureTheory](https://github.com/cscherrer/MeasureTheory.jl.git) are interesting alternatives still under development.
- Many of these packages benefit from extensive development in Julia on topics like Automatic Differentiation.
- Some basic statistics tools under the [JuliaStats group](https://github.com/JuliaStats/) include [Distributions](https://github.com/JuliaStats/Distributions.jl.git), [StatsModels](https://github.com/JuliaStats/StatsModels.jl.git), and [StatsBase](https://github.com/JuliaStats/StatsBase.jl.git)
    - the organization here is not as good as I would like it to be but it is dependent on volunteers
"""

# ╔═╡ Cell order:
# ╟─f64b7897-2e6e-4654-bb84-43bf3807f771
# ╟─4752d352-f864-11eb-3cdc-eff915ddae6e


text:
function types(typed::DataType)
  println(typed)
  while typed != Any
    typed = supertype(typed)
      println(typed)
  end
end

types(Float16)

max_stars_repo_name: wchmielarz/NobodyDoesItLikeJuliet

max_issues_repo_name: wchmielarz/NobodyDoesItLikeJuliet

max_forks_repo_name: wchmielarz/NobodyDoesItLikeJuliet


text:
import DataAPI: leftjoin, innerjoin

# A set of kwargs that can be provided by the user.
# Used for deciding whether to use the `DTables` join implementation directly
# or to attempt using an external join function.
const JOINKWARGS = Set([
    :l_sorted,
    :r_sorted,
    :r_unique,
    :lookup,
])

"""
    leftjoin(d1::DTable, d2; on=nothing, l_sorted=false, r_sorted=false, r_unique=false, lookup=nothing)

Perform a left join of `d1` with any `Tables.jl` compatible table type.
Returns a `DTable` with the result.

If the underlying table type happens to have a `leftjoin` implementation
and none of the below `DTable` related kwargs will be provided the specialized function will be used.
A good example of that is calling `leftjoin` on a `DTable` with a `DataFrame` underlying type
and a `d2` of `DataFrame` type.

# Keyword arguments

- `on`: Column symbols to join on. Can be provided as a symbol or a pair of symbols in case the column names differ. For joins on multiple columns a vector of the previously mentioned can be provided.
- `l_sorted`: To indicate the left table is sorted - only useful if the `r_sorted` is set to `true` as well.
- `r_sorted`: To indicate the right table is sorted.
- `r_unique`: To indicate the right table only contains unique keys.
- `lookup`: To provide a dict-like structure that will allow for direct matching of inner rows. The structure needs to contain keys in form of a `Tuple` and values in form of type `Vector{UInt}` containing the related row indices.
"""
function leftjoin(d1::DTable, d2; kwargs...)
    f = (l, r, ks) -> _leftjoinwrapper(l, r; ks...)
    v = [Dagger.@spawn f(c, d2, kwargs) for c in d1.chunks]
    DTable(v, d1.tabletype)
end

function leftjoin(d1::DTable, d2::DTable; kwargs...)
    f = (l, r, ks) -> _join(:leftjoin, l, r; ks...)
    v = [Dagger.@spawn f(c, d2, kwargs) for c in d1.chunks]
    DTable(v, d1.tabletype)
end

function leftjoin(d1::GDTable, d2; kwargs...)
    d = leftjoin(d1.dtable, d2; kwargs...)
    GDTable(d, d1.cols, d1.index)
end

function _leftjoinwrapper(l, r; kwargs...)
    if !any(k in JOINKWARGS for k in keys(kwargs)) && use_dataframe_join(typeof(l), typeof(r))
        leftjoin(l, r; kwargs...)
    else
        _join(:leftjoin, l, r; kwargs...)
    end
end

"""
    innerjoin(d1::DTable, d2; on=nothing, l_sorted=false, r_sorted=false, r_unique=false, lookup=nothing)

Perform an inner join of `d1` with any `Tables.jl` compatible table type.
Returns a `DTable` with the result.

If the underlying table type happens to have a `innerjoin` implementation
and none of the below `DTable` related kwargs will be provided the specialized function will be used.
A good example of that is calling `innerjoin` on a `DTable` with a `DataFrame` underlying type
and a `d2` of `DataFrame` type.

# Keyword arguments

- `on`: Column symbols to join on. Can be provided as a symbol or a pair of symbols in case the column names differ. For joins on multiple columns a vector of the previously mentioned can be provided.
- `l_sorted`: To indicate the left table is sorted - only useful if the `r_sorted` is set to `true` as well.
- `r_sorted`: To indicate the right table is sorted.
- `r_unique`: To indicate the right table only contains unique keys.
- `lookup`: To provide a dict-like structure that will allow for direct matching of inner rows. The structure needs to contain keys in form of a `Tuple` and values in form of type `Vector{UInt}` containing the related row indices.
"""
function innerjoin(d1::DTable, d2; kwargs...)
    f = (l, r, ks) -> _innerjoinwrapper(l, r; ks...)
    v = [Dagger.@spawn f(c, d2, kwargs) for c in d1.chunks]
    DTable(v, d1.tabletype)
end

function innerjoin(d1::DTable, d2::DTable; kwargs...)
    f = (l, r, ks) -> _join(:innerjoin, l, r; ks...)
    v = [Dagger.@spawn f(c, d2, kwargs) for c in d1.chunks]
    DTable(v, d1.tabletype)
end

function innerjoin(d1::GDTable, d2; kwargs...)
    d = innerjoin(d1.dtable, d2; kwargs...)
    GDTable(d, d1.cols, d1.index)
end

function _innerjoinwrapper(l, r; kwargs...)
    if !any(k in JOINKWARGS for k in keys(kwargs)) && use_dataframe_join(typeof(l), typeof(r))
        innerjoin(l, r; kwargs...)
    else
        _join(:innerjoin, l, r; kwargs...)
    end
end

"""
    match_inner_indices(l, r, cmp_l, cmp_r, lookup, r_sorted, l_sorted, r_unique)

Function responsible for picking the optimal method of joining inner indices depending on the
additional information about the tables provided by the user.
"""
function match_inner_indices(l, r, cmp_l, cmp_r, lookup, r_sorted, l_sorted, r_unique)
    if lookup !== nothing
        match_inner_indices_lookup(l, lookup, cmp_l) # uses the `lookup` to find indices
    elseif r_sorted && l_sorted
        match_inner_indices_lsorted_rsorted(l, r, cmp_l, cmp_r, r_unique) # loop through r once
    elseif r_unique
        match_inner_indices_runique(l, r, cmp_l, cmp_r) # break on first match
    elseif r_sorted
        match_inner_indices_rsorted(l, r, cmp_l, cmp_r) # break on last match
    else # generic fallback, no optimization
        match_inner_indices(l, r, cmp_l, cmp_r)
    end
end

"""
    _join(type::Symbol, l_chunk, r; kwargs...)

Low level join method for `DTable` joins using the generic implementation.
It joins an `l_chunk` with `r` assuming `r` is a continuous table.
"""
function _join(
        type::Symbol,
        l_chunk,
        r;
        on=nothing,
        l_sorted=false,
        r_sorted=false,
        r_unique=false,
        lookup=nothing
    )

    names, _, other_r, cmp_l, cmp_r = resolve_colnames(l_chunk, r, on)

    inner_l, inner_r = match_inner_indices(l_chunk, r, cmp_l, cmp_r, lookup, r_sorted, l_sorted, r_unique)

    outer_l = type == :innerjoin ? Set{UInt}() : find_outer_indices(l_chunk, inner_l)
    build_joined_table(type, names, l_chunk, r, inner_l, inner_r, outer_l, other_r)
end

"""
    _join(type::Symbol, l_chunk, r::DTable; kwargs...)

Low level join method for `DTable` joins using the generic implementation.
It joins an `l_chunk` with `r` assuming `r` is a `DTable`.
In this case the join is split into multiple joins of `l_chunk` with each chunk of `r` and a final merge operation.
"""
function _join(
        type::Symbol,
        l_chunk,
        r::DTable;
        on=nothing,
        l_sorted=false,
        r_sorted=false,
        r_unique=false,
        lookup=nothing
    )

    names, _, other_r, cmp_l, cmp_r = resolve_colnames(l_chunk, r, on)

    process_one_chunk = (type, l, r, cmp_l, cmp_r, other_r, lookup, r_sorted, l_sorted, r_unique) -> begin
        inner_l, inner_r = match_inner_indices(l, r, cmp_l, cmp_r, lookup, r_sorted, l_sorted, r_unique)
        inner_chunk = Dagger.tochunk(build_joined_table(type, names, l, r, inner_l, inner_r, Set{UInt}(), other_r))
        outer_l = type == :innerjoin ? Set{UInt}() : find_outer_indices(l, inner_l)
        return inner_chunk, outer_l
    end

    vs = [Dagger.@spawn process_one_chunk(type, l_chunk, chunk, cmp_l, cmp_r, other_r, lookup, r_sorted, l_sorted, r_unique) for chunk in r.chunks]

    to_merge = Vector{Chunk}()
    sizehint!(to_merge, length(r.chunks))

    v = fetch.(vs)
    append!(to_merge, getindex.(v, 1))

    if type == :leftjoin
        outer_l = intersect(getindex.(v, 2)...)
        inner_l = inner_r = Vector{UInt}()
        outer = Dagger.tochunk(build_joined_table(type, names, l_chunk, r, inner_l, inner_r, outer_l, other_r))
        push!(to_merge, outer)
    end

    merge_chunks(Tables.materializer(l_chunk), to_merge)
end

"""
    use_dataframe_join(d1type, d2type)

Determines whether to use the DataAPI join function, which leads to usage of DataFrames join function if both types are `DataFrame`.
Remove this function and it's usage once a generic Tables.jl compatible join function becomes available.
Porting the Dagger join functions to TableOperations is an option to achieve that.
"""
function use_dataframe_join(d1type, d2type)
    :DataFrame == d1type.name.name == d2type.name.name
end