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

using Markdown
using InteractiveUtils

# ╔═╡ 32ae409c-9ed9-11eb-3ce8-2bef48bebcae
using CSV, DataFrames, Plots

# ╔═╡ 6779c065-2fa6-4630-8263-65ce31f1488f
html"
Franck-Hertz Experiment

Calculating the first excited energy of mercury

Benjamin Jiron

April 2021

Junior Lab

"

# ╔═╡ c04d9cd2-684f-466a-af91-d61bea5761dc
begin
df = CSV.File("FranckHertzData.csv") |> DataFrame
df[1:5, :]
end

# ╔═╡ 8c24fdc4-77d7-4d97-a617-f1cf12c9c469
let
plot(df[!, "Voltage"], df[!, "180C"], label="T=180°C", lw=2)
plot!(df[!, "Voltage"], df[!, "195C"], label="T=195°C", lw=2)
plot!(df[!, "Voltage"], df[!, "210C"], label="T=210°C", lw=2)
plot!(title="Anode Current vs. Voltage", xlabel="Voltage (V)", ylabel="Current (A)",
	legend=:topleft)
end

# ╔═╡ 5ddc7b4a-4060-4b12-96c9-3e8ab707b5a7
html"
There's a lot of noise before Voltage=10V so let's just drop that part and analyze the rest:
"

# ╔═╡ fc7f0733-7d93-474c-9688-73a02c6a51cc
begin
T = df[df[!, "Voltage"] .> 10, :] 	# Filtering
T₁ = select(T, "Voltage", "180C") 	# ⇒ :lowtemp
T₂ = select(T, "Voltage", "195C") 	# ⇒ :medtemp
T₃ = select(T, "Voltage", "210C") 	# ⇒ :hightemp
rename!(T₁, "180C" => "Current") 	# Renaming columns
rename!(T₂, "195C" => "Current")
rename!(T₃, "210C" => "Current")
T₁[1:5, :], T₂[1:5, :], T₃[1:5, :]
end

# ╔═╡ d30c9e49-0d28-4882-93fd-9354dce62432
html"

We expect the energy difference between the first and second excited states of mercury to be about equal to the mean ΔE between the peaks of the curves.

Let's set up some naïve filtering functions to try and pick out maximum and minimum currents, and two corrective functions in case we miss:
"

# ╔═╡ 5e945e11-869f-4040-af71-c79dd0351f60
function argminima(v::Vector)::Vector
	N = length(v)
	minima = []
	for i ∈ 2:(N-1)
		(v[i] ≤ v[i-1] && v[i] < v[i+1]) && push!(minima, i)
	end
	minima
end

# ╔═╡ f5b8520b-0d99-4e8a-8460-32dc1420deb3
function argmaxima(v::Vector)::Vector
	N = length(v)
	maxima = []
	for i ∈ 2:(N-1)
		(v[i] ≥ v[i-1] && v[i] > v[i+1]) && push!(maxima, i)
	end
	maxima
end

# ╔═╡ fb4b04e8-5901-4a7c-bd6d-19c4201fbf01
function correctminima!(iminima::Vector, v::Vector)
	corrections = []
	k = 10 		# We're going to check k nearest neighbors for local minima
	for i ∈ iminima
		nearbypoints = v[i-k:i+k]
		correction = argmin(nearbypoints) - (k+1) # set correction = 0 if not needed
		push!(corrections, correction)
	end
	iminima .+= corrections
	unique!(iminima)
end

# ╔═╡ 4b6207e2-e835-4c26-a252-734986b31f14
function correctmaxima!(imaxima::Vector, v::Vector)
	corrections = []
	k = 10 		# k nearest neighbors on each side
	for i ∈ imaxima
		nearbypoints = v[i-k:i+k]
		correction = argmax(nearbypoints) - (k+1) # correction = 0 if not needed
		push!(corrections, correction)
	end
	imaxima .+= corrections
	unique!(imaxima)
end

# ╔═╡ 9053e773-18fc-43f4-b54b-872a4d22dcb8
html"
We'll use these two Dicts, allminima and allmaxima to keep track of our extrema:
"

# ╔═╡ c30ee767-87cd-44e4-a7a2-a752828c80d4
begin
allminima = Dict()
allmaxima = Dict()
end

# ╔═╡ 8e681f02-2f7b-49f5-b10b-fd448141dd93
let
# Get indices of maxima using our function
imaxima = argmaxima(T₁.Current)
# Make minor corrections to the indices found
correctmaxima!(imaxima, T₁.Current)
# Get corresponding voltages
maxvoltages = T₁[imaxima, :Voltage]
# Repeat for the minima
iminima = argminima(T₁.Current)
correctminima!(iminima, T₁.Current)
minvoltages = T₁[iminima, :Voltage]
# Record unique minima/maxima in the dictionaries
allminima[:lowtemp] = minvoltages
allmaxima[:lowtemp] = maxvoltages
# Plotting
plot(T₁.Voltage, T₁.Current*10^9, label=false, legend=:topleft, lw=2)
vline!(minvoltages, linestyle=:dash, label="Calculated Minima", lw=2)
vline!(maxvoltages, linestyle=:dot, label="Calculated Maxima", lw=2)
plot!(title="Current vs Voltage (180°C)", xlabel="Volts", ylabel="Nanoamperes")
end

# ╔═╡ 0d8d30ab-20dd-4d53-9d8e-ddaa2d1d6182
html"
That naïve filtering function worked really well on the first dataset, but the next two need some adjusting:
"

# ╔═╡ 332f7f3e-af9b-4d2b-8e22-551d3cb49213
let
imaxima = argmaxima(T₂.Current)
correctmaxima!(imaxima, T₂.Current)
maxvoltages = T₂[imaxima, :Voltage]
iminima = argminima(T₂.Current)
correctminima!(iminima, T₂.Current)
minvoltages = T₂[iminima, :Voltage]
# New step: filter out bad data (the region below 18V is pretty noisy)
filter!(>(18), minvoltages) 	# some trial & error to find these filters
filter!(≠(35.3), minvoltages) 	# if you miss TOO badly the correction function fails
filter!(>(15), maxvoltages)
allminima[:medtemp] = minvoltages
allmaxima[:medtemp] = maxvoltages
plot(T₂.Voltage, T₂.Current*10^9, label=false, legend=:topleft, lw=2, c=:orange)
vline!(minvoltages, linestyle=:dash, label="Calculated Minima", lw=2)
vline!(maxvoltages, linestyle=:dot, label="Calculated Maxima", lw=2)
plot!(title="Current vs Voltage (195°C)", xlabel="Volts", ylabel="Nanoamperes")
end

# ╔═╡ 2dd9a2c7-5279-4c5f-9c02-9f006f8fff3e
let
imaxima = argmaxima(T₃.Current)
correctmaxima!(imaxima, T₃.Current)
maxvoltages = T₃[imaxima, :Voltage]
iminima = argminima(T₃.Current)
correctminima!(iminima, T₃.Current)
minvoltages = T₃[iminima, :Voltage]
filter!(>(20), minvoltages) 			# a bit more trial and error
filter!(∉([21.4, 23.6]), minvoltages) 	# Specific to this dataset & filtering func
filter!(>(20.2), maxvoltages)
allminima[:hightemp] = minvoltages
allmaxima[:hightemp] = maxvoltages
plot(T₃.Voltage, T₃.Current*10^9, label=false, legend=:topleft, lw=2, c=:green)
vline!(minvoltages, linestyle=:dash, label="Calculated Minima", lw=2)
vline!(maxvoltages, linestyle=:dot, label="Calculated Maxima", lw=2)
plot!(title="Current vs Voltage (210°C)", xlabel="Volts", ylabel="Nanoamperes")
end

# ╔═╡ fc26a4e5-1050-4d4d-97cb-2e2b4262d9cd
html"
Finding Eₐ (energy (eV) of the first excited energy of mercury)

From above: We expect the energy difference between the first and second excited states of mercury to be about equal to the mean ΔE between peaks:"

# ╔═╡ 77f54490-4e74-4ba8-866a-9447a21c5b41
let
mins = allminima[:lowtemp]
plot(T₁.Voltage, T₁.Current*10^9, lw=2, legend=false, xlabel="Volts", ylabel="nA")
vline!(mins[2:5], linestyle=:dash)
plot!(mins[3:4], [20, 20], lw=4) 	# short horizontal line
annotate!(mins[3]+2.3, 22, "ΔE")
end

# ╔═╡ d20116e4-20e8-4f69-b965-924fd5e69027
results = DataFrame(:Temperature => ["180°C", "195°C", "210°C"], :Eₐ => 0.0);

# ╔═╡ 64961f92-6ab1-49d2-9eb2-53bf92f986df
mean(X::Vector) = sum(X) / length(X)

# ╔═╡ 7f7ab7a9-39a9-4954-bbd0-7ad4210c6c1a
html"ΔE for 180°C → ( [ΔEminima], [ΔEmaxima] )"

# ╔═╡ 5fe21528-e4f8-4377-ae03-c0b3dd64ed2f
let
# α and β are dummy variables which represent our min/max voltages, respectively
α = allminima[:lowtemp]
β = allmaxima[:lowtemp]
# Find differences in voltages between peaks
ΔE₁ = [α[i+1] - α[i] for i ∈ 1:(length(α)-1)]
ΔE₂ = [β[i+1] - β[i] for i ∈ 1:(length(β)-1)]
# Recording results (Eₐ is estimated to be the mean of means of ΔE₁, ΔE₂)
results[1, :Eₐ] = mean([mean(ΔE₁), mean(ΔE₂)])
ΔE₁, ΔE₂
end

# ╔═╡ c0453f52-c9d6-4361-be6e-0981f6e62246
html"195°C"

# ╔═╡ a45de7a9-0fa3-44f0-abc6-3b6d4380d6ee
let
α = allminima[:medtemp]
β = allmaxima[:medtemp]
ΔE₁ = [α[i+1] - α[i] for i ∈ 1:(length(α)-1)]
ΔE₂ = [β[i+1] - β[i] for i ∈ 1:(length(β)-1)]
results[2, :Eₐ] = mean([mean(ΔE₁), mean(ΔE₂)])
ΔE₁, ΔE₂
end

# ╔═╡ 0cd63574-2afc-43e1-a02f-4ec8eb3d70bc
html"210°C"

# ╔═╡ 39b549f1-7ef6-41e8-9211-ea1266b347c5
let
α = allminima[:hightemp]
β = allmaxima[:hightemp]
ΔE₁ = [α[i+1] - α[i] for i ∈ 1:(length(α)-1)]
ΔE₂ = [β[i+1] - β[i] for i ∈ 1:(length(β)-1)]
results[3, :Eₐ] = mean([mean(ΔE₁), mean(ΔE₂)])
ΔE₁, ΔE₂
end

# ╔═╡ 78669699-861e-4e3b-b405-88088b303285
html"
Thus:
"

# ╔═╡ b869706a-ea91-4e2d-884a-dba0a84334e2
html"An accepted value of Eₐ is 4.86 eV:"

# ╔═╡ b41620cf-1f7f-41d5-b086-25024fab8598
let
Eₐ = 4.86 	# eV
results[!, :AbsoluteError] = results.Eₐ .- Eₐ
results[!, :PercentError] = 100results.AbsoluteError / Eₐ
results
end

# ╔═╡ Cell order:
# ╟─6779c065-2fa6-4630-8263-65ce31f1488f
# ╠═32ae409c-9ed9-11eb-3ce8-2bef48bebcae
# ╠═c04d9cd2-684f-466a-af91-d61bea5761dc
# ╠═8c24fdc4-77d7-4d97-a617-f1cf12c9c469
# ╟─5ddc7b4a-4060-4b12-96c9-3e8ab707b5a7
# ╠═fc7f0733-7d93-474c-9688-73a02c6a51cc
# ╟─d30c9e49-0d28-4882-93fd-9354dce62432
# ╠═5e945e11-869f-4040-af71-c79dd0351f60
# ╠═f5b8520b-0d99-4e8a-8460-32dc1420deb3
# ╠═fb4b04e8-5901-4a7c-bd6d-19c4201fbf01
# ╠═4b6207e2-e835-4c26-a252-734986b31f14
# ╟─9053e773-18fc-43f4-b54b-872a4d22dcb8
# ╠═c30ee767-87cd-44e4-a7a2-a752828c80d4
# ╠═8e681f02-2f7b-49f5-b10b-fd448141dd93
# ╟─0d8d30ab-20dd-4d53-9d8e-ddaa2d1d6182
# ╠═332f7f3e-af9b-4d2b-8e22-551d3cb49213
# ╠═2dd9a2c7-5279-4c5f-9c02-9f006f8fff3e
# ╟─fc26a4e5-1050-4d4d-97cb-2e2b4262d9cd
# ╠═77f54490-4e74-4ba8-866a-9447a21c5b41
# ╠═d20116e4-20e8-4f69-b965-924fd5e69027
# ╠═64961f92-6ab1-49d2-9eb2-53bf92f986df
# ╟─7f7ab7a9-39a9-4954-bbd0-7ad4210c6c1a
# ╠═5fe21528-e4f8-4377-ae03-c0b3dd64ed2f
# ╟─c0453f52-c9d6-4361-be6e-0981f6e62246
# ╠═a45de7a9-0fa3-44f0-abc6-3b6d4380d6ee
# ╟─0cd63574-2afc-43e1-a02f-4ec8eb3d70bc
# ╠═39b549f1-7ef6-41e8-9211-ea1266b347c5
# ╟─78669699-861e-4e3b-b405-88088b303285
# ╟─b869706a-ea91-4e2d-884a-dba0a84334e2
# ╠═b41620cf-1f7f-41d5-b086-25024fab8598

max_stars_repo_name: benharmonics/GeneralPhysics

max_issues_repo_name: benharmonics/GeneralPhysics

max_forks_repo_name: benharmonics/GeneralPhysics


text:
using PyCall

import QuickPOMDPs: preprocess

function preprocess(v::PyObject)
    if pybuiltin("callable")(v)
        return convert(Function, v)
    end

    return v
end

# TODO
# Actually convert the callables
# Maybe we can g

ext: jl

max_stars_repo_name: logankilpatrick/quickpomdps

max_issues_repo_name: logankilpatrick/quickpomdps

max_forks_repo_name: logankilpatrick/quickpomdps


text:
using DynamicPolynomials

@testset "Input-output DynamicPolynomials" begin

    @polyvar x y
    fs = [x, y]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs)
    @test fsfs == [Polynomial(x), Polynomial(y)]

    fs = [x^2, 2y]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs)
    @test fsfs == [Polynomial(x^2), Polynomial(2y)]

    fs = [(3//4)x, (2//9)y, -y]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs)
    @test fsfs == [Polynomial((3//4)x), Polynomial((2//9)y), Polynomial(-y)]

    fs = [x^2*y + 3//4, (2^31 - 5)*x - (2^31 - 4)*y]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs)
    @test fsfs == fs

    fs = [x^2 + 1, (2^31 - 5)*x - (2^31 - 4)*y]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs)
    @test fsfs == fs

    #= DynamicPolynomials does not work with AbstractAlgebra coefficients
    fs = [x^2*y + ff(1), ff(2)*x - (2^31 - 4)*y]
    ring, exps, cfs = Groebner.convert_to_internal(fs, :input)
    fsfs = Groebner.convert_to_output(ring, fs, exps, cfs)
    @test fsfs == fs
    =#

end

@testset "DynamicPolynomials io consistency" begin
    @polyvar x y

    fs = [x, y]
    @test Groebner.groebner(fs) isa Vector{Polynomial{true, Int64}}

    fs = [UInt16(2)x + UInt16(3), UInt16(2)y]
    @test Groebner.groebner(fs) isa Vector{Polynomial{true, UInt16}}

    fs = [BigInt(1)x + BigInt(20), y]
    @test Groebner.groebner(fs) isa Vector{Polynomial{true, BigInt}}

    fs = [x + BigInt(20)//BigInt(3), y]
    @test Groebner.groebner(fs) isa Vector{Polynomial{true, Rational{BigInt}}}

    fs = [34343343433x*y^2 + 3431234567833, 3434343434x*y - 342343242342]
    @test_throws DomainError Groebner.groebner(fs)

    fs = [BigInt(34343343433)x*y^2 + 3431234567833, 3434343434x*y - 342343242342]
    @test Groebner.groebner(fs) isa Vector{Polynomial{true, BigInt}}
end


text:
using Transparency

"""
    function B_ν(ν, T)

Planck's law! Radiation in LTE. Takes frequency and temperature, returns
specific intensity
"""
function B_ν(ν, T)::AbstractFloat
    return 2*h*ν^3/c_0^2 * 1/(exp(h*ν/(k_B*T)) - 1)
end

"""
    function B_ν(λ, T)

Planck's law! Radiation in LTE. Takes wavelength and temperature, returns
specific intensity
"""
function B_λ(λ::Unitful.Length, T::Unitful.Temperature)
    return 2*h*c_0^2/λ^5 * 1/(exp(h*c_0/(λ*k_B*T)) - 1)
end

"""
    function α_absorption(λ::Unitful.Length, temperature::Unitful.Temperature,
               electron_density::NumberDensity, h_ground_density::NumberDensity,
               proton_density::NumberDensity)

Extinction from photon destruction processes.
"""
function α_absorption(λ::Unitful.Length, temperature::Unitful.Temperature,
                      electron_density::NumberDensity, h_neutral_density::NumberDensity,
                      proton_density::NumberDensity)

    # α = max(0u"m^-1", Transparency.hminus_ff_stilley(λ, temperature, h_ground_density, electron_density))
    α = hminus_ff(λ, temperature, h_neutral_density, electron_density; recipe="stilley")
    α += hminus_bf(λ, temperature, h_neutral_density, electron_density; recipe="geltman")
    α += hydrogenic_ff(c_0 / λ, temperature, electron_density, proton_density, 1)
    # hydrogenic_bf ??? why return return 0 * αbf_const * species_densitys
    α += h2plus_ff(λ, temperature, h_neutral_density, proton_density)
    α += h2plus_bf(λ, temperature, h_neutral_density, proton_density)
    return α
end

"""
    α_scattering(λ::Unitful.Length,
                 electron_density::NumberDensity,
                 h_ground_density::NumberDensity)

Extinction from scattering processes.
"""
function α_scattering(λ::Unitful.Length,
                      electron_density::NumberDensity,
                      h_ground_density::NumberDensity)

   α = thomson(electron_density)
   α += rayleigh_h(λ, h_ground_density) # Gives the edge... (1217.7 Å)
   return α
end

# function α_absorption(λ::Unitful.Length, temperature::Unitful.Temperature,
                      # electron_density::NumberDensity, h_ground_density::NumberDensity,
                      # proton_density::NumberDensity)
#
    # α = max(0u"m^-1", Transparency.hminus_ff_stilley(λ, temperature, h_ground_density, electron_density))
    # α += Transparency.hminus_bf_wbr(λ, temperature, h_ground_density, electron_density)
    # α += hydrogenic_ff(c_0 / λ, temperature, electron_density, proton_density, 1)
    # α += h2plus_ff(λ, temperature, h_ground_density, proton_density)
    # α += h2plus_bf(λ, temperature, h_ground_density, proton_density)
    # return α
# end

max_stars_repo_name: meudnaes/VoronoiRT

max_issues_repo_name: meudnaes/VoronoiRT

max_forks_repo_name: meudnaes/VoronoiRT


text:
"""
    x = ComponentArray(nt::NamedTuple)
    x = ComponentArray(;kwargs...)
    x = ComponentArray(data::AbstractVector, ax)
    x = ComponentArray{T}(args...; kwargs...) where T

Array type that can be accessed like an arbitrary nested mutable struct.

# Examples

```jldoctest
julia> using ComponentArrays

julia> x = ComponentArray(a=1, b=[2, 1, 4], c=(a=2, b=[1, 2]))
ComponentVector{Int64}(a = 1, b = [2, 1, 4], c = (a = 2, b = [1, 2]))

julia> x.c.a = 400; x
ComponentVector{Int64}(a = 1, b = [2, 1, 4], c = (a = 400, b = [1, 2]))

julia> x[5]
400

julia> collect(x)
7-element Array{Int64,1}:
   1
   2
   1
   4
 400
   1
   2
```
"""
struct ComponentArray{T,N,A<:AbstractArray{T,N},Axes<:Tuple{Vararg{AbstractAxis}}} <: DenseArray{T,N}
    data::A
    axes::Axes
end

# Entry from type (used for broadcasting)
ComponentArray{Axes}(data) where Axes = ComponentArray(data, getaxes(Axes)...)
ComponentArray(::UndefInitializer, ax::Axes) where Axes<:Tuple =
    ComponentArray(similar(Array{Float64}, last_index.(ax)), ax...)
ComponentArray{A}(::UndefInitializer, ax::Axes) where {A<:AbstractArray,Axes<:Tuple} =
    ComponentArray(similar(A, last_index.(ax)), ax...)
ComponentArray{T}(::UndefInitializer, ax::Axes) where {T,Axes<:Tuple} =
    ComponentArray(similar(Array{T}, last_index.(ax)), ax...)

# Entry from data array and AbstractAxis types dispatches to correct shapes and partitions
# then packs up axes into a tuple for inner constructor
ComponentArray(data, ::FlatAxis...) = data
ComponentArray(data, ax::NotShapedOrPartitionedAxis...) = ComponentArray(data, ax)
ComponentArray(data, ax::NotPartitionedAxis...) = ComponentArray(maybe_reshape(data, ax...), unshape.(ax)...)
function ComponentArray(data, ax::AbstractAxis...)
    part_axs = filter_by_type(PartitionedAxis, ax...)
    part_data = partition(data, size.(part_axs)...)
    axs = Axis.(ax)
    return LazyArray(ComponentArray(x, axs) for x in part_data)
end

# Entry from NamedTuple, Dict, or kwargs
ComponentArray{T}(nt::NamedTuple) where T = ComponentArray(make_carray_args(T, nt)...)
ComponentArray(nt::NamedTuple) = ComponentArray(make_carray_args(nt)...)
ComponentArray(d::AbstractDict) = ComponentArray(NamedTuple{Tuple(keys(d))}(values(d)))
ComponentArray{T}(;kwargs...) where T = ComponentArray{T}((;kwargs...))
ComponentArray(;kwargs...) = ComponentArray((;kwargs...))

ComponentArray(x::ComponentArray) = x
ComponentArray{T}(x::ComponentArray) where {T} = T.(x)
(CA::Type{<:ComponentArray{T,N,A,Ax}})(x::ComponentArray) where {T,N,A,Ax} = ComponentArray(T.(getdata(x)), getaxes(x))

## Some aliases
"""
    x = ComponentVector(nt::NamedTuple)
    x = ComponentVector(;kwargs...)
    x = ComponentVector(data::AbstractVector, ax)
    x = ComponentVector{T}(args...; kwargs...) where T

A `ComponentVector` is an alias for a one-dimensional `ComponentArray`.
"""
const ComponentVector{T,A,Axes} = ComponentArray{T,1,A,Axes}
ComponentVector(nt) = ComponentArray(nt)
ComponentVector{T}(nt) where {T} = ComponentArray{T}(nt)
ComponentVector(;kwargs...) = ComponentArray(;kwargs...)
ComponentVector{T}(;kwargs...) where {T} = ComponentArray{T}(;kwargs...)
ComponentVector{T}(::UndefInitializer, ax) where {T} = ComponentArray{T}(undef, ax)
ComponentVector(data::AbstractVector, ax) = ComponentArray(data, ax)
ComponentVector(data::AbstractArray, ax) = throw(DimensionMismatch("A `ComponentVector` must be initialized with a 1-dimensional array. This array is $(ndims(data))-dimensional."))

# Add new fields to component Vector
function ComponentArray(x::ComponentVector; kwargs...)
    return foldl((x1, kwarg) -> _maybe_add_field(x1, kwarg), (kwargs...,); init=x)
end
ComponentVector(x::ComponentVector; kwargs...) = ComponentArray(x; kwargs...)

ComponentVector{T}(x::ComponentVector) where {T} = T.(x)

"""
    x = ComponentMatrix(data::AbstractMatrix, ax...)
    x = ComponentMatrix{T}(data::AbstractMatrix, ax...) where T

A `ComponentMatrix` is an alias for a two-dimensional `ComponentArray`.
"""
const ComponentMatrix{T,A,Axes} = ComponentArray{T,2,A,Axes}
ComponentMatrix{T}(::UndefInitializer, ax...) where {T} = ComponentArray{T}(undef, ax...)
ComponentMatrix(data::AbstractMatrix, ax...) = ComponentArray(data, ax...)
ComponentMatrix(data::AbstractArray, ax...) = throw(DimensionMismatch("A `ComponentMatrix` must be initialized with a 2-dimensional array. This array is $(ndims(data))-dimensional."))

ComponentMatrix(x::ComponentMatrix) = x
ComponentMatrix{T}(x::ComponentMatrix) where {T} = T.(x)

const CArray = ComponentArray
const CVector = ComponentVector
const CMatrix = ComponentMatrix

const AdjOrTrans{T, A} = Union{Adjoint{T, A}, Transpose{T, A}}
const AdjOrTransComponentArray{T, A} = Union{Adjoint{T, A}, Transpose{T, A}} where A<:ComponentArray

## Constructor helpers
# For making ComponentArrays from named tuples
function make_carray_args(nt)
    data, ax = make_carray_args(Vector, nt)
    data = length(data)==1 ? [data[1]] : reduce(vcat, data)
    return (data, ax)
end
make_carray_args(T::Type, nt) = make_carray_args(Vector{T}, nt)
function make_carray_args(A::Type{<:AbstractArray}, nt)
    data, idx = make_idx([], nt, 0)
    return (A(data), Axis(idx))
end

# Builds up data vector and returns appropriate AbstractAxis type for each input type
function make_idx(data, nt::NamedTuple, last_val)
    len = recursive_length(nt)
    kvs = []
    lv = 0
    for (k,v) in zip(keys(nt), nt)
        (_,val) = make_idx(data, v, lv)
        push!(kvs, k => val)
        lv = val
    end
    return (data, ViewAxis(last_index(last_val) .+ (1:len), (;kvs...)))
end
function make_idx(data, pair::Pair, last_val)
    data, ax = make_idx(data, pair.second, last_val)
    return (data, ViewAxis(last_val:(last_val+len-1), Axis(pair.second)))
end
make_idx(data, x, last_val) = (
    push!(data, x),
    ViewAxis(last_index(last_val) + 1)
)
make_idx(data, x::ComponentVector, last_val) = (
    pushcat!(data, x),
    ViewAxis(
        last_index(last_val) .+ (1:length(x)),
        getaxes(x)[1]
    )
)
function make_idx(data, x::AbstractArray, last_val)
    pushcat!(data, x)
    out = last_index(last_val) .+ (1:length(x))
    return (data, ViewAxis(out, ShapedAxis(size(x))))
end
function make_idx(data, x::A, last_val) where A<:AbstractArray{<:Union{NamedTuple, ComponentArray}}
    len = recursive_length(x)
    if eltype(x) |> isconcretetype
        out = ()
        for elem in x
            (_,out) = make_idx(data, elem, last_val)
        end
        return (
            data,
            ViewAxis(
                last_index(last_val) .+ (1:len),
                PartitionedAxis(
                    len ÷ length(x),
                    indexmap(out)
                )
            )
        )
    else
        error("Only homogeneous arrays of inner ComponentArrays are allowed.")
    end
end

#TODO: Make all internal function names start with underscores
_maybe_add_field(x, pair) = haskey(x, pair.first) ? _update_field(x, pair) : _add_field(x, pair)
function _add_field(x, pair)
    data = copy(getdata(x))
    new_data, new_ax = make_idx(data, pair.second, length(data))
    new_ax = Axis(NamedTuple{tuple(pair.first)}(tuple(new_ax)))
    new_ax = merge(getaxes(x)[1], new_ax)
    return ComponentArray(new_data, new_ax)
end
function _update_field(x, pair)
    x_copy = copy(x)
    x_copy[pair.first] = pair.second
    return x_copy
end

pushcat!(a, b) = reduce((x1,x2) -> push!(x1,x2), b; init=a)

# Reshape ComponentArrays with ShapedAxis axes
maybe_reshape(data, ::NotShapedOrPartitionedAxis...) = data
function maybe_reshape(data, axs::AbstractAxis...)
    shapes = filter_by_type(ShapedAxis, axs...) .|> size
    shapes = reduce((tup, s) -> (tup..., s...), shapes)
    return reshape(data, shapes)
end

# Recurse through nested ViewAxis types to find the last index
last_index(x) = last(x)
last_index(x::ViewAxis) = last_index(viewindex(x))
last_index(x::AbstractAxis) = last_index(last(indexmap(x)))

# Reduce singleton dimensions
remove_nulls() = ()
remove_nulls(x1, args...) = (x1, remove_nulls(args...)...)
remove_nulls(::NullAxis, args...) = (remove_nulls(args...)...,)

## Attributes
"""
    getdata(x::ComponentArray)

Access ```.data``` field of a ```ComponentArray```, which contains the array that ```ComponentArray``` wraps.
"""
@inline getdata(x::ComponentArray) = getfield(x, :data)
@inline getdata(x) = x
@inline getdata(x::Adjoint) = getdata(x.parent)'
@inline getdata(x::Transpose) = transpose(getdata(x.parent))

"""
    getaxes(x::ComponentArray)

Access ```.axes``` field of a ```ComponentArray```. This is different than ```axes(x::ComponentArray)```, which
    returns the axes of the contained array.

# Examples

```jldoctest
julia> using ComponentArrays

julia> ax = Axis(a=1:3, b=(4:6, (a=1, b=2:3)))
Axis{(a = 1:3, b = (4:6, (a = 1, b = 2:3)))}()

julia> A = zeros(6,6);

julia> ca = ComponentArray(A, (ax, ax))
6×6 ComponentMatrix{Float64} with axes Axis(a = 1:3, b = (4:6, (a = 1, b = 2:3))) × Axis(a = 1:3, b = (4:6, (a = 1, b = 2:3)))
 0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0

julia> getaxes(ca)
(Axis(a = 1:3, b = (4:6, (a = 1, b = 2:3))), Axis(a = 1:3, b = (4:6, (a = 1, b = 2:3))))
```
"""
@inline getaxes(x::ComponentArray) = getfield(x, :axes)
@inline getaxes(x::AdjOrTrans{T, <:ComponentVector}) where T = (FlatAxis(), getaxes(x.parent)[1])
@inline getaxes(x::AdjOrTrans{T, <:ComponentMatrix}) where T = reverse(getaxes(x.parent))

@inline getaxes(::Type{<:ComponentArray{T,N,A,<:Axes}}) where {T,N,A,Axes} = map(x->x(), (Axes.types...,))
@inline getaxes(::Type{<:AdjOrTrans{T,<:CA}}) where {T,CA<:ComponentVector} = (FlatAxis(), getaxes(CA)[1]) |> typeof
@inline getaxes(::Type{<:AdjOrTrans{T,<:CA}}) where {T,CA<:ComponentMatrix} = reverse(getaxes(CA)) |> typeof

Base.parent(x::ComponentArray) = getfield(x, :data)

Base.size(x::ComponentArray) = size(getdata(x))

Base.elsize(x::Type{<:ComponentArray{T,N,A,Axes}}) where {T,N,A,Axes} = Base.elsize(A)

Base.axes(x::ComponentArray) = axes(getdata(x))

Base.reinterpret(::Type{T}, x::ComponentArray, args...) where T = ComponentArray(reinterpret(T, getdata(x), args...), getaxes(x))

Base.propertynames(x::ComponentVector) = propertynames(indexmap(getaxes(x)[1]))

Base.keys(x::ComponentVector) = keys(indexmap(getaxes(x)[1]))

"""
    valkeys(x::ComponentVector)

Returns `Val`-wrapped keys of `ComponentVector` for fast iteration over component keys.

# Examples

```julia-repl
julia> using ComponentArrays

julia> ca = ComponentArray(a=1, b=[1,2,3], c=(a=4,))
ComponentVector{Int64}(a = 1, b = [1, 2, 3], c = (a = 4))

julia> [ca[k] for k in valkeys(ca)]
3-element Array{Any,1}:
 1
  [1, 2, 3]
  ComponentVector{Int64,SubArray...}(a = 4)

julia> sum(prod(ca[k]) for k in valkeys(ca))
11
```
"""
@generated function valkeys(ca::ComponentVector)
    idxmap = ComponentArrays.indexmap(getaxes(ca)[1])
    k = Val.(keys(idxmap))
    return :($k)
end

Base.haskey(x::ComponentVector, s::Symbol) = haskey(indexmap(getaxes(x)[1]), s)

function Base.permutedims(x::ComponentArray, dims)
    axs = getaxes(x)
    return ComponentArray(permutedims(getdata(x), dims), map(i->axs[i], dims)...)
end

Base.IndexStyle(::Type{<:ComponentArray{T,N,<:A,<:Axes}}) where {T,N,A,Axes} = IndexStyle(A)

ArrayInterface.parent_type(::Type{ComponentArray{T,N,A,Axes}}) where {T,N,A,Axes} = A
ArrayInterface.size(A::ComponentArray) = ArrayInterface.size(parent(A))
ArrayInterface.strides(A::ComponentArray) = ArrayInterface.strides(parent(A))
for f in [:device, :stride_rank, :contiguous_axis, :contiguous_batch_size, :dense_dims]
    @eval ArrayInterface.$f(::Type{ComponentArray{T,N,A,Axes}}) where {T,N,A,Axes} = ArrayInterface.$f(A)
end


text:
using DataDeps
using MAT

function register_burgers()
    register(DataDep(
        "Burgers",
        """
        Burgers' equation dataset from
        [fourier_neural_operator](https://github.com/zongyi-li/fourier_neural_operator)
        """,
        "http://www.med.cgu.edu.tw/NeuralOperators/Burgers_R10.zip",
        "9cbbe5070556c777b1ba3bacd49da5c36ea8ed138ba51b6ee76a24b971066ecd",
        post_fetch_method=unpack
    ))
end

function get_data(; n=2048, Δsamples=2^3, grid_size=div(2^13, Δsamples), T=Float32)
    file = matopen(joinpath(datadep"Burgers", "burgers_data_R10.mat"))
    x_data = T.(collect(read(file, "a")[1:n, 1:Δsamples:end]'))
    y_data = T.(collect(read(file, "u")[1:n, 1:Δsamples:end]'))
    close(file)

    x_loc_data = Array{T, 3}(undef, 2, grid_size, n)
    x_loc_data[1, :, :] .= reshape(repeat(LinRange(0, 1, grid_size), n), (grid_size, n))
    x_loc_data[2, :, :] .= x_data

    return x_loc_data, y_data
end

function get_data_don(; n=2048, Δsamples=2^3, grid_size=div(2^13, Δsamples))
    file = matopen(joinpath(datadep"Burgers", "burgers_data_R10.mat"))
    x_data = collect(read(file, "a")[1:n, 1:Δsamples:end])
    y_data = collect(read(file, "u")[1:n, 1:Δsamples:end])
    close(file)

    return x_data, y_data
end

function get_dataloader(; n_train=1800, n_test=200, batchsize=100)
    𝐱, 𝐲 = get_data(n=2048)

    𝐱_train, 𝐲_train = 𝐱[:, :, 1:n_train], 𝐲[:, 1:n_train]
    loader_train = Flux.DataLoader((𝐱_train, 𝐲_train), batchsize=batchsize, shuffle=true)

    𝐱_test, 𝐲_test = 𝐱[:, :, end-n_test+1:end], 𝐲[:, end-n_test+1:end]
    loader_test = Flux.DataLoader((𝐱_test, 𝐲_test), batchsize=batchsize, shuffle=false)

    return loader_train, loader_test
end


text:
using LinearAlgebra
using Printf

function optimal_design_fw_aa(
    items,
    measures,
    max_iter = 50000,
    early_stopping = 1e-5,
    drop = false,
    fix_threshold = 1e-5,
)
    K = length(items[1, :])
    L = length(measures[:, 1])
    design = ones(length(items[:, 1]))
    design /= sum(design)
    Vinv = Matrix{Float64}(I, K, K)

    for count = 1:max_iter
        idx = randmin([maximum([transpose(measures[i, :]) *
                                sherman_morrison(Vinv, items[j, :]) * measures[i, :] for i = 1:L]) for j = 1:K])

        Vinv = sherman_morrison(Vinv, items[idx, :])

        γ = 1 / (count + 1)
        design_update = -γ * design
        design_update[idx] += γ

        relative = norm(design_update) / norm(design)

        design += design_update

        if relative < early_stopping
            println("early stopped")
            break
        end
    end

    if drop
        drop_total = sum(filter(x -> x < fix_threshold, design))
        if drop_total > 0
            design = replace(x -> x < fix_threshold ? 0 : x, design)
            design[argmax(design)] += drop_total
        end
    end

    return design
end

function optimal_design_fw_ab(
    items,
    measures,
    μ,
    max_iter = 10000,
    early_stopping = 1e-5,
    drop = false,
    fix_threshold = 1e-5,
)
    K = length(items[1, :])
    L = length(measures[:, 1])
    design = ones(length(items[:, 1]))
    design /= sum(design)
    measure_counts = ones(L)
    Vinv = Matrix{Float64}(I, K, K)

    for count = 1:max_iter
        ida = argmin([sum([measure_counts[i] * (-2) *
                           (transpose(items[j, :]) * Vinv * measures[i, :]) .^ 2 for i = 1:L]) for j = 1:K])
        # r = [measure_counts[i] * (-2) * transpose(items[1]) * Vinv * measures[i, :] for i = 1:L]
        # @show r
        idb = argmax([transpose(measures[i, :]) * Vinv * measures[i, :] for i = 1:L])
        # s = [transpose(measures[i, :]) * Vinv * measures[i, :] for i = 1:L]
        # @show r, idb, ida

        Vinv = sherman_morrison(Vinv, items[ida, :])
        measure_counts[idb] += 1

        γ = 1 / (count + 1)
        design_update = -γ * design
        design_update[ida] += γ

        relative = norm(design_update) / norm(design)

        design += design_update

        if relative < early_stopping
            println("early stopped")
            break
        end
    end

    if drop
        drop_total = sum(filter(x -> x < fix_threshold, design))
        if drop_total > 0
            design = replace(x -> x < fix_threshold ? 0 : x, design)
            design[argmax(design)] += drop_total
        end
    end

    # compute the complexity
    AB = maximum([transpose(measures[i, :]) * Vinv * measures[i, :] for i = 1:L])
    value = maximum([2 * transpose(items[1, :] - items[i, :]) * Vinv *
                     (items[1, :] - items[i, :]) / (items[1, :] - items[i, :])'μ for i = 2:K])

    return design, AB, value
end

function optimal_design_fw(
    items,
    measures,
    μ,
    max_iter = 10000,
    early_stopping = 1e-5,
    drop = false,
    fix_threshold = 1e-5,
)
    K = length(items[1])
    L = length(measures)
    design = ones(length(items))
    design /= sum(design)
    measure_counts = ones(L)
    Ainv = Matrix{Float64}(I, K, K)

    for count = 1:max_iter
        ida = argmin([sum([measure_counts[i] * (-2) *
                           (transpose(items[j]) * Ainv * measures[i]) .^ 2 for i = 1:L]) for j = 1:K])
        idb = argmax([transpose(measures[i]) * Ainv * measures[i] for i = 1:L])

        Ainv = sherman_morrison(Ainv, items[ida])
        measure_counts[idb] += 1

        γ = 1 / (count + 1)
        design_update = -γ * design
        design_update[ida] += γ

        relative = norm(design_update) / norm(design)

        design += design_update

        if relative < early_stopping
            #println("early stopped")
            break
        end
    end

    if drop
        drop_total = sum(filter(x -> x < fix_threshold, design))
        if drop_total > 0
            design = replace(x -> x < fix_threshold ? 0 : x, design)
            design[argmax(design)] += drop_total
        end
    end

    # compute the complexity
    AB = maximum([transpose(measures[i]) * Ainv * measures[i] for i = 1:L])
    value = maximum([2 * transpose(items[1] - items[i]) * Ainv *
                     (items[1] - items[i]) / (items[1] - items[i])'μ for i = 2:K])

    return design, AB, value
end

function build_XY_ab(items)
    num_items = size(items)[1]
    num_features = size(items)[2]
    Y = zeros((num_items * num_items, num_features))

    for i = 1:num_items
        Y[(num_items*(i-1)+1):num_items*i, :] = items .- items[i, :]'
    end

    return Y
end

function build_XY(items)
    num_items = length(pep.arms)
    num_features = length(pep.arms[1])

    Y = Vector{Float64}[]
    for i = 1:num_items
        for j = 1:num_items
            v = items[j] - items[i]
            push!(Y, v)
        end
    end

    return Y
end

function build_T_ab(items, xstar, θ)
    num_items = size(items)[1]
    num_features = size(items)[2]
    Y = zeros((num_items - 1, num_features))

    for i = 1:(num_items-1)
        Y[i, :] = sqrt(2) * (xstar - items[i+1, :]) / (xstar - items[i+1, :])'θ
    end

    return Y
end

function build_T(items, xstar, θ)
    num_items = length(items)
    num_features = length(items[1])

    Y = Vector{Float64}[]
    for i = 1:num_items
        if xstar == items[i]
            continue
        else
            v = sqrt(2) * (xstar - items[i]) / (xstar - items[i])'θ
            push!(Y, v)
        end
    end

    return Y
end

"""
Tests
"""
# dist = Gaussian();
# dim = 4
# μ = zeros(dim);
# µ[1] = 0.5; μ[2] = 0.45; μ[3] = 0.43; μ[4] = 0.4
#
# arms = Vector{Float64}[]
# for k = 1:dim
#     v = zeros(dim)
#     v[k] = 1.0
#     push!(arms, v)
# end
# # ω = pi / 6
# # v = zeros(dim);
# # v[1] = cos(ω);
# # v[2] = sin(ω);
# # push!(arms, v)
#
# pep = LinearBestArm(dist, arms);
# designAA, _, _ = optimal_design_fw(pep.arms, pep.arms, μ)
# designABdir, _, _ = optimal_design_fw(pep.arms, build_XY(pep.arms), μ)
# designABstar, _, _ = optimal_design_fw(pep.arms, build_T(pep.arms, pep.arms[1], μ), μ, 50000)
# @show designABstar


text:
#sgu.jl
#Code to implement Schmitt-Grohe and Uribe (2003, JIE) Model 2
#Shamelessly copied from Dynare .mod file of Cesa-Bianchi (2012)
#https://sites.google.com/site/ambropo/dynarecodes
#and from Dynare .mod file of Pfeifer (2019)
#https://github.com/JohannesPfeifer/DSGE_mod/tree/master/SGU_2003

function sgu()
    ∞ = Inf
    # Parameters
    @variables γ, ω, ρ, σe, δ, ψ, α, ϕ, β, r_w, d_bar
    # x and y
    @variables t::Integer, d(..), c(..), h(..), GDP(..), i(..), k(..), a(..), λ(..), tb(..),
               ca(..), riskpremium(..), r(..), Ld(..), Lr(..), Lk(..)

    x = [a, Ld, Lk, Lr]
    y = [d, c, h, GDP, i, k, λ, tb, ca, riskpremium, r]
    p = [γ, ω, ρ, σe, δ, ψ, α, ϕ, β, r_w, d_bar]

    # Model equations
    # Rewritten from Dynare to Schmitt-Grohe and Uribe timing conventions
    H = [d(t) - (1 + exp(Lr(t))) * Ld(t) + exp(GDP(t)) - exp(c(t)) - exp(i(t)) -
         (ϕ / 2) * (exp(k(t)) - exp(Lk(t)))^2, #Debt evolution
         exp(GDP(t)) - exp(a(t)) * (exp(Lk(t))^α) * (exp(h(t))^(1 - α)), #Production function
         exp(k(t)) - exp(i(t)) - (1 - δ) * exp(Lk(t)), #Capital evolution
         exp(λ(t)) - β * (1 + exp(r(t))) * exp(λ(t + 1)), #Euler equation
         (exp(c(t)) - ((exp(h(t))^ω) / ω))^(-γ) - exp(λ(t)), #Marginal utility
         ((exp(c(t)) - ((exp(h(t))^ω) / ω))^(-γ)) * (exp(h(t))^(ω - 1)) -
         exp(λ(t)) * (1 - α) * exp(GDP(t)) / exp(h(t)), #Labor FOC
         exp(λ(t)) * (1 + ϕ * (exp(k(t)) - exp(Lk(t)))) -
         β *
         exp(λ(t + 1)) *
         (α * exp(GDP(t + 1)) / exp(k(t)) + 1 - δ + ϕ * (exp(k(t + 1)) - exp(k(t)))), #Investment FOC
         a(t + 1) - ρ * a(t), #TFP evolution
         exp(r(t)) - r_w - riskpremium(t), #Interest rate
         riskpremium(t) - ψ * (exp(d(t) - d_bar) - 1), #Risk premium
         tb(t) - 1 +
         ((exp(c(t)) + exp(i(t)) + (ϕ / 2) * (exp(k(t)) - exp(Lk(t)))^2) / exp(GDP(t))), #Trade Balance
         ca(t) - (1 / exp(GDP(t))) * (Ld(t) - d(t)), #Current Account
         Ld(t + 1) - d(t), #Auxiliary lagged variables
         Lk(t + 1) - k(t), Lr(t + 1) - r(t)]

    #Define representations for steady state

    hstar = ((1 - α) * (α / (r_w + δ))^(α / (1 - α)))^(1 / (ω - 1))
    kstar = hstar / (((r_w + δ) / α)^(1 / (1 - α)))
    istar = δ * kstar
    GDPstar = (kstar^α) * (hstar^(1 - α))
    cstar = GDPstar - istar - r_w * d_bar
    tbstar = GDPstar - cstar - istar
    λstar = (cstar - ((hstar^ω) / ω))^(-γ)

    #Steady state values

    steady_states = [a(∞) ~ 0, Ld(∞) ~ d_bar, Lk(∞) ~ log(kstar), Lr(∞) ~ log((1 - β) / β),
                     d(∞) ~ d_bar, c(∞) ~ log(cstar), h(∞) ~ log(hstar),
                     GDP(∞) ~ log(GDPstar), i(∞) ~ log(istar), k(∞) ~ log(kstar),
                     λ(∞) ~ log(λstar), tb(∞) ~ log(tbstar), ca(∞) ~ 0, riskpremium(∞) ~ 0,
                     r(∞) ~ log((1 - β) / β)]
    n_ϵ = 1
    Γ = reshape([σe], n_ϵ, n_ϵ)
    η = reshape([-1, 0, 0, 0], length(x), n_ϵ) # η is n_x * n_ϵ matrix

    # Add some observation matrix for Q and Ω?
    return H, (; t, x, y, p, steady_states, Γ, η), "sgu"
end


text:
using Oceananigans.BoundaryConditions: PBC, ZFBC, NFBC

function instantiate_boundary_function(B, X1, X2, func)
    boundary_function = BoundaryFunction{B, X1, X2}(func)
    return true
end

function instantiate_tracer_boundary_condition(bctype, B, func)
    boundary_condition = TracerBoundaryCondition(bctype, B, func)
    return true
end

function instantiate_u_boundary_condition(bctype, B, func)
    boundary_condition = UVelocityBoundaryCondition(bctype, B, func)
    return true
end

function instantiate_v_boundary_condition(bctype, B, func)
    boundary_condition = VVelocityBoundaryCondition(bctype, B, func)
    return true
end

function instantiate_w_boundary_condition(bctype, B, func)
    boundary_condition = WVelocityBoundaryCondition(bctype, B, func)
    return true
end

@testset "Boundary conditions" begin
    @info "Testing boundary conditions..."

    @testset "Boundary functions" begin
        @info "  Testing boundary functions..."

        simple_bc(ξ, η, t) = exp(ξ) * cos(η) * sin(t)
        for B in (:x, :y, :z)
            for X1 in (:Face, :Cell)
                @test instantiate_boundary_function(B, X1, Cell, simple_bc)
            end

            for bctype in (Value, Gradient)
                @test instantiate_tracer_boundary_condition(bctype, B, simple_bc)
                @test instantiate_u_boundary_condition(bctype, B, simple_bc)
                @test instantiate_v_boundary_condition(bctype, B, simple_bc)
                @test instantiate_w_boundary_condition(bctype, B, simple_bc)
            end
        end
    end

    @testset "Field boundary conditions" begin
        @info "  Testing field boundary functions..."

        # Triply periodic
        ppp_topology = (Periodic, Periodic, Periodic)
        ppp_grid = RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1), topology=ppp_topology)

        u_bcs = UVelocityBoundaryConditions(ppp_grid)
        v_bcs = VVelocityBoundaryConditions(ppp_grid)
        w_bcs = WVelocityBoundaryConditions(ppp_grid)
        T_bcs = TracerBoundaryConditions(ppp_grid)

        @test u_bcs isa FieldBoundaryConditions
        @test u_bcs.x.left  isa PBC
        @test u_bcs.x.right isa PBC
        @test u_bcs.y.left  isa PBC
        @test u_bcs.y.right isa PBC
        @test u_bcs.z.left  isa PBC
        @test u_bcs.z.right isa PBC

        @test v_bcs isa FieldBoundaryConditions
        @test v_bcs.x.left  isa PBC
        @test v_bcs.x.right isa PBC
        @test v_bcs.y.left  isa PBC
        @test v_bcs.y.right isa PBC
        @test v_bcs.z.left  isa PBC
        @test v_bcs.z.right isa PBC

        @test w_bcs isa FieldBoundaryConditions
        @test w_bcs.x.left  isa PBC
        @test w_bcs.x.right isa PBC
        @test w_bcs.y.left  isa PBC
        @test w_bcs.y.right isa PBC
        @test w_bcs.z.left  isa PBC
        @test w_bcs.z.right isa PBC

        @test T_bcs isa FieldBoundaryConditions
        @test T_bcs.x.left  isa PBC
        @test T_bcs.x.right isa PBC
        @test T_bcs.y.left  isa PBC
        @test T_bcs.y.right isa PBC
        @test T_bcs.z.left  isa PBC
        @test T_bcs.z.right isa PBC

        # Doubly periodic. Engineers call this a "Channel geometry".
        ppb_topology = (Periodic, Periodic, Bounded)
        ppb_grid = RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1), topology=ppb_topology)

        u_bcs = UVelocityBoundaryConditions(ppb_grid)
        v_bcs = VVelocityBoundaryConditions(ppb_grid)
        w_bcs = WVelocityBoundaryConditions(ppb_grid)
        T_bcs = TracerBoundaryConditions(ppb_grid)

        @test u_bcs isa FieldBoundaryConditions
        @test u_bcs.x.left  isa PBC
        @test u_bcs.x.right isa PBC
        @test u_bcs.y.left  isa PBC
        @test u_bcs.y.right isa PBC
        @test u_bcs.z.left  isa ZFBC
        @test u_bcs.z.right isa ZFBC

        @test v_bcs isa FieldBoundaryConditions
        @test v_bcs.x.left  isa PBC
        @test v_bcs.x.right isa PBC
        @test v_bcs.y.left  isa PBC
        @test v_bcs.y.right isa PBC
        @test v_bcs.z.left  isa ZFBC
        @test v_bcs.z.right isa ZFBC

        @test w_bcs isa FieldBoundaryConditions
        @test w_bcs.x.left  isa PBC
        @test w_bcs.x.right isa PBC
        @test w_bcs.y.left  isa PBC
        @test w_bcs.y.right isa PBC
        @test w_bcs.z.left  isa NFBC
        @test w_bcs.z.right isa NFBC

        @test T_bcs isa FieldBoundaryConditions
        @test T_bcs.x.left  isa PBC
        @test T_bcs.x.right isa PBC
        @test T_bcs.y.left  isa PBC
        @test T_bcs.y.right isa PBC
        @test T_bcs.z.left  isa ZFBC
        @test T_bcs.z.right isa ZFBC

        # Singly periodic. Oceanographers call this a "Channel", engineers call it a "Pipe"
        pbb_topology = (Periodic, Bounded, Bounded)
        pbb_grid = RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1), topology=pbb_topology)

        u_bcs = UVelocityBoundaryConditions(pbb_grid)
        v_bcs = VVelocityBoundaryConditions(pbb_grid)
        w_bcs = WVelocityBoundaryConditions(pbb_grid)
        T_bcs = TracerBoundaryConditions(pbb_grid)

        @test u_bcs isa FieldBoundaryConditions
        @test u_bcs.x.left  isa PBC
        @test u_bcs.x.right isa PBC
        @test u_bcs.y.left  isa ZFBC
        @test u_bcs.y.right isa ZFBC
        @test u_bcs.z.left  isa ZFBC
        @test u_bcs.z.right isa ZFBC

        @test v_bcs isa FieldBoundaryConditions
        @test v_bcs.x.left  isa PBC
        @test v_bcs.x.right isa PBC
        @test v_bcs.y.left  isa NFBC
        @test v_bcs.y.right isa NFBC
        @test v_bcs.z.left  isa ZFBC
        @test v_bcs.z.right isa ZFBC

        @test w_bcs isa FieldBoundaryConditions
        @test w_bcs.x.left  isa PBC
        @test w_bcs.x.right isa PBC
        @test w_bcs.y.left  isa ZFBC
        @test w_bcs.y.right isa ZFBC
        @test w_bcs.z.left  isa NFBC
        @test w_bcs.z.right isa NFBC

        @test T_bcs isa FieldBoundaryConditions
        @test T_bcs.x.left  isa PBC
        @test T_bcs.x.right isa PBC
        @test T_bcs.y.left  isa ZFBC
        @test T_bcs.y.right isa ZFBC
        @test T_bcs.z.left  isa ZFBC
        @test T_bcs.z.right isa ZFBC

        # Triply bounded. Oceanographers call this a "Basin", engineers call it a "Box"
        bbb_topology = (Bounded, Bounded, Bounded)
        bbb_grid = RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1), topology=bbb_topology)

        u_bcs = UVelocityBoundaryConditions(bbb_grid)
        v_bcs = VVelocityBoundaryConditions(bbb_grid)
        w_bcs = WVelocityBoundaryConditions(bbb_grid)
        T_bcs = TracerBoundaryConditions(bbb_grid)

        @test u_bcs isa FieldBoundaryConditions
        @test u_bcs.x.left  isa NFBC
        @test u_bcs.x.right isa NFBC
        @test u_bcs.y.left  isa ZFBC
        @test u_bcs.y.right isa ZFBC
        @test u_bcs.z.left  isa ZFBC
        @test u_bcs.z.right isa ZFBC

        @test v_bcs isa FieldBoundaryConditions
        @test v_bcs.x.left  isa ZFBC
        @test v_bcs.x.right isa ZFBC
        @test v_bcs.y.left  isa NFBC
        @test v_bcs.y.right isa NFBC
        @test v_bcs.z.left  isa ZFBC
        @test v_bcs.z.right isa ZFBC

        @test w_bcs isa FieldBoundaryConditions
        @test w_bcs.x.left  isa ZFBC
        @test w_bcs.x.right isa ZFBC
        @test w_bcs.y.left  isa ZFBC
        @test w_bcs.y.right isa ZFBC
        @test w_bcs.z.left  isa NFBC
        @test w_bcs.z.right isa NFBC

        @test T_bcs isa FieldBoundaryConditions
        @test T_bcs.x.left  isa ZFBC
        @test T_bcs.x.right isa ZFBC
        @test T_bcs.y.left  isa ZFBC
        @test T_bcs.y.right isa ZFBC
        @test T_bcs.z.left  isa ZFBC
        @test T_bcs.z.right isa ZFBC
    end
end


text:
"""
"""
module IO
end


text:
using Test, WAVI

@testset "Melting" begin
    @testset "test binary input file melt rate construction and update" begin

        @info "Testing binary input file construction and update"
        grid = Grid()
        bed_elevation = -900 .* ones(grid.nx, grid.ny) #low bed so we can make it float everywhere
        initial_conditions = InitialConditions(initial_thickness = 100.0*ones(grid.nx, grid.ny)) #initialize to 100m thickn

        #write a binary file
        m = ones(grid.nx, grid.ny)
        #m .= hton.(m)
        filename =  "melt_test_file.bin"
        mfileID =  open(filename,"w")
          write(mfileID, m[:,:])
        close(mfileID)

        #create a melt rate model from this input
        binfile_melt_rate = BinfileMeltRate(input_filename =filename)
        @test binfile_melt_rate isa BinfileMeltRate

        #create model
        model = Model(grid = grid,
                      bed_elevation = bed_elevation,
                      initial_conditions = initial_conditions,
                      melt_rate = binfile_melt_rate,
                      solver_params = SolverParams(maxiter_picard = 1))
        @test model.melt_rate isa BinfileMeltRate

        #check that melt rate read OK
        update_state!(model)
        @test all(model.fields.gh.basal_melt .== 1)

        #change the binary file and check we can update
        m = 2 .* ones(grid.nx, grid.ny)
        #m .= hton.(m)
        mfileID =  open(filename,"w")
          write(mfileID, m[:,:])
        close(mfileID)
        update_state!(model)
        @test all(model.fields.gh.basal_melt .== 2)

        #change the binary file to one of the wrong size, and check we return an error when updating
        m = 2 .* ones(grid.nx+1, grid.ny+1)
        m .= hton.(m)
        mfileID =  open(filename,"w")
          write(mfileID, m[:,:])
        close(mfileID)
        @test_throws DimensionMismatch update_state!(model)

        #remove the file we just made
        rm(filename)

    end

end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to the orthonormalization of DCMs.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/dcm.jl
# ==================

# Functions: orthonormalize
# -------------------------

@testset "DCM orthonormalization" begin
    for T in (Float32, Float64)
        α = 1 + rand(T)
        D = rand(DCM{T})

        v₁ = D[:, 1]
        v₂ = D[:, 2]
        v₃ = D[:, 3]

        D1 = DCM(hcat(v₁ * α, v₂,     v₃    ))
        D2 = DCM(hcat(v₁,     v₂ * α, v₃    ))
        D3 = DCM(hcat(v₁,     v₂,     v₃ * α))
        D4 = DCM(hcat(v₁ * α, v₂ * α, v₃    ))
        D5 = DCM(hcat(v₁,     v₂ * α, v₃ * α))
        D6 = DCM(hcat(v₁ * α, v₂,     v₃ * α))

        @test norm(D - orthonormalize(D1)) ≈ 0 atol = 10 * eps(T)
        @test norm(D - orthonormalize(D2)) ≈ 0 atol = 10 * eps(T)
        @test norm(D - orthonormalize(D3)) ≈ 0 atol = 10 * eps(T)
        @test norm(D - orthonormalize(D4)) ≈ 0 atol = 10 * eps(T)
        @test norm(D - orthonormalize(D5)) ≈ 0 atol = 10 * eps(T)
        @test norm(D - orthonormalize(D6)) ≈ 0 atol = 10 * eps(T)

        Do = orthonormalize(D1)
        @test eltype(Do) === T
    end
end


text:
using Shapefile
using Test
using RemoteFiles
import DBFTables
import Tables
import DataFrames

url_physical = "https://github.com/nvkelso/natural-earth-vector/raw/v4.1.0/110m_physical"
url_cultural = "https://github.com/nvkelso/natural-earth-vector/raw/v4.1.0/110m_cultural"
datadir = joinpath(@__DIR__, "data")

@RemoteFileSet natural_earth "Natural Earth 110m physical" begin
    # polygon
    ne_land_shp = @RemoteFile joinpath(url_physical, "ne_110m_land.shp") dir = datadir
    ne_land_shx = @RemoteFile joinpath(url_physical, "ne_110m_land.shx") dir = datadir
    ne_land_dbf = @RemoteFile joinpath(url_physical, "ne_110m_land.dbf") dir = datadir
    # linestring
    ne_coastline_shp = @RemoteFile joinpath(url_physical, "ne_110m_coastline.shp") dir = datadir
    ne_coastline_shx = @RemoteFile joinpath(url_physical, "ne_110m_coastline.shx") dir = datadir
    ne_coastline_dbf = @RemoteFile joinpath(url_physical, "ne_110m_coastline.dbf") dir = datadir
    # point
    ne_cities_shp = @RemoteFile joinpath(
        url_cultural,
        "ne_110m_populated_places_simple.shp",
    ) dir = datadir
    ne_cities_shx = @RemoteFile joinpath(
        url_cultural,
        "ne_110m_populated_places_simple.shx",
    ) dir = datadir
    ne_cities_dbf = @RemoteFile joinpath(
        url_cultural,
        "ne_110m_populated_places_simple.dbf",
    ) dir = datadir
end

download(natural_earth)

@testset "Tables interface" begin

@test isfile(natural_earth)

# create tables that we can use in the tests below
ne_land = Shapefile.Table(path(natural_earth, "ne_land_shp"))
ne_coastline = Shapefile.Table(path(natural_earth, "ne_coastline_shp"))
ne_cities = Shapefile.Table(path(natural_earth, "ne_cities_shp"))

@testset "Create from parts" begin
    ne_land_shp = open(path(natural_earth, "ne_land_shp")) do io
        read(io, Shapefile.Handle)
    end
    ne_land_dbf = open(path(natural_earth, "ne_land_dbf")) do io
        DBFTables.Table(io)
    end
    ne_land_parts = Shapefile.Table(ne_land_shp, ne_land_dbf)

    @test ne_land_shp isa Shapefile.Handle{Union{Shapefile.Polygon,Missing}}
    @test ne_land_dbf isa DBFTables.Table
    @test ne_land_parts isa Shapefile.Table
end

# without .shp extension it should also work
@test Shapefile.Table(splitext(path(natural_earth, "ne_land_shp"))[1]) isa Shapefile.Table

@testset "ne_land" begin
    @test ne_land isa Shapefile.Table{Union{Shapefile.Polygon,Missing}}
    @test length(ne_land) == 127
    @test count(true for r in ne_land) == length(ne_land)
    @test propertynames(ne_land) == [:featurecla, :scalerank, :min_zoom]
    @test propertynames(first(ne_land)) == [:featurecla, :scalerank, :min_zoom]
    @test ne_land.featurecla isa Vector{String}
    @test length(ne_land.scalerank) == length(ne_land)
    @test sum(ne_land.scalerank) == 58
    @test Shapefile.shapes(ne_land) isa Vector{Union{Shapefile.Polygon,Missing}}
    @test Tables.istable(ne_land)
    @test Tables.rows(ne_land) === ne_land
    @test Tables.columns(ne_land) === ne_land
    @test Tables.rowaccess(ne_land)
    @test Tables.columnaccess(ne_land)
    @test Tables.schema(ne_land) == Tables.Schema(
        (:featurecla, :scalerank, :min_zoom),
        (Union{String,Missing}, Union{Int,Missing}, Union{Float64,Missing}),
    )
    for r in ne_land
        @test Shapefile.shape(r) isa Shapefile.Polygon
        @test r.featurecla === "Land"
    end
    df_land = DataFrames.DataFrame(ne_land)
    @test size(df_land) == (127, 3)
    @test names(df_land) == [:featurecla, :scalerank, :min_zoom]
    df_land.featurecla isa Vector{String}
end

@testset "ne_coastline" begin
    @test ne_coastline isa Shapefile.Table{Union{Shapefile.Polyline,Missing}}
    @test length(ne_coastline) == 134
    @test count(true for r in ne_coastline) == length(ne_coastline)
    @test propertynames(ne_coastline) == [:scalerank, :featurecla, :min_zoom]
    @test propertynames(first(ne_coastline)) == [:scalerank, :featurecla, :min_zoom]
    @test ne_coastline.featurecla isa Vector{String}
    @test length(ne_coastline.scalerank) == length(ne_coastline)
    @test sum(ne_coastline.scalerank) == 59
    @test Shapefile.shapes(ne_coastline) isa Vector{Union{Shapefile.Polyline,Missing}}
    @test Tables.istable(ne_coastline)
    @test Tables.rows(ne_coastline) === ne_coastline
    @test Tables.columns(ne_coastline) === ne_coastline
    @test Tables.rowaccess(ne_coastline)
    @test Tables.columnaccess(ne_coastline)
    @test Tables.schema(ne_coastline) == Tables.Schema(
        (:scalerank, :featurecla, :min_zoom),
        (Union{Int,Missing}, Union{String,Missing}, Union{Float64,Missing}),
    )
    for r in ne_coastline
        @test Shapefile.shape(r) isa Shapefile.Polyline
        @test r.featurecla in ("Coastline", "Country")
    end
    df_coastline = DataFrames.DataFrame(ne_coastline)
    @test size(df_coastline) == (134, 3)
    @test names(df_coastline) == [:scalerank, :featurecla, :min_zoom]
    df_coastline.featurecla isa Vector{String}
end

@testset "ne_cities" begin
    @test ne_cities isa Shapefile.Table{Union{Shapefile.Point,Missing}}
    @test length(ne_cities) == 243
    @test count(true for r in ne_cities) == length(ne_cities)
    colnames = [
        :scalerank, :natscale, :labelrank, :featurecla, :name, :namepar, :namealt, :diffascii,
        :nameascii, :adm0cap, :capalt, :capin, :worldcity, :megacity, :sov0name, :sov_a3,
        :adm0name, :adm0_a3, :adm1name, :iso_a2, :note, :latitude, :longitude, :changed,
        :namediff, :diffnote, :pop_max, :pop_min, :pop_other, :rank_max, :rank_min, :geonameid,
        :meganame, :ls_name, :ls_match, :checkme, :min_zoom, :ne_id,
    ]
    @test propertynames(ne_cities) == colnames
    @test propertynames(first(ne_cities)) == colnames
    @test ne_cities.featurecla isa Vector{String}
    @test length(ne_cities.scalerank) == length(ne_cities)
    @test sum(ne_cities.scalerank) == 612
    @test Shapefile.shapes(ne_cities) isa Vector{Union{Shapefile.Point,Missing}}
    @test Tables.istable(ne_cities)
    @test Tables.rows(ne_cities) === ne_cities
    @test Tables.columns(ne_cities) === ne_cities
    @test Tables.rowaccess(ne_cities)
    @test Tables.columnaccess(ne_cities)
    classes = ["Admin-0 capital", "Admin-0 capital alt", "Populated place",
        "Admin-1 capital", "Admin-1 region capital", "Admin-0 region capital"]
    @test unique(ne_cities.featurecla) == classes
    for r in ne_cities
        @test Shapefile.shape(r) isa Shapefile.Point
        @test r.featurecla in classes
    end
    show_result = "Shapefile.Table{Union{Missing, Shapefile.Point}} with 243 rows and the following 38 columns:\n\t\nscalerank, natscale, labelrank, featurecla, name, namepar, namealt, diffascii, nameascii, adm0cap, capalt, capin, worldcity, megacity, sov0name, sov_a3, adm0name, adm0_a3, adm1name, iso_a2, note, latitude, longitude, changed, namediff, diffnote, pop_max, pop_min, pop_other, rank_max, rank_min, geonameid, meganame, ls_name, ls_match, checkme, min_zoom, ne_id\n"
    if VERSION < v"1.1"
        show_result = replace(show_result, "Shapefile.Point" => "Point")
    end
    @test sprint(
            show,
            ne_cities,
        ) === show_result
    df_cities = DataFrames.DataFrame(ne_cities)
    @test size(df_cities) == (243, 38)
    @test names(df_cities) == colnames
    df_cities.featurecla isa Vector{String}
end

# no need to use shx in Shapefile.Tables since we read the shapes into a Vector and can thus index them
# but since we have these files, we may as well try reading them
ne_shx = [path(rfs) for rfs in files(natural_earth) if endswith(path(rfs), ".shx")]
@test length(ne_shx) == 3
for shx_path in ne_shx
    open(path(natural_earth, "ne_land_shx")) do io
        read(io, Shapefile.IndexHandle)
    end
end

end  # testset "Tables interface"


text:
# Create BP Synthetic data w/ 20 Hz peak frequency and using
# the same geometry as the original data set, but w/o surface
# related multiples.
# Author: Philipp Witte, pwitte.slim@gmail.com
# Date: May 2018
#

# TO DO
# Set up path where data will be saved
data_path = "/path/to/data/"

using JUDI, SegyIO, JLD, PyPlot

# Load velocity
if !isfile("bp_synthetic_2004_true_velocity.jld")
    run(`wget ftp://slim.gatech.edu/data/SoftwareRelease/Imaging.jl/CompressiveLSRTM/bp_synthetic_2004_true_velocity.jld`)
end
vp = load(join([pwd(), "/bp_synthetic_2004_true_velocity.jld"]))["vp"] / 1f3

# Load density
if !isfile("bp_synthetic_2004_density.jld")
    run(`wget ftp://slim.gatech.edu/data/SoftwareRelease/Imaging.jl/CompressiveLSRTM/bp_synthetic_2004_density.jld`)
end
rho = load(join([pwd(), "/bp_synthetic_2004_density.jld"]))["rho"]

# Load geometry of original BP data set
if !isfile("bp_synthetic_2004_header_geometry.jld")
    run(`wget ftp://slim.gatech.edu/data/SoftwareRelease/Imaging.jl/CompressiveLSRTM/bp_synthetic_2004_header_geometry.jld`)
end
geometry = load(join([pwd(), "/bp_synthetic_2004_header_geometry.jld"]))

# Set up model structure
d = (6.25, 6.25)
o = (0., 0.)
m0 = (1f0 ./ vp).^2
n = size(m0)
model0 = Model(n, d, o, m0; rho=rho)

# Load header geometry from original data set
src_geometry = geometry["src_geometry"]
rec_geometry = geometry["rec_geometry"]
nsrc = geometry["nsrc"]

# Set up source
wavelet = ricker_wavelet(src_geometry.t[1], src_geometry.dt[1], 0.020)  # 27 Hz peak frequency
q = judiVector(src_geometry, wavelet)
###################################################################################################

# Save data to disk
opt = Options(limit_m = true,
              space_order = 16,
              buffer_size = 4000f0,
              save_data_to_disk = true,
              file_path = data_path,
              file_name = "bp_observed_data_")

# Setup operators
F = judiModeling(model0, q.geometry, rec_geometry; options=opt)

# Model data (write to disk)
F*q


text:
# # [A classification task when labels are known - determining the country of origin of cars given the cars characteristics](@id classification_tutorial)

# In this exercise we have some car technical characteristics (mpg, horsepower,weight, model year...) and the country of origin and we would like to create a model such that the country of origin can be accurately predicted given the technical characteristics.
# As the information to predict is a multi-class one, this is a _[classification]_(https://en.wikipedia.org/wiki/Statistical_classification) task.
# It is a challenging exercise due to the simultaneous presence of three factors: (1) presence of missing data; (2) unbalanced data - 254 out of 406 cars are US made; (3) small dataset.

#
# Data origin:
# - dataset description: [https://archive.ics.uci.edu/ml/datasets/auto+mpg](https://archive.ics.uci.edu/ml/datasets/auto+mpg)
#src Also useful: https://www.rpubs.com/dksmith01/cars
# - data source we use here: [https://archive.ics.uci.edu/ml/machine-learning-databases/auto-mpg/auto-mpg.data](https://archive.ics.uci.edu/ml/machine-learning-databases/auto-mpg/auto-mpg.data-original)

# Field description:

# 1. mpg:           _continuous_
# 2. cylinders:     _multi-valued discrete_
# 3. displacement:  _continuous_
# 4. horsepower:    _continuous_
# 5. weight:        _continuous_
# 6. acceleration:  _continuous_
# 7. model year:    _multi-valued discrete_
# 8. origin:        _multi-valued discrete_
# 9. car name:      _string (unique for each instance)_ - not used here

cd(@__DIR__)
using Pkg
Pkg.activate(".")
Pkg.add(["HTTP","CSV","DataFrames","Pipe","BenchmarkTools", "BetaML"])
Pkg.instantiate()

# ## Library and data loading

# We load a buch of packages that we'll use during this tutorial..
using Random, HTTP, CSV, DataFrames, BenchmarkTools,BetaML
#import DecisionTree, Flux
import Pipe: @pipe

# To load the data from the internet our workflow is
# (1) Retrieve the data --> (2) Clean it --> (3) Load it --> (4) Output it as a DataFrame.

# For step (1) we use `HTTP.get()`, for step (2) we use `replace!`, for steps (3) and (4) we uses the `CSV` package, and we use the "pip" `|>` operator to chain these operations:

urlDataOriginal = "https://archive.ics.uci.edu/ml/machine-learning-databases/auto-mpg/auto-mpg.data-original"
data = @pipe HTTP.get(urlDataOriginal).body                                                |>
             replace!(_, UInt8('\t') => UInt8(' '))                                        |>
             CSV.File(_, delim=' ', missingstring="NA", ignorerepeated=true, header=false) |>
             DataFrame;

# This results in a table where the rows are the observations (the various cars) and the column the fields. All BetaML models expect this layout.

# As the dataset is ordered, we randomly shuffle the data. Note that we pass to shuffle `copy(FIXEDRNG)` as the random nuber generator in order to obtain reproducible output ( [`FIXEDRNG`](@ref BetaML.Utils.FIXEDRNG) is nothing else than an istance of `StableRNG(123)` defined in the [`BetaML.Utils`](@ref utils_module) sub-module, but you can choose of course your own "fixed" RNG). See the [Dealing with stochasticity](@ref dealing_with_stochasticity) section in the [Getting started](@ref getting_started) tutorial for details.
data[shuffle(copy(FIXEDRNG),axes(data, 1)), :]
describe(data)

# Columns 1 to 7 contain  characteristics of the car, while column 8 encodes the country or origin ("1" -> US, "2" -> EU, "3" -> Japan). That's the variable we want to be able to predict.

# Columns 9 contains the car name, but we are not going to use this information in this tutorial.
# Note also that some fields have missing data.

# Our first step is hence to divide the dataset in features (the x) and the labels (the y) we want to predict. The `x` is then a Julia standard `Matrix` of 406 rows by 7 columns and the `y` is a vector of the 406 observations:
x     = Matrix{Union{Missing,Float64}}(data[:,1:7]);
y     = Vector{Int64}(data[:,8]);

# Some algorithms that we will use today don't work with missing data, so we need to _impute_ them. We use the [`predictMissing`](@ref) function provided by the [`BetaML.Clustering`](@ref clustering_module) sub-module. Internally the function uses a Gaussian Mixture Model to assign to the missing walue of a given record an average of the values of the non-missing records weighted for how close they are to our specific record.
# Note that the same function (`predictMissing`) can be used for Collaborative Filtering / recomendation systems. Using GMM has the advantage over traditional algorithms as k-nearest neighbors (KNN) that GMM can "detect" the hidden structure of the observed data, where some observation can be similar to a certain pool of other observvations for a certain characteristic, but similar to an other pool of observations for other characteristics.

xFull = predictMissing(x,rng=copy(FIXEDRNG)).X̂;

# Further, some models don't work with categorical data as such, so we need to represent our `y` as a matrix with a separate column for each possible categorical value (the so called "one-hot" representation).
# For example, within a three classes field, the individual value `2` (or `"Europe"` for what it matters) would be represented as the vector `[0 1 0]`, while `3` (or `"Japan"`) would become the vector `[0 0 1]`.
# To encode as one-hot we use the function [`oneHotEncoder`](@ref) in [`BetaML.Utils`](@ref utils_module)
y_oh  = oneHotEncoder(y);

# In supervised machine learning it is good practice to partition the available data in a _training_, _validation_, and _test_ subsets, where the first one is used to train the ML algorithm, the second one to train any eventual "hyper-parameters" of the algorithm and the _test_ subset is finally used to evaluate the quality of the algorithm.
# Here, for brevity, we use only the _train_ and the _test_ subsets, implicitly assuming we already know the best hyper-parameters. Please refer to the [regression tutorial](@ref regression_tutorial) for examples of how to use the validation subset to train the hyper-parameters, or even better the [clustering tutorial](@ref clustering_tutorial) for an example of using the [`crossValidation`](@ref) function.

# We use then the [`partition`](@ref) function in [BetaML.Utils](@ref utils_module), where we can specify the different data to partition (each matrix or vector to partition must have the same number of observations) and the shares of observation that we want in each subset. Here we keep 80% of observations for training (`xtrain`, `xTrainFull` and `ytrain`) and we use 20% of them for testing (`xtest`, `xTestFull` and `ytest`):
((xtrain,xtest),(xtrainFull,xtestFull),(ytrain,ytest),(ytrain_oh,ytest_oh)) = partition([x,xFull,y,y_oh],[0.8,1-0.8],rng=copy(FIXEDRNG));

# ## Random Forests

# We are now ready to use our first model, the Random Forests (in the [`BetaML.Trees`](@ref trees_module) sub-module). Random Forests build a "forest" of decision trees models and then average their predictions in order to make an overall prediction out of a feature vector.

# To "build" the forest model (i.e. to "train" it) we need to give the model the training feature matrix and the associated "true" training labels, and we need to specify the number of trees to employ (this is an example of hyper-parameters). Here we use 30 individual decision trees.

# As the labels are encoded using integers,  we need also to specify the parameter `forceClassification=true`, otherwise the model would undergo a _regression_ job instead.

myForest       = buildForest(xtrain,ytrain,30, rng=copy(FIXEDRNG),forceClassification=true);

# To obtain the predicted values, we can simply use the function [`BetaML.Trees.predict`](@ref)
#src [`predict`](@ref BetaML.Trees.predict)  [`predict`](@ref forest_prediction)
# with our `myForest` model and either the training or the testing data.
ŷtrain,ŷtest   = predict.(Ref(myForest), [xtrain,xtest],rng=copy(FIXEDRNG));
# Finally we can measure the _accuracy_ of our predictions with the [`accuracy`](@ref) function, with the sidenote that we need first to "parse" the ŷs as forcing the classification job transformed automatically them to strings (they originally were integers):
trainAccuracy,testAccuracy  = accuracy.([parse.(Int64,mode(ŷtrain,rng=copy(FIXEDRNG))),parse.(Int64,mode(ŷtest,rng=copy(FIXEDRNG)))],[ytrain,ytest])
#src (0.9969230769230769,0.8024691358024691)

@test testAccuracy > 0.8 #src

# The predictions are quite good, for the training set the algoritm predicted almost all cars' origins correctly, while for the testing set (i.e. those records that has **not** been used to train the algorithm), the correct prediction level is still quite high, at 80%

# While accuracy can sometimes suffice, we may often want to better understand which categories our model has trouble to predict correctly.
# We can investigate the output of a multi-class classifier more in-deep with a [`ConfusionMatrix`](@ref) where the true values (`y`) are given in rows and the predicted ones (`ŷ`) in columns, together to some per-class metrics like the _precision_ (true class _i_ over predicted in class _i_), the _recall_ (predicted class _i_ over the true class _i_) and others.

# We fist build the [`ConfusionMatrix`](@ref BetaML.Utils.ConfusionMatrix) object between `ŷ` and `y` and then we print it (we do it here for the test subset):

cm = ConfusionMatrix(parse.(Int64,mode(ŷtest,rng=copy(FIXEDRNG))),ytest,classes=[1,2,3],labels=["US","EU","Japan"])
print(cm;what="all")

# Because the printing of the confusion matrix employs `display`, the confusion matrix is printed on the script building this documentation pages rather than on the documentation pages themselves. Here it is:

#
# ```text
# 4×4 Matrix{Any}:
#  "Labels"    "US"    "EU"   "Japan"
#  "US"      43       3      3
#  "EU"       2      13      1
#  "Japan"    5       2      9
# 4×4 Matrix{Any}:
#  "Labels"   "US"      "EU"       "Japan"
#  "US"      0.877551  0.0612245  0.0612245
#  "EU"      0.125     0.8125     0.0625
#  "Japan"   0.3125    0.125      0.5625
# ```

# From the report we can see that Japanese cars have more trouble in being correctly classified, and in particular many Japanease cars are classified as US ones. This is likely a result of the class imbalance of the data set, and could be solved by balancing the dataset with various sampling tecniques before training the model.

# When we benchmark the resourse used (time and memory) we find that Random Forests remain pretty fast, expecially when we compare them with neural networks (see later)
@btime buildForest(xtrain,ytrain,30, rng=copy(FIXEDRNG),forceClassification=true);
#src   134.096 ms (781027 allocations: 196.30 MiB)

# ### Comparision with DecisionTree.jl

# DecisionTrees.jl random forests are similar in usage: we first "build" (train) the forest and we then make predictions out of the trained model.
# The main difference is that the model requires data with nonmissing values, so we are going to use the `xtrainFull` and `xtestFull` feature labels we created earlier:
## We train the model...
model = DecisionTree.build_forest(ytrain, xtrainFull,-1,30,rng=123)
## ..and we generate predictions and measure their error
(ŷtrain,ŷtest) = DecisionTree.apply_forest.([model],[xtrainFull,xtestFull]);
(trainAccuracy,testAccuracy) = accuracy.([ŷtrain,ŷtest],[ytrain,ytest])
#src (0.9969230769230769, 0.7530864197530864)

#src cm = ConfusionMatrix(ŷtest,ytest,classes=[1,2,3],labels=["US","EU","Japan"])
#src println(cm)
@test testAccuracy > 0.75 #src

# While the accuracy on the training set is exactly the same as for `BetaML` random forets, `DecisionTree.jl` random forests are slighly less accurate in the testing sample.
# Where however `DecisionTrees.jl` excell is in the efficiency: they are extremelly fast and memory thrifty, even if to this benchmark we should add the resources needed to impute the missing values.

# Also, one of the reasons DecisionTrees are such efficient is that internally they sort the data to avoid repeated comparision, but in this way they work only with features that are sortable, while BetaML random forests accept virtually any kind of input without the need of adapt it.
@btime  DecisionTree.build_forest(ytrain, xtrainFull,-1,30,rng=123);
#src 1.431 ms (10875 allocations: 1.52 MiB)

# ### Neural network

# Neural networks (NN) can be very powerfull, but have two "inconvenients" compared with random forests: first, are a bit "picky". We need to do a bit of work to provide data in specific format. Note that this is _not_ feature engineering. One of the advantages on neural network is that for the most this is not needed for neural networks. However we still need to "clean" the data. One issue is that NN don't like missing data. So we need to provide them with the feature matrix "clean" of missing data. Secondly, they work only with numerical data. So we need to use the one-hot encoding we saw earlier.
# Further, they work best if the features are scaled such that each feature has mean zero and standard deviation 1. We can achieve it with the function [`scale`](@ref) or, as in this case, [`getScaleFactors`](@ref).

xScaleFactors   = getScaleFactors(xtrainFull)
D               = size(xtrainFull,2)
classes         = unique(y)
nCl             = length(classes)

# The second "inconvenient" of NN is that, while not requiring feature engineering, they stil lneed a bit of practice on the way to build the network. It's not as simple as `train(model,x,y)`. We need here to specify how we want our layers, _chain_ the layers together and then decide a _loss_ overall function. Only when we done these steps, we have the model ready for training.
# Here we define 2 [`DenseLayer`](@ref) where, for each of them, we specify the number of neurons in input (the first layer being equal to the dimensions of the data), the output layer (for a classification task, the last layer output size beying equal to the number of classes) and an _activation function_ for each layer (default the `identity` function).
ls   = 50
l1   = DenseLayer(D,ls,f=relu,rng=copy(FIXEDRNG))
l2   = DenseLayer(ls,nCl,f=relu,rng=copy(FIXEDRNG))

# For a classification the last layer is a [`VectorFunctionLayer`](@ref) that has no learnable parameters but whose activation function is applied to the ensemble of the neurons, rather than individually on each neuron. In particular, for classification we pass the [`BetaML.Utils.softmax`](@ref) function whose output has the same size as the input (and the number of classes to predict), but we can use the `VectorFunctionLayer` with any function, including the [`pool1d`](@ref) function to create a "pooling" layer (using maximum, mean or whatever other subfunction we pass to `pool1d`)

l3   = VectorFunctionLayer(nCl,f=softmax) ## Add a (parameterless) layer whose activation function (softMax in this case) is defined to all its nodes at once

# Finally we _chain_ the layers and assign a loss function with [`buildNetwork`](@ref):
mynn = buildNetwork([l1,l2,l3],squaredCost,name="Multinomial logistic regression Model Cars") ## Build the NN and use the squared cost (aka MSE) as error function (crossEntropy could also be used)

# Now we can train our network using the function [`train!`](@ref). It has many options, have a look at the documentation for all the possible arguments.
# Note that we train the network based on the scaled feature matrix.
res  = train!(mynn,scale(xtrainFull,xScaleFactors),ytrain_oh,epochs=500,batchSize=8,optAlg=ADAM(),rng=copy(FIXEDRNG)) ## Use optAlg=SGD() to use Stochastic Gradient Descent instead

# Once trained, we can predict the label. As the trained was based on the scaled feature matrix, so must be for the predictions
(ŷtrain,ŷtest)  = predict.(Ref(mynn),[scale(xtrainFull,xScaleFactors),scale(xtestFull,xScaleFactors)])
trainAccuracy   = accuracy(ŷtrain,ytrain,rng=copy(FIXEDRNG))

#-

testAccuracy    = accuracy(ŷtest,ytest,rng=copy(FIXEDRNG))
#src (0.8953846153846153, 0.7654320987654321

@test testAccuracy > 0.76 #src

#src accuracy(mode(ŷtest,rng=copy(FIXEDRNG)),ytest)

cm = ConfusionMatrix(ŷtest,ytest,classes=[1,2,3],labels=["US","EU","Japan"],rng=copy(FIXEDRNG))
print(cm)

# 4×4 Matrix{Any}:
#  "Labels"    "US"    "EU"   "Japan"
#  "US"      44       0      5
#  "EU"       3      10      3
#  "Japan"    6       2      8
# 4×4 Matrix{Any}:
#  "Labels"   "US"      "EU"   "Japan"
#  "US"      0.897959  0.0    0.102041
#  "EU"      0.1875    0.625  0.1875
#  "Japan"   0.375     0.125  0.5

# We see a bit the limits of neural networks in this example. While NN can be extremelly performant in many domains, they also require lot of data and computational power, expecially considering the many possible hyper-parameters and hence its large space in the hyper-parameter tuning.
# In this example we arrive short to the performance of random forests, yet with a significant numberof neurons.

@btime train!(mynn,scale(xtrainFull),ytrain_oh,epochs=300,batchSize=8,rng=copy(FIXEDRNG),verbosity=NONE);
#src  11.841 s (62860672 allocations: 4.21 GiB)

# ### Comparisons with Flux

# In Flux the input must be in the form (fields, observations), so we transpose our original matrices
xtrainT, ytrain_ohT = transpose.([scale(xtrainFull,xScaleFactors), ytrain_oh])
xtestT, ytest_ohT   = transpose.([scale(xtestFull,xScaleFactors), ytest_oh])

# We define the Flux neural network model in a similar way than BetaML and load it with data, we train it, predict and measure the accuracies on the training and the test sets:

#src function poolForFlux(x,wsize=5)
#src     hcat([pool1d(x[:,i],wsize;f=maximum) for i in 1:size(x,2)]...)
#src end
Random.seed!(123)

l1         = Flux.Dense(D,ls,Flux.relu)
#l2         = Flux.Dense(ls,ls,Flux.relu)
l3         = Flux.Dense(ls,nCl,Flux.relu)
Flux_nn    = Flux.Chain(l1,l3)
loss(x, y) = Flux.logitcrossentropy(Flux_nn(x), y)
ps         = Flux.params(Flux_nn)
nndata     = Flux.Data.DataLoader((xtrainT, ytrain_ohT), batchsize=8,shuffle=true)
begin for i in 1:500  Flux.train!(loss, ps, nndata, Flux.ADAM()) end end
ŷtrain     = Flux.onecold(Flux_nn(xtrainT),1:3)
ŷtest      = Flux.onecold(Flux_nn(xtestT),1:3)
trainAccuracy =  accuracy(ŷtrain,ytrain)

#-

testAccuracy  = accuracy(ŷtest,ytest)
#src 0.9384615384615385, 0.7407407407407407
# While the train accuracy is little bit higher that BetaML, the test accuracy remains comparable

@test testAccuracy > 0.74 #src

# However the time is again lower than BetaML, even if here for "just" a factor 2
@btime begin for i in 1:500 Flux.train!(loss, ps, nndata, Flux.ADAM()) end end;
#src  5.665 s (8943640 allocations: 1.07 GiB)

# ## Summary

# This is the summary of the results we had trying to predict the country of origin of the cars, based on their technical characteristics:

# | Model                | Train acc     | Test Acc |  Training time (ms)* | Training mem (MB) *|
# |:-------------------- |:-------------:| --------:| ------------------- | ----------------- |
# | RF                   | 0.9969        | 0.8025   | 134                 | 196               |
# | RF (DecisionTree.jl) | 0.9969        | 0.7531   | 1.43                | 1.5               |
# | NN                   | 0.895         | 0.765    | 11841               | 4311              |
# | NN (Flux.jl)         | 0.938         | 0.741    | 5665                | 1096              |

# * on a Intel Core i5-8350U laptop

# We warn that this table just provides a rought idea of the various algorithms performances. Indeed there is a large amount of stochasticity both in the sampling of the data used for training/testing and in the initial settings of the parameters of the algorithm. For a statistically significant comparision we would have to repeat the analysis with multiple sampling (e.g. by cross-validation, see the [clustering tutorial](@ref clustering_tutorial) for an example) and initial random parameters.

# Neverthless the table above shows that, when we compare BetaML with the algorithm-specific leading packages, we found similar results in terms of accuracy, but often the leading packages are better optimised and run more efficiently (but sometimes at the cost of being less verstatile).
# Also, for this dataset, Random Forests seems to remain marginally more accurate than Neural Network, altought of course this depends on the hyper-parameters and, with a single run of the models, we don't know if this difference is significant.

max_stars_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_issues_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_forks_repo_name: sylvaticus/IntroSPMLJuliaCourse


text:
function compute_sdf(masks, out_shape)
    normalized_sdf = zeros(out_shape)
    for batch in 1:size(out_shape)[5]
        posmask = masks[:,:,:,:, batch] .!= 0
        if any(posmask)

        end
    end
end

max_stars_repo_name: skoc/julia-project

max_issues_repo_name: skoc/julia-project

max_forks_repo_name: skoc/julia-project


text:
module DexApiClients
using gRPCClient

include("api.jl")
using .api

import Base: show

# begin service: api.Dex

export DexBlockingClient, DexClient

struct DexBlockingClient
    controller::gRPCController
    channel::gRPCChannel
    stub::DexBlockingStub

    function DexBlockingClient(api_base_url::String; kwargs...)
        controller = gRPCController(; kwargs...)
        channel = gRPCChannel(api_base_url)
        stub = DexBlockingStub(channel)
        new(controller, channel, stub)
    end
end

struct DexClient
    controller::gRPCController
    channel::gRPCChannel
    stub::DexStub

    function DexClient(api_base_url::String; kwargs...)
        controller = gRPCController(; kwargs...)
        channel = gRPCChannel(api_base_url)
        stub = DexStub(channel)
        new(controller, channel, stub)
    end
end

show(io::IO, client::DexBlockingClient) = print(io, "DexBlockingClient(", client.channel.baseurl, ")")
show(io::IO, client::DexClient) = print(io, "DexClient(", client.channel.baseurl, ")")

import .api: CreateClient
"""
    CreateClient

- input: api.CreateClientReq
- output: api.CreateClientResp
"""
CreateClient(client::DexBlockingClient, inp::api.CreateClientReq) = CreateClient(client.stub, client.controller, inp)
CreateClient(client::DexClient, inp::api.CreateClientReq, done::Function) = CreateClient(client.stub, client.controller, inp, done)

import .api: UpdateClient
"""
    UpdateClient

- input: api.UpdateClientReq
- output: api.UpdateClientResp
"""
UpdateClient(client::DexBlockingClient, inp::api.UpdateClientReq) = UpdateClient(client.stub, client.controller, inp)
UpdateClient(client::DexClient, inp::api.UpdateClientReq, done::Function) = UpdateClient(client.stub, client.controller, inp, done)

import .api: DeleteClient
"""
    DeleteClient

- input: api.DeleteClientReq
- output: api.DeleteClientResp
"""
DeleteClient(client::DexBlockingClient, inp::api.DeleteClientReq) = DeleteClient(client.stub, client.controller, inp)
DeleteClient(client::DexClient, inp::api.DeleteClientReq, done::Function) = DeleteClient(client.stub, client.controller, inp, done)

import .api: CreatePassword
"""
    CreatePassword

- input: api.CreatePasswordReq
- output: api.CreatePasswordResp
"""
CreatePassword(client::DexBlockingClient, inp::api.CreatePasswordReq) = CreatePassword(client.stub, client.controller, inp)
CreatePassword(client::DexClient, inp::api.CreatePasswordReq, done::Function) = CreatePassword(client.stub, client.controller, inp, done)

import .api: UpdatePassword
"""
    UpdatePassword

- input: api.UpdatePasswordReq
- output: api.UpdatePasswordResp
"""
UpdatePassword(client::DexBlockingClient, inp::api.UpdatePasswordReq) = UpdatePassword(client.stub, client.controller, inp)
UpdatePassword(client::DexClient, inp::api.UpdatePasswordReq, done::Function) = UpdatePassword(client.stub, client.controller, inp, done)

import .api: DeletePassword
"""
    DeletePassword

- input: api.DeletePasswordReq
- output: api.DeletePasswordResp
"""
DeletePassword(client::DexBlockingClient, inp::api.DeletePasswordReq) = DeletePassword(client.stub, client.controller, inp)
DeletePassword(client::DexClient, inp::api.DeletePasswordReq, done::Function) = DeletePassword(client.stub, client.controller, inp, done)

import .api: ListPasswords
"""
    ListPasswords

- input: api.ListPasswordReq
- output: api.ListPasswordResp
"""
ListPasswords(client::DexBlockingClient, inp::api.ListPasswordReq) = ListPasswords(client.stub, client.controller, inp)
ListPasswords(client::DexClient, inp::api.ListPasswordReq, done::Function) = ListPasswords(client.stub, client.controller, inp, done)

import .api: GetVersion
"""
    GetVersion

- input: api.VersionReq
- output: api.VersionResp
"""
GetVersion(client::DexBlockingClient, inp::api.VersionReq) = GetVersion(client.stub, client.controller, inp)
GetVersion(client::DexClient, inp::api.VersionReq, done::Function) = GetVersion(client.stub, client.controller, inp, done)

import .api: ListRefresh
"""
    ListRefresh

- input: api.ListRefreshReq
- output: api.ListRefreshResp
"""
ListRefresh(client::DexBlockingClient, inp::api.ListRefreshReq) = ListRefresh(client.stub, client.controller, inp)
ListRefresh(client::DexClient, inp::api.ListRefreshReq, done::Function) = ListRefresh(client.stub, client.controller, inp, done)

import .api: RevokeRefresh
"""
    RevokeRefresh

- input: api.RevokeRefreshReq
- output: api.RevokeRefreshResp
"""
RevokeRefresh(client::DexBlockingClient, inp::api.RevokeRefreshReq) = RevokeRefresh(client.stub, client.controller, inp)
RevokeRefresh(client::DexClient, inp::api.RevokeRefreshReq, done::Function) = RevokeRefresh(client.stub, client.controller, inp, done)

import .api: VerifyPassword
"""
    VerifyPassword

- input: api.VerifyPasswordReq
- output: api.VerifyPasswordResp
"""
VerifyPassword(client::DexBlockingClient, inp::api.VerifyPasswordReq) = VerifyPassword(client.stub, client.controller, inp)
VerifyPassword(client::DexClient, inp::api.VerifyPasswordReq, done::Function) = VerifyPassword(client.stub, client.controller, inp, done)

# end service: api.Dex

end # module DexApiClients


text:
include("1.2.4.1 (a x^q+b x^n+c x^(2 n-q))^p.jl")
include("1.2.4.2 (d x)^m (a x^q+b x^n+c x^(2 n-q))^p.jl")
include("1.2.4.3 (d+e x^(n-q)) (a x^q+b x^n+c x^(2 n-q))^p.jl")
include("1.2.4.4 (f x)^m (d+e x^(n-q)) (a x^q+b x^n+c x^(2 n-q))^p.jl")

max_stars_repo_path: src/rules/1 Algebraic functions/1.2 Trinomial products/1.2.4 Improper/.jl


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Call Julia's builtin flisp-based parser. `offset` is 0-based offset into the
# byte buffer or string.
function fl_parse(text::Union{Core.SimpleVector,String},
                  filename::String, offset, options)
    if text isa Core.SimpleVector
        # Will be generated by C entry points jl_parse_string etc
        text, text_len = text
    else
        text_len = sizeof(text)
    end
    ccall(:jl_fl_parse, Any, (Ptr{UInt8}, Csize_t, Any, Csize_t, Any),
          text, text_len, filename, offset, options)
end

function fl_parse(text::AbstractString, filename::AbstractString, offset, options)
    fl_parse(String(text), String(filename), offset, options)
end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
using Test
using Distances
using Clustering
include("test_helpers.jl")

@testset "kmedoids() (k-medoids)" begin

@testset "Argument checks" begin
    Random.seed!(34568)
    @test_throws ArgumentError kmedoids(randn(2, 3), 1)
    @test_throws ArgumentError kmedoids(randn(2, 3), 4)
    dist = inv.(max.(pairwise(Euclidean(), randn(2, 3), dims=2), 0.1))
    @test kmedoids(dist, 2) isa KmedoidsResult
    @test_throws ArgumentError kmedoids(dist, 2, display=:mylog)
    for disp in keys(Clustering.DisplayLevels)
        @test kmedoids(dist, 2, display=disp) isa KmedoidsResult
    end
end

Random.seed!(34568)

d = 3
n = 200
k = 10

X = rand(d, n)
dist = pairwise(SqEuclidean(), X, dims=2)
@assert size(dist) == (n, n)

Random.seed!(34568)  # reset seed again to known state
R = kmedoids(dist, k)
@test isa(R, KmedoidsResult)
@test nclusters(R) == k
@test length(R.medoids) == length(unique(R.medoids))
@test all(a -> 1 <= a <= k, R.assignments)
@test R.assignments[R.medoids] == 1:k # Every medoid should belong to its own cluster
@test sum(counts(R)) == n
@test wcounts(R) == counts(R)
@test R.costs == dist[LinearIndices((n, n))[CartesianIndex.(R.medoids[R.assignments], 1:n)]]
@test isapprox(sum(R.costs), R.totalcost)
@test R.converged

@testset "Support for arrays other than Matrix{T}" begin
    @testset "$(typeof(M))" for M in equivalent_matrices(dist)
        Random.seed!(34568)  # restore seed as kmedoids is not determantistic
        R2 = kmedoids(M, k)
        @test R2.assignments == R.assignments
    end
end

# k=1 and k=n cases
x = pairwise(SqEuclidean(), [1 2 3; .1 .2 .3; 4 5.6 7], dims=2)
kmed1 = kmedoids(x, 1)
@test nclusters(kmed1) == 1
@test assignments(kmed1) == [1, 1, 1]
@test kmed1.medoids == [2]
kmed3 = kmedoids(x, 3)
@test nclusters(kmed3) == 3
@test sort(assignments(kmed3)) == [1, 2, 3]
@test sort(kmed3.medoids) == [1, 2, 3]

# this data set has three obvious groups:
# group 1: [1, 3, 4], values: [1, 2, 3]
# group 2: [2, 5, 7], values: [6, 7, 8]
# group 3: [6, 8, 9], values: [21, 20, 22]
#

X = reshape(map(Float64, [1, 6, 2, 3, 7, 21, 8, 20, 22]), 1, 9)
dist = pairwise(SqEuclidean(), X, dims=2)

R = kmedoids!(dist, [1, 2, 6])
@test isa(R, KmedoidsResult)
@test nclusters(R) == 3
@test R.medoids == [3, 5, 6]
@test R.assignments == [1, 2, 1, 1, 2, 3, 2, 3, 3]
@test counts(R) == [3, 3, 3]
@test wcounts(R) == counts(R)
@test R.costs ≈ [1, 1, 0, 1, 0, 0, 1, 1, 1]
@test R.totalcost ≈ 6.0
@test R.converged

end


text:
import MappedArrays
import LightGraphs

using MappedArrays: ReadonlyMappedArray

"""
    LightAutomaton{GT, ET} <: AbstractAutomaton

A hybrid automaton that uses the `LightGraphs` backend. See the constructor
[`LightAutomaton(::Int)`](@ref).

###  Fields

- `G` -- graph of type `GT` whose vertices determine the states
- `Σ` -- dictionary mapping the edges to their labels
"""
mutable struct LightAutomaton{GT, ET} <: AbstractAutomaton
    G::GT
    Σ::Dict{ET, Dict{Int, Int}}
    next_id::Int
    nt::Int
end

"""
    LightAutomaton(n::Int)

Creates a `LightAutomaton` with `n` states 1, 2, ..., `n`.
The automaton is initialized without any transitions,
use [`add_transition!`](@ref) to add transitions.

## Examples

To create an automaton with 2 nodes 1, 2, self-loops of labels 1, a transition
from 1 to 2 with label 2 and transition from 2 to 1 with label 3, do the
following:
```jldoctest
julia> a = LightAutomaton(2);

julia> add_transition!(a, 1, 1, 1) # Add a self-loop of label 1 for state 1
HybridSystems.LightTransition{LightGraphs.SimpleGraphs.SimpleEdge{Int64}}(Edge 1 => 1, 1)

julia> add_transition!(a, 2, 2, 1) # Add a self-loop of label 1 for state 2
HybridSystems.LightTransition{LightGraphs.SimpleGraphs.SimpleEdge{Int64}}(Edge 2 => 2, 2)

julia> add_transition!(a, 1, 2, 2) # Add a transition from state 1 to state 2 with label 2
HybridSystems.LightTransition{LightGraphs.SimpleGraphs.SimpleEdge{Int64}}(Edge 1 => 2, 3)

julia> add_transition!(a, 2, 1, 3) # Add a transition from state 2 to state 1 with label 3
HybridSystems.LightTransition{LightGraphs.SimpleGraphs.SimpleEdge{Int64}}(Edge 2 => 1, 4)
```
"""
function LightAutomaton(n::Int)
    G = LightGraphs.DiGraph(n)
    Σ = Dict{LightGraphs.edgetype(G), Dict{Int, Int}}()
    LightAutomaton(G, Σ, 0, 0)
end

struct LightTransition{ET} <: AbstractTransition
    edge::ET
    id::Int
end

"""
    struct LightTransitionIterator{GT, ET, VT}
        automaton::LightAutomaton{GT, ET}
        edge_iterator::VT
    end

Iterate over the transitions of `automaton` by iterating over the edges `edge`
of `edge_iterator` and the ids `id` of `automaton.Σ[edge]` for each one. Its
elements are `LightTransition(edge, id)`.
"""
struct LightTransitionIterator{GT, ET, VT}
    automaton::LightAutomaton{GT, ET}
    edge_iterator::VT
end
Base.eltype(::LightTransitionIterator{GT, ET}) where {GT, ET} = LightTransition{ET}
function Base.length(tit::LightTransitionIterator)
    eit = tit.edge_iterator
    # empty iterator must be handled separately (see #29)
    return isempty(eit) ? 0 : sum(edge -> length(tit.automaton.Σ[edge]), eit)
end
function new_id_iterate(tit::LightTransitionIterator, edge, edge_state, ids, ::Nothing)
    return new_edge_iterate(tit, iterate(tit.edge_iterator, edge_state))
end
function new_id_iterate(tit::LightTransitionIterator, edge, edge_state, ids, id_item_state)
    idσ, id_state = id_item_state
    return LightTransition(edge, idσ[1]), (edge, edge_state, ids, id_state)
end
new_edge_iterate(::LightTransitionIterator, ::Nothing) = nothing
function new_edge_iterate(tit::LightTransitionIterator, edge_item_state)
    edge, edge_state = edge_item_state
    ids = tit.automaton.Σ[edge]
    return new_id_iterate(tit, edge, edge_state, ids, iterate(ids))
end
function Base.iterate(tit::LightTransitionIterator)
    return new_edge_iterate(tit, iterate(tit.edge_iterator))
end
function Base.iterate(tit::LightTransitionIterator, edge_id_state)
    edge, edge_state, ids, id_state = edge_id_state
    return new_id_iterate(tit, edge, edge_state, ids, iterate(ids, id_state))
end

states(A::LightAutomaton) = LightGraphs.vertices(A.G)
nstates(A::LightAutomaton) = LightGraphs.nv(A.G)

transitiontype(A::LightAutomaton) = LightTransition{LightGraphs.edgetype(A.G)}
function transitions(A::LightAutomaton)
    return LightTransitionIterator(A, LightGraphs.edges(A.G))
end
ntransitions(A::LightAutomaton) = A.nt

function edge_object(A::LightAutomaton, q, r)
    @assert 1 <= q <= nstates(A)
    @assert 1 <= r <= nstates(A)
    return LightGraphs.Edge(q, r)
end

function add_transition!(A::LightAutomaton, q, r, σ)
    edge = edge_object(A, q, r)
    A.next_id += 1
    A.nt += 1
    id = A.next_id
    ids = get(A.Σ, edge, nothing)
    new_ids = ids === nothing
    if new_ids
        ids = Dict{Int, Int}()
    end
    ids[id] = σ
    if new_ids
        LightGraphs.add_edge!(A.G, edge)
        A.Σ[edge] = ids
    end
    return LightTransition(edge, id)
end
function has_transition(A::LightAutomaton, q, r)
    return LightGraphs.has_edge(A.G, edge_object(A, q, r))
end
function has_transition(A::LightAutomaton, t::LightTransition)
    return LightGraphs.has_edge(A.G, t.edge) && haskey(A.Σ[t.edge], t.id)
end
function rem_transition!(A::LightAutomaton, t::LightTransition)
    ids = A.Σ[t.edge]
    delete!(ids, t.id)
    A.nt -= 1
    if isempty(ids)
        LightGraphs.rem_edge!(A.G, t.edge)
        delete!(A.Σ, t.edge)
    end
end

function rem_state!(A::LightAutomaton, st)
    # We need to collect as we delete them in the loop
    for t in collect(in_transitions(A, st))
        rem_transition!(A, t)
    end
    for t in collect(out_transitions(A, st))
        rem_transition!(A, t)
    end
    return LightGraphs.rem_vertex!(A.G, st)
end

function add_state!(A::LightAutomaton)
    return LightGraphs.add_vertex!(A.G)
end

source(::LightAutomaton, t::LightTransition) = t.edge.src
event(A::LightAutomaton, t::LightTransition) = A.Σ[t.edge][t.id]
target(::LightAutomaton, t::LightTransition) = t.edge.dst

function in_transitions(A::LightAutomaton{GT, ET}, s) where {GT, ET}
    f(src) = edge_object(A, src, s)
    inneigh = LightGraphs.inneighbors(A.G, s)
    edges = ReadonlyMappedArray{ET, 1, typeof(inneigh), typeof(f)}(f, inneigh)
    LightTransitionIterator(A, edges)
end
function out_transitions(A::LightAutomaton{GT, ET}, s) where {GT, ET}
    f(dst) = edge_object(A, s, dst)
    outneigh = LightGraphs.outneighbors(A.G, s)
    edges = ReadonlyMappedArray{ET, 1, typeof(outneigh), typeof(f)}(f, outneigh)
    LightTransitionIterator(A, edges)
end

struct LightStateProperty{GT, ET, T} <: StateProperty{T}
    automaton::LightAutomaton{GT, ET}
    value::Vector{T}
end
function state_property_type(::Type{LightAutomaton{GT, ET}},
                             T::Type) where {GT, ET}
    return LightStateProperty{GT, ET, T}
end
function state_property(automaton::LightAutomaton, T::Type)
    return LightStateProperty(automaton, Vector{T}(undef, nstates(automaton)))
end
function Base.getindex(p::LightStateProperty,
                       s::Int)
    p.value[s]
end
function Base.setindex!(p::LightStateProperty, value,
                        s::Int)
    p.value[s] = value
end

struct LightTransitionProperty{GT, ET, T} <: TransitionProperty{T}
    automaton::LightAutomaton{GT, ET}
    value::Vector{T}
end
function transition_property_type(::Type{LightAutomaton{GT, ET}}, T::Type) where {GT, ET}
    return LightTransitionProperty{GT, ET, T}
end
function transition_property(automaton::LightAutomaton, T::Type)
    return LightTransitionProperty(automaton, Vector{T}(undef, automaton.next_id))
end
function Base.getindex(p::LightTransitionProperty,
                       t::LightTransition)
    p.value[t.id]
end
function Base.setindex!(p::LightTransitionProperty, value,
                        t::LightTransition)
    p.value[t.id] = value
end