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

using Markdown
using InteractiveUtils

# ╔═╡ 61f4e85b-3b4e-4906-85a6-fd4a08056532
using LinearAlgebra

# ╔═╡ 898608b9-09c8-4e6d-a0ee-29454d5ce0c9
using FFTW

# ╔═╡ 5ce3d306-5a65-4f07-86eb-c417d1959d4e
using BlockArrays

# ╔═╡ 7a6a3ee1-4fc8-4b1f-ba8a-5312fb2b4a1a
using Plots

# ╔═╡ 46dda2a7-9ca8-4c2e-b00e-e6aa95761f65
md"
# 3D problem

This debugging notebook compares Lorenzo's solver for 2D unit cells (in progress)
against Luke's solver for 1D unit cells.
"

# ╔═╡ aefce721-b5f9-4008-963a-898ac77aec66
md"
## Incident wave parameters
"

# ╔═╡ b75a1f3b-e5f6-4cb5-b2a9-550e07976f44
kᵢ = 10.0 # wavenumber k²=ω²μ₀ϵ₀ of incident wave

# ╔═╡ 0c5978ca-6b44-4bf0-b6e7-467bfa98ec60
θⁱ = π/2.0 # incidence angle (measured with respect to the x axis in xz plane)

# ╔═╡ a07150dd-1230-425e-b773-76cc0b03ba58
αᵢ = kᵢ * cos(θⁱ) # component of kᵢ along x (sheet)

# ╔═╡ a25e6d0f-eedd-40e9-ba76-984b72b19563
βᵢ = kᵢ * sin(θⁱ) # component of kᵢ along y (transmiision axis)

# ╔═╡ d0138354-2ce6-4d8c-b222-992f5df51b57
md"
## Metasurface impedance parameters

Using the metasurface parameters described in
[Pérez-Arancibia et al., 2018, Sideways adiabaticity: beyond ray optics for slowly varying metasurfaces](https://doi.org/10.1364/OE.26.030202)
which result in an incident wave being transmitted at an angle of choice,
while the reflection is attenuated.
"

# ╔═╡ 517389c5-8ea6-4406-9726-06d6438d9ceb
k = kᵢ / 1.1 # wavenumber designed for metasurface
# do not choose k = kᵢ (a.k.a Wood anomaly)

# ╔═╡ 2f1ab012-2b32-4cff-a6e7-1154a07a890a
θᵗ = -π/3    # transmitted field angle

# ╔═╡ e166637c-5160-41bc-b4d0-924d3da59985
d  = cos(θᵗ)-cos(θⁱ) # metasurface design parameter

# ╔═╡ 863dca2b-d32a-4cb1-8ba2-bb051bcd0b6b
L = (2*π)/(k*abs(d))  # Unit cell width

# ╔═╡ c2ea899f-29cf-44d2-b5f0-108b1b9bbb4c
ρₑˣˣ(x⃗) = 0.5sin(abs(θⁱ)) * (1 - exp.(im*k*d*x⃗[1]))

# ╔═╡ b551098c-c1c3-45c0-9129-99f555acbf4f
ρₑˣʸ(x⃗) = zero(eltype(x⃗[1]))

# ╔═╡ 34c99ebe-8429-401a-a101-c18d70ac2543
ρₑʸˣ(x⃗) = zero(eltype(x⃗[1]))

# ╔═╡ 4dc98adb-729d-4458-a6cc-9ee2bc2b73e6
ρₑʸʸ(x⃗) = zero(eltype(x⃗[1]))

# ╔═╡ dea64f1e-e74f-452e-9b04-2e8de9e0199c
σₘˣˣ(x⃗) = zero(eltype(x⃗[1]))

# ╔═╡ 5d492cf3-a5e9-4f17-8dee-ed14393ef003
σₘˣʸ(x⃗) = zero(eltype(x⃗[1]))

# ╔═╡ d6a00833-a98e-4f58-bb72-d12f5c0bf7d3
σₘʸˣ(x⃗) = zero(eltype(x⃗[1]))

# ╔═╡ ead74c7d-b9cf-40fc-b2d8-9f6357599d63
σₘʸʸ(x⃗) = 2sin(abs(θⁱ)) * (1 + exp(im*k*d*x⃗[1]))

# ╔═╡ 5f608c76-a381-468a-824b-5be3f1a1c009
md"
## Discretization parameters
"

# ╔═╡ 513288b8-a444-44f1-9e64-1eeff0b0da6c
Nx = 25 # number of discretization points / modes along unit cell

# ╔═╡ dfdcfa71-2876-4e3f-be1a-27730fcfc47c
Ny = 1

# ╔═╡ 4216696f-8958-4691-9f77-c9b77998ea55
Nxy = Nx * Ny

# ╔═╡ b6ee3a1c-790a-46d6-877f-6fd17a52c0f9
Lx = L # Unit cell width (integer multiple of period of metasurface)

# ╔═╡ b0869bb2-25a7-437a-900a-17ccef7162ba
Ly = L

# ╔═╡ 952ad83c-39ae-4f2c-beed-6124b1b8f3ac
dx = Lx/Nx

# ╔═╡ 6328f2b3-71eb-4a14-8868-36dee395b562
dy = Ly/Ny

# ╔═╡ 7a0d472a-705c-4699-8097-86894476c147
x⃗ = range(0, step=dx, length=Nx)

# ╔═╡ dbf104bd-4baf-4319-8aab-d9c1cffc1369
y⃗ = range(0, step=dy, length=Ny)

# ╔═╡ 4025b58f-f538-4492-8ce7-0b1ce5a41095
r⃗ = Iterators.product(x⃗, y⃗)

# ╔═╡ 18b42581-4df0-48ef-b2b2-c74332421b05
z⃗ = range(-L, L, length=2Nx)

# ╔═╡ 50a8c6d5-e45e-4393-8cb7-88235603908a
z⃗₁ = z⃗[z⃗ .< 0]

# ╔═╡ 88a871b9-4533-49f0-ab10-2544463c11f8
z⃗₂ = z⃗[z⃗ .> 0]

# ╔═╡ b6ca58f3-a916-4c80-b188-58fd07db8461
k⃗ˣ = 2*pi*fftfreq(Nx, 1/dx)

# ╔═╡ 16f10984-4df6-4c92-925e-a220aae156a3
k⃗ʸ = 2*pi*fftfreq(Ny, 1/dy)

# ╔═╡ 39e29cf4-e8ba-4350-b382-d9ee32e3963e
k⃗ = Iterators.product(k⃗ˣ, k⃗ʸ)

# ╔═╡ 202161a1-8f80-4d74-aeda-8433fbfe7378
β⃗ˣ = @. sqrt(Complex(kᵢ^2 - k⃗ˣ^2))

# ╔═╡ 2f205787-8ea0-419b-8972-a6aed1557e63
β⃗ˣʸ = map(e⃗ -> sqrt(Complex(kᵢ^2 - mapreduce(abs2, +, e⃗, init=zero(kᵢ)))), k⃗)

# ╔═╡ 7dcb06d9-8e27-4182-97b2-77eecf851e6e
i = findfirst(k⃗ˣ .- αᵢ .< eps()) # mode index of incident wave

# ╔═╡ 56a1ed22-f069-4926-abd9-1ef4c7e4e0ae
@assert i == findfirst(real.(β⃗ˣ .- βᵢ) .< eps())

# ╔═╡ 545c10ff-fa96-4807-80a8-1b065fb5e244
md"
## Helper functions
"

# ╔═╡ 68f9cc5b-d55b-42e7-a31f-b85b8ebc3820
"""
    diagFT(A::AbstractArray{<:Number, N}) where N
Accepts a diagonal endomorphism in N-D position space and returns the Fourier
transformed rank-2N tensor reshaped to a matrix with flattened indices
"""
function diagFT(A::AbstractArray{<:Number, ndim}) where ndim
    n⃗ = size(A)
    n = prod(n⃗)
    Ã = Array(reshape(Diagonal(reshape(A, n)), (n⃗..., n⃗...)))
    Ã = fft(ifft(Ã, 1:ndim), (ndim+1):2ndim)
    reshape(Ã, (n, n))
end

# ╔═╡ 5020f47a-61f9-4622-82a6-80bbd5e1c833
md"
## Luke's solver
"

# ╔═╡ d1a60c40-d46f-4ac9-9293-cd09946e4209
begin
	β̂ˣ = Diagonal(β⃗ˣ)
    ρ̃ₑˣˣ = diagFT(ρₑˣˣ.(Iterators.product(x⃗)))
    σ̃ₘʸʸ = diagFT(σₘʸʸ.(Iterators.product(x⃗)))
    A = [
        0.5β̂ˣ+kᵢ*ρ̃ₑˣˣ    -0.5β̂ˣ-kᵢ*ρ̃ₑˣˣ;
        2β̂ˣ+kᵢ*σ̃ₘʸʸ     2β̂ˣ+kᵢ*σ̃ₘʸʸ
    ]
    B = [
        0.5β̂ˣ-kᵢ*ρ̃ₑˣˣ    -0.5β̂ˣ+kᵢ*ρ̃ₑˣˣ;
        2β̂ˣ-kᵢ*σ̃ₘʸʸ     2β̂ˣ-kᵢ*σ̃ₘʸʸ
    ]
	S = A\B
	if θⁱ < 0
		I₁ = zeros(Nx)
		I₂ = [n==i ? 1 : 0 for n in 1:Nx]
	else
		I₁ = [n==i ? 1 : 0 for n in 1:Nx]
		I₂ = zeros(Nx)
	end
	inc = [I₁; I₂]
	scat = S*inc
	O₁ = scat[1:Nx]
	O₂ = scat[Nx+1:2Nx]
end;

# ╔═╡ eb3dd186-3741-42ee-9515-ace5b7126a86
begin
	# quick plotting
	BO₁ = rotr90(bfft(exp.(-β⃗ˣ * transpose(im * z⃗₁)) .* O₁, 1))
	BI₂ = rotr90(bfft(exp.(-β⃗ˣ * transpose(im * z⃗₂)) .* I₂, 1))
	BI₁ = rotr90(bfft(exp.( β⃗ˣ * transpose(im * z⃗₁)) .* I₁, 1))
	BO₂ = rotr90(bfft(exp.( β⃗ˣ * transpose(im * z⃗₂)) .* O₂, 1))
	fields = cat(BI₁ + BO₁, BI₂ + BO₂; dims=1)
	heatmap(x⃗, z⃗, real.(fields), xguide="x", yguide="z", aspect_ratio=:equal,color=:RdBu,clim=(-1.0,1.0))
end

# ╔═╡ 47bde8a6-2106-4baa-b472-8636e035baef
md"
## Lorenzo's solver
"

# ╔═╡ 8006c599-9451-4d85-b18e-adf2cd4eca56
R = [
	0I(Nxy)	-I(Nxy);
	1I(Nxy)	0I(Nxy)
]

# ╔═╡ 36e2b2e0-c5f4-4d9c-92f9-779828c44e2f
K = [
	Diagonal(vec(map(e⃗ -> e⃗[1]^2, k⃗)))	Diagonal(vec(map(e⃗ -> e⃗[1]*e⃗[2], k⃗)));
	Diagonal(vec(map(e⃗ -> e⃗[1]*e⃗[2], k⃗)))	Diagonal(vec(map(e⃗ -> e⃗[2]^2, k⃗)))
]

# ╔═╡ edf27dee-85ea-42e6-a489-a8fea95d3075
ρ̃ₑ = [
	diagFT(ρₑˣˣ.(r⃗))	diagFT(ρₑˣʸ.(r⃗));
	diagFT(ρₑʸˣ.(r⃗))	diagFT(ρₑʸʸ.(r⃗))
]

# ╔═╡ 9159265a-1858-42c8-8ab0-e48c3cfa86eb
ρ̃′ₑ = R' * ρ̃ₑ * R

# ╔═╡ 78eb0b10-c309-4da4-8bfc-75d1642b90f7
σ̃ₘ = [
	diagFT(σₘˣˣ.(r⃗))	diagFT(σₘˣʸ.(r⃗));
	diagFT(σₘʸˣ.(r⃗))	diagFT(σₘʸʸ.(r⃗))
]

# ╔═╡ b34d564a-7c80-4fa9-aa6d-14d30887d894
β̂ˣʸ = [
	Diagonal(vec(β⃗ˣʸ))	0.0im*I(Nxy);
	0.0im*I(Nxy) Diagonal(vec(β⃗ˣʸ))
]

# ╔═╡ 97bc0901-86ba-4f1d-b217-f25069977395
k⃗² = kᵢ^2*I(2Nxy)

# ╔═╡ da289e3e-283b-41b3-82cc-8342b7417cf4
AA = [
	(-(K - k⃗²) * ρ̃′ₑ + 0.5kᵢ*β̂ˣʸ)	((K - k⃗²) * ρ̃′ₑ - 0.5kᵢ*β̂ˣʸ);
	(-(K - k⃗²) * σ̃ₘ + 2kᵢ*β̂ˣʸ)		(-(K - k⃗²) * σ̃ₘ + 2kᵢ*β̂ˣʸ)
	# (-(R * K * R' - k⃗²) + 0.5kᵢ*β̂ˣʸ * σ̃ₘ)		((R * K * R' - k⃗²) - 0.5kᵢ*β̂ˣʸ * σ̃ₘ)
]

# ╔═╡ 111d9973-91bb-481c-b33e-067cfe11d0e1
BB = [
	((K - k⃗²) * ρ̃′ₑ + 0.5kᵢ*β̂ˣʸ)	(-(K - k⃗²) * ρ̃′ₑ - 0.5kᵢ*β̂ˣʸ);
	((K - k⃗²) * σ̃ₘ + 2kᵢ*β̂ˣʸ)		((K - k⃗²) * σ̃ₘ + 2kᵢ*β̂ˣʸ)
	# ((R * K * R' - k⃗²) + 0.5kᵢ*β̂ˣʸ * σ̃ₘ)		(-(R * K * R' - k⃗²) - 0.5kᵢ*β̂ˣʸ * σ̃ₘ)
]

# ╔═╡ 8a686e69-6b69-4e68-8b79-23140b5ee4dd
SS = AA\BB

# ╔═╡ 1d3b118a-a8b0-4dad-aecb-a75e9ea3161b
begin
	### Define equivalent 2D inputs
	I₁ˣ = [n==1 ? 0 : 0 for n in 1:Nxy]
	I₂ˣ = [n==1 ? 0 : 0 for n in 1:Nxy]
	I₁ʸ = [n==1 ? 1 : 0 for n in 1:Nxy]
	I₂ʸ = [n==1 ? 0 : 0 for n in 1:Nxy]
end;

# ╔═╡ e2892425-ca03-4823-89be-8ebc4439206b
out = SS * vcat(I₁ˣ, I₁ʸ, I₂ˣ, I₂ʸ)

# ╔═╡ 877eca8b-9e16-4a7d-ba56-8d44569ffaa3
begin
	O₁ˣ = out[1:Nxy]
	O₁ʸ = out[Nxy+1:2Nxy]
	O₂ˣ = out[2Nxy+1:3Nxy]
	O₂ʸ = out[3Nxy+1:4Nxy]
end;

# ╔═╡ 209f25c0-aecd-4983-888e-ddaa4b513e0e
begin
	# quick plotting
	CO₁ = rotr90(bfft(exp.(-β⃗ˣ * transpose(im * z⃗₁)) .* O₁ʸ[1:Nx], 1))
	CI₂ = rotr90(bfft(exp.(-β⃗ˣ * transpose(im * z⃗₂)) .* I₂ʸ[1:Nx], 1))
	CI₁ = rotr90(bfft(exp.( β⃗ˣ * transpose(im * z⃗₁)) .* I₁ʸ[1:Nx], 1))
	CO₂ = rotr90(bfft(exp.( β⃗ˣ * transpose(im * z⃗₂)) .* O₂ʸ[1:Nx], 1))
	heatmap(x⃗, z⃗, real.(cat(CI₁ + CO₁, CI₂ + CO₂; dims=1)), xguide="x", yguide="z", aspect_ratio=:equal,color=:RdBu,clim=(-1.0,1.0))
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

[compat]
BlockArrays = "~0.16.7"
FFTW = "~1.4.5"
Plots = "~1.22.2"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[AbstractFFTs]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "485ee0867925449198280d4af84bdb46a2a404d0"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.0.1"

[[Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.1"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[ArrayLayouts]]
deps = ["FillArrays", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "623a32b87ef0b85d26320a8cc7e57ded707aef64"
uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
version = "0.7.5"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[BlockArrays]]
deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra"]
git-tree-sha1 = "57a7b84893e2e11739e345ecfef14e43d7ed1143"
uuid = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
version = "0.16.7"

[[Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.8+0"

[[Cairo_jll]]
deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "f2202b55d816427cd385a9a4f3ffb226bee80f99"
uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
version = "1.16.1+0"

[[ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "9995eb3977fbf67b86d0a0a0508e83017ded03f2"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.14.0"

[[ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "1a90210acd935f222ea19657f143004d2c2a1117"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.38.0"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[Contour]]
deps = ["StaticArrays"]
git-tree-sha1 = "9f02045d934dc030edad45944ea80dbd1f0ebea7"
uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
version = "0.5.7"

[[DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[EarCut_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "3f3a2501fa7236e9b911e0f7a588c657e822bb6d"
uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
version = "2.2.3+0"

[[Expat_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b3bfd02e98aedfa5cf885665493c5598c350cd2f"
uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
version = "2.2.10+0"

[[FFMPEG]]
deps = ["FFMPEG_jll"]
git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8"
uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
version = "0.4.1"

[[FFMPEG_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "Pkg", "Zlib_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
git-tree-sha1 = "d8a578692e3077ac998b50c0217dfd67f21d1e5f"
uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
version = "4.4.0+0"

[[FFTW]]
deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"]
git-tree-sha1 = "463cb335fa22c4ebacfd1faba5fde14edb80d96c"
uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
version = "1.4.5"

[[FFTW_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c6033cc3892d0ef5bb9cd29b7f2f0331ea5184ea"
uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
version = "3.3.10+0"

[[FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "7f6ad1a7f4621b4ab8e554133dade99ebc6e7221"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.12.5"

[[FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[Fontconfig_jll]]
deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "21efd19106a55620a188615da6d3d06cd7f6ee03"
uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
version = "2.13.93+0"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[FreeType2_jll]]
deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "87eb71354d8ec1a96d4a7636bd57a7347dde3ef9"
uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
version = "2.10.4+0"

[[FriBidi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91"
uuid = "559328eb-81f9-559d-9380-de523a88c83c"
version = "1.0.10+0"

[[GLFW_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libglvnd_jll", "Pkg", "Xorg_libXcursor_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll"]
git-tree-sha1 = "0c603255764a1fa0b61752d2bec14cfbd18f7fe8"
uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89"
version = "3.3.5+1"

[[GR]]
deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Printf", "Random", "Serialization", "Sockets", "Test", "UUIDs"]
git-tree-sha1 = "c2178cfbc0a5a552e16d097fae508f2024de61a3"
uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
version = "0.59.0"

[[GR_jll]]
deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt5Base_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "ef49a187604f865f4708c90e3f431890724e9012"
uuid = "d2c73de3-f751-5644-a686-071e5b155ba9"
version = "0.59.0+0"

[[GeometryBasics]]
deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
git-tree-sha1 = "58bcdf5ebc057b085e58d95c138725628dd7453c"
uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
version = "0.4.1"

[[Gettext_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046"
uuid = "78b55507-aeef-58d4-861c-77aaff3498b1"
version = "0.21.0+0"

[[Glib_jll]]
deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "74ef6288d071f58033d54fd6708d4bc23a8b8972"
uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
version = "2.68.3+1"

[[Graphite2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "344bf40dcab1073aca04aa0df4fb092f920e4011"
uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472"
version = "1.3.14+0"

[[Grisu]]
git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe"
version = "1.0.2"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "60ed5f1643927479f845b0135bb369b031b541fa"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.14"

[[HarfBuzz_jll]]
deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"]
git-tree-sha1 = "129acf094d168394e80ee1dc4bc06ec835e510a3"
uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566"
version = "2.8.1+1"

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[IntelOpenMP_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d979e54b71da82f3a65b62553da4fc3d18c9004c"
uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
version = "2018.0.3+2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[IterTools]]
git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.3.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[JpegTurbo_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d735490ac75c5cb9f1b00d8b5509c11984dc6943"
uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
version = "2.1.0+0"

[[LAME_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c"
uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d"
version = "3.100.1+0"

[[LZO_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6"
uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
version = "2.10.1+0"

[[LaTeXStrings]]
git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.2.1"

[[Latexify]]
deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
git-tree-sha1 = "a4b12a1bd2ebade87891ab7e36fdbce582301a92"
uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
version = "0.15.6"

[[LazyArtifacts]]
deps = ["Artifacts", "Pkg"]
uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[Libffi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "0b4a5d71f3e5200a7dff793393e09dfc2d874290"
uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
version = "3.2.2+1"

[[Libgcrypt_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"]
git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae"
uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
version = "1.8.7+0"

[[Libglvnd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll", "Xorg_libXext_jll"]
git-tree-sha1 = "7739f837d6447403596a75d19ed01fd08d6f56bf"
uuid = "7e76a0d4-f3c7-5321-8279-8d96eeed0f29"
version = "1.3.0+3"

[[Libgpg_error_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9"
uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
version = "1.42.0+0"

[[Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

[[Libmount_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73"
uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
version = "2.35.0+0"

[[Libtiff_jll]]
deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "340e257aada13f95f98ee352d316c3bed37c8ab9"
uuid = "89763e89-9b03-5906-acba-b20f662cd828"
version = "4.3.0+0"

[[Libuuid_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066"
uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
version = "2.36.0+0"

[[LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[MKL_jll]]
deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
git-tree-sha1 = "5455aef09b40e5020e1520f551fa3135040d4ed0"
uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
version = "2021.1.1+2"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "5a5bc6bf062f0f95e62d0fe0a2d99699fed82dd9"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.8"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Measures]]
git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f"
uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
version = "0.3.1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NaNMath]]
git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.5"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[Ogg_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7937eda4681660b4d6aeeecc2f7e1c81c8ee4e2f"
uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
version = "1.3.5+0"

[[OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[OpenSSL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "15003dcb7d8db3c6c857fda14891a539a8f2705a"
uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
version = "1.1.10+0"

[[Opus_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720"
uuid = "91d4177d-7536-5919-b921-800302f37372"
version = "1.3.2+0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[PCRE_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b2a7af664e098055a7529ad1a900ded962bca488"
uuid = "2f80f16e-611a-54ab-bc61-aa92de5b98fc"
version = "8.44.0+0"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "9d8c00ef7a8d110787ff6f170579846f776133a9"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.4"

[[Pixman_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b4f5d02549a10e20780a24fce72bea96b6329e29"
uuid = "30392449-352a-5448-841d-b1acce4e97dc"
version = "0.40.1+0"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlotThemes]]
deps = ["PlotUtils", "Requires", "Statistics"]
git-tree-sha1 = "a3a964ce9dc7898193536002a6dd892b1b5a6f1d"
uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
version = "2.0.1"

[[PlotUtils]]
deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "Statistics"]
git-tree-sha1 = "2537ed3c0ed5e03896927187f5f2ee6a4ab342db"
uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
version = "1.0.14"

[[Plots]]
deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "JSON", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs"]
git-tree-sha1 = "457b13497a3ea4deb33d273a6a5ea15c25c0ebd9"
uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
version = "1.22.2"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Qt5Base_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "xkbcommon_jll"]
git-tree-sha1 = "ad368663a5e20dbb8d6dc2fddeefe4dae0781ae8"
uuid = "ea2cea3b-5b76-57ae-a6ef-0a8af62496e1"
version = "5.15.3+0"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[RecipesBase]]
git-tree-sha1 = "44a75aa7a527910ee3d1751d1f0e4148698add9e"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.1.2"

[[RecipesPipeline]]
deps = ["Dates", "NaNMath", "PlotUtils", "RecipesBase"]
git-tree-sha1 = "7ad0dfa8d03b7bcf8c597f59f5292801730c55b8"
uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c"
version = "0.4.1"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Scratch]]
deps = ["Dates"]
git-tree-sha1 = "0b4b7f1393cff97c33891da2a0bf69c6ed241fda"
uuid = "6c6a2e73-6563-6170-7368-637461726353"
version = "1.1.0"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Showoff]]
deps = ["Dates", "Grisu"]
git-tree-sha1 = "91eddf657aca81df9ae6ceb20b959ae5653ad1de"
uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
version = "1.0.3"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3240808c6d463ac46f1c1cd7638375cd22abbccb"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.12"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StatsAPI]]
git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.0.0"

[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "8cbbc098554648c84f79a463c9ff0fd277144b6c"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.10"

[[StructArrays]]
deps = ["Adapt", "DataAPI", "StaticArrays", "Tables"]
git-tree-sha1 = "2ce41e0d042c60ecd131e9fb7154a3bfadbf50d3"
uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
version = "0.6.3"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "1162ce4a6c4b7e31e0e6b14486a6986951c73be9"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.5.2"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Wayland_jll]]
deps = ["Artifacts", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "3e61f0b86f90dacb0bc0e73a0c5a83f6a8636e23"
uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89"
version = "1.19.0+0"

[[Wayland_protocols_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll"]
git-tree-sha1 = "2839f1c1296940218e35df0bbb220f2a79686670"
uuid = "2381bf8a-dfd0-557d-9999-79630e7b1b91"
version = "1.18.0+4"

[[XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[XSLT_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"]
git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a"
uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
version = "1.1.34+0"

[[Xorg_libX11_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
git-tree-sha1 = "5be649d550f3f4b95308bf0183b82e2582876527"
uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
version = "1.6.9+4"

[[Xorg_libXau_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4e490d5c960c314f33885790ed410ff3a94ce67e"
uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
version = "1.0.9+4"

[[Xorg_libXcursor_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXfixes_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "12e0eb3bc634fa2080c1c37fccf56f7c22989afd"
uuid = "935fb764-8cf2-53bf-bb30-45bb1f8bf724"
version = "1.2.0+4"

[[Xorg_libXdmcp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fe47bd2247248125c428978740e18a681372dd4"
uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
version = "1.1.3+4"

[[Xorg_libXext_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3"
uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
version = "1.3.4+4"

[[Xorg_libXfixes_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "0e0dc7431e7a0587559f9294aeec269471c991a4"
uuid = "d091e8ba-531a-589c-9de9-94069b037ed8"
version = "5.0.3+4"

[[Xorg_libXi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXfixes_jll"]
git-tree-sha1 = "89b52bc2160aadc84d707093930ef0bffa641246"
uuid = "a51aa0fd-4e3c-5386-b890-e753decda492"
version = "1.7.10+4"

[[Xorg_libXinerama_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll"]
git-tree-sha1 = "26be8b1c342929259317d8b9f7b53bf2bb73b123"
uuid = "d1454406-59df-5ea1-beac-c340f2130bc3"
version = "1.1.4+4"

[[Xorg_libXrandr_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "34cea83cb726fb58f325887bf0612c6b3fb17631"
uuid = "ec84b674-ba8e-5d96-8ba1-2a689ba10484"
version = "1.5.2+4"

[[Xorg_libXrender_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96"
uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
version = "0.9.10+4"

[[Xorg_libpthread_stubs_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "6783737e45d3c59a4a4c4091f5f88cdcf0908cbb"
uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
version = "0.1.0+3"

[[Xorg_libxcb_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
git-tree-sha1 = "daf17f441228e7a3833846cd048892861cff16d6"
uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
version = "1.13.0+3"

[[Xorg_libxkbfile_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "926af861744212db0eb001d9e40b5d16292080b2"
uuid = "cc61e674-0454-545c-8b26-ed2c68acab7a"
version = "1.1.0+4"

[[Xorg_xcb_util_image_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "0fab0a40349ba1cba2c1da699243396ff8e94b97"
uuid = "12413925-8142-5f55-bb0e-6d7ca50bb09b"
version = "0.4.0+1"

[[Xorg_xcb_util_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll"]
git-tree-sha1 = "e7fd7b2881fa2eaa72717420894d3938177862d1"
uuid = "2def613f-5ad1-5310-b15b-b15d46f528f5"
version = "0.4.0+1"

[[Xorg_xcb_util_keysyms_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "d1151e2c45a544f32441a567d1690e701ec89b00"
uuid = "975044d2-76e6-5fbe-bf08-97ce7c6574c7"
version = "0.4.0+1"

[[Xorg_xcb_util_renderutil_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "dfd7a8f38d4613b6a575253b3174dd991ca6183e"
uuid = "0d47668e-0667-5a69-a72c-f761630bfb7e"
version = "0.3.9+1"

[[Xorg_xcb_util_wm_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "e78d10aab01a4a154142c5006ed44fd9e8e31b67"
uuid = "c22f9ab0-d5fe-5066-847c-f4bb1cd4e361"
version = "0.4.1+1"

[[Xorg_xkbcomp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxkbfile_jll"]
git-tree-sha1 = "4bcbf660f6c2e714f87e960a171b119d06ee163b"
uuid = "35661453-b289-5fab-8a00-3d9160c6a3a4"
version = "1.4.2+4"

[[Xorg_xkeyboard_config_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xkbcomp_jll"]
git-tree-sha1 = "5c8424f8a67c3f2209646d4425f3d415fee5931d"
uuid = "33bec58e-1273-512f-9401-5d533626f822"
version = "2.27.0+4"

[[Xorg_xtrans_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "79c31e7844f6ecf779705fbc12146eb190b7d845"
uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
version = "1.4.0+3"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "cc4bf3fdde8b7e3e9fa0351bdeedba1cf3b7f6e6"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.5.0+0"

[[libass_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "5982a94fcba20f02f42ace44b9894ee2b140fe47"
uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"
version = "0.15.1+0"

[[libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[libfdk_aac_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "daacc84a041563f965be61859a36e17c4e4fcd55"
uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280"
version = "2.0.2+0"

[[libpng_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
version = "1.6.38+0"

[[libvorbis_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"]
git-tree-sha1 = "c45f4e40e7aafe9d086379e5578947ec8b95a8fb"
uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a"
version = "1.3.7+0"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"

[[x264_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fea590b89e6ec504593146bf8b988b2c00922b2"
uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a"
version = "2021.5.5+0"

[[x265_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "ee567a171cce03570d77ad3a43e90218e38937a9"
uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76"
version = "3.5.0+0"

[[xkbcommon_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll", "Wayland_protocols_jll", "Xorg_libxcb_jll", "Xorg_xkeyboard_config_jll"]
git-tree-sha1 = "ece2350174195bb31de1a63bea3a41ae1aa593b6"
uuid = "d8fb68d0-12a3-5cfd-a85a-d49703b185fd"
version = "0.9.1+5"
"""

# ╔═╡ Cell order:
# ╟─46dda2a7-9ca8-4c2e-b00e-e6aa95761f65
# ╠═61f4e85b-3b4e-4906-85a6-fd4a08056532
# ╠═898608b9-09c8-4e6d-a0ee-29454d5ce0c9
# ╠═5ce3d306-5a65-4f07-86eb-c417d1959d4e
# ╠═7a6a3ee1-4fc8-4b1f-ba8a-5312fb2b4a1a
# ╟─aefce721-b5f9-4008-963a-898ac77aec66
# ╠═b75a1f3b-e5f6-4cb5-b2a9-550e07976f44
# ╠═0c5978ca-6b44-4bf0-b6e7-467bfa98ec60
# ╠═a07150dd-1230-425e-b773-76cc0b03ba58
# ╠═a25e6d0f-eedd-40e9-ba76-984b72b19563
# ╟─d0138354-2ce6-4d8c-b222-992f5df51b57
# ╠═517389c5-8ea6-4406-9726-06d6438d9ceb
# ╠═2f1ab012-2b32-4cff-a6e7-1154a07a890a
# ╠═e166637c-5160-41bc-b4d0-924d3da59985
# ╠═863dca2b-d32a-4cb1-8ba2-bb051bcd0b6b
# ╠═c2ea899f-29cf-44d2-b5f0-108b1b9bbb4c
# ╠═b551098c-c1c3-45c0-9129-99f555acbf4f
# ╠═34c99ebe-8429-401a-a101-c18d70ac2543
# ╠═4dc98adb-729d-4458-a6cc-9ee2bc2b73e6
# ╠═dea64f1e-e74f-452e-9b04-2e8de9e0199c
# ╠═5d492cf3-a5e9-4f17-8dee-ed14393ef003
# ╠═d6a00833-a98e-4f58-bb72-d12f5c0bf7d3
# ╠═ead74c7d-b9cf-40fc-b2d8-9f6357599d63
# ╟─5f608c76-a381-468a-824b-5be3f1a1c009
# ╠═513288b8-a444-44f1-9e64-1eeff0b0da6c
# ╠═dfdcfa71-2876-4e3f-be1a-27730fcfc47c
# ╠═4216696f-8958-4691-9f77-c9b77998ea55
# ╠═b6ee3a1c-790a-46d6-877f-6fd17a52c0f9
# ╠═b0869bb2-25a7-437a-900a-17ccef7162ba
# ╠═952ad83c-39ae-4f2c-beed-6124b1b8f3ac
# ╠═6328f2b3-71eb-4a14-8868-36dee395b562
# ╠═7a0d472a-705c-4699-8097-86894476c147
# ╠═dbf104bd-4baf-4319-8aab-d9c1cffc1369
# ╠═4025b58f-f538-4492-8ce7-0b1ce5a41095
# ╠═18b42581-4df0-48ef-b2b2-c74332421b05
# ╠═50a8c6d5-e45e-4393-8cb7-88235603908a
# ╠═88a871b9-4533-49f0-ab10-2544463c11f8
# ╠═b6ca58f3-a916-4c80-b188-58fd07db8461
# ╠═16f10984-4df6-4c92-925e-a220aae156a3
# ╠═39e29cf4-e8ba-4350-b382-d9ee32e3963e
# ╠═202161a1-8f80-4d74-aeda-8433fbfe7378
# ╠═2f205787-8ea0-419b-8972-a6aed1557e63
# ╠═7dcb06d9-8e27-4182-97b2-77eecf851e6e
# ╠═56a1ed22-f069-4926-abd9-1ef4c7e4e0ae
# ╟─545c10ff-fa96-4807-80a8-1b065fb5e244
# ╠═68f9cc5b-d55b-42e7-a31f-b85b8ebc3820
# ╟─5020f47a-61f9-4622-82a6-80bbd5e1c833
# ╠═d1a60c40-d46f-4ac9-9293-cd09946e4209
# ╠═eb3dd186-3741-42ee-9515-ace5b7126a86
# ╟─47bde8a6-2106-4baa-b472-8636e035baef
# ╠═8006c599-9451-4d85-b18e-adf2cd4eca56
# ╠═36e2b2e0-c5f4-4d9c-92f9-779828c44e2f
# ╠═edf27dee-85ea-42e6-a489-a8fea95d3075
# ╠═9159265a-1858-42c8-8ab0-e48c3cfa86eb
# ╠═78eb0b10-c309-4da4-8bfc-75d1642b90f7
# ╠═b34d564a-7c80-4fa9-aa6d-14d30887d894
# ╠═97bc0901-86ba-4f1d-b217-f25069977395
# ╠═da289e3e-283b-41b3-82cc-8342b7417cf4
# ╠═111d9973-91bb-481c-b33e-067cfe11d0e1
# ╠═8a686e69-6b69-4e68-8b79-23140b5ee4dd
# ╠═1d3b118a-a8b0-4dad-aecb-a75e9ea3161b
# ╠═e2892425-ca03-4823-89be-8ebc4439206b
# ╠═877eca8b-9e16-4a7d-ba56-8d44569ffaa3
# ╠═209f25c0-aecd-4983-888e-ddaa4b513e0e
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using PyCall

try
    pyimport("wandb")
    @info "Using Pre-Installed Wandb Version"
catch e
    try
        run(`$(PyCall.pyprogramname) -m pip install wandb`)
    catch ee
        if !(typeof(ee) <: PyCall.PyError)
            rethrow(ee)
        end
        @warn "Python dependencies not installed.\n" *
              "Either\n" *
              "- Rebuild `PyCall` to use Conda by running the following in Julia REPL " *
              "- `ENV[PYTHON]=\"\"; using Pkg; Pkg.build(\"PyCall\"); Pkg.build(\"Wandb\")\n" *
              "- Or install the dependencies by running `pip` - `pip install wandb>=0.11`"
    end
end

ext: jl


text:
module FakeTerminals

import REPL

mutable struct FakeTerminal <: REPL.Terminals.UnixTerminal
    in_stream::Base.IO
    out_stream::Base.IO
    err_stream::Base.IO
    hascolor::Bool
    raw::Bool
    FakeTerminal(stdin,stdout,stderr,hascolor=true) =
        new(stdin,stdout,stderr,hascolor,false)
end

REPL.Terminals.hascolor(t::FakeTerminal) = t.hascolor
REPL.Terminals.raw!(t::FakeTerminal, raw::Bool) = t.raw = raw
REPL.Terminals.size(t::FakeTerminal) = (24, 80)

end


text:
using Compat.Test
using Compat.LinearAlgebra
using Compat.SparseArrays
using MathProgBase

function linprogtest(solver; objtol = 1e-7, primaltol = 1e-6)
    @testset "Testing linprog and subfunctions with $solver" begin
        # min -x
        # s.t. 2x + y <= 1.5
        # x,y >= 0
        # solution is (0.75,0) with objval -0.75

        @testset "LP1" begin
            @testset "with buildlp and solvelp" begin
                m = buildlp([-1,0],[2 1],'<',1.5,solver)
                sol = solvelp(m)
                @test sol.status == :Optimal
                @test isapprox(sol.objval, -0.75, atol=objtol)
                @test isapprox(norm(sol.sol - [0.75,0.0]), 0.0, atol=primaltol)
            end

            @testset "with linprog" begin
                sol = linprog([-1,0],[2 1],'<',1.5,solver)
                @test sol.status == :Optimal
                @test isapprox(sol.objval, -0.75, atol=objtol)
                @test isapprox(norm(sol.sol - [0.75,0.0]), 0.0, atol=primaltol)
            end
        end

        @testset "LP2" begin
            @testset "with buildlp and solvelp" begin
                m = buildlp([-1,0],[2 1],'<',1.5,solver)
                sol = solvelp(m)
                @test sol.status == :Optimal
                @test isapprox(sol.objval, -0.75, atol=objtol)
                @test isapprox(norm(sol.sol - [0.75,0.0]), 0.0, atol=primaltol)
            end

            @testset "with linprog" begin
                sol = linprog([-1,0],sparse([2 1]),'<',1.5,solver)
                @test sol.status == :Optimal
                @test isapprox(sol.objval, -0.75, atol=objtol)
                @test isapprox(norm(sol.sol - [0.75,0.0]), 0.0, atol=primaltol)
            end
        end

        # test infeasible problem:
        # min x
        # s.t. 2x+y <= -1
        # x,y >= 0

        @testset "LP3 infeasible" begin
            @testset "with buildlp and solvelp" begin
                m = buildlp([1,0],[2 1],'<',-1,solver)
                sol = solvelp(m)
                @test sol.status == :Infeasible

                r = sol.attrs[:infeasibilityray]
                @test length(r) == 1
                r[1] = sol.attrs[:infeasibilityray][1]
                @test isapprox(r[1]/abs(r[1]), -1.0)
            end

            @testset "with linprog" begin
                sol = linprog([1,0],[2 1],'<',-1,solver)
                @test sol.status == :Infeasible

                r = sol.attrs[:infeasibilityray]
                @test length(r) == 1
                r[1] = sol.attrs[:infeasibilityray][1]
                @test isapprox(r[1]/abs(r[1]), -1.0)
            end
        end

        # test unbounded problem:
        # min -x-y
        # s.t. -x+2y <= 0
        # x,y >= 0

        @testset "LP4 unbounded" begin
            @testset "with buildlp and solvelp" begin
                m = buildlp([-1,-1],[-1 2],'<',[0],solver)
                sol = solvelp(m)
                @test sol.status == :Unbounded
            end

            @testset "with linprog" begin
                sol = linprog([-1,-1],[-1 2],'<',[0],solver)
                @test sol.status == :Unbounded
            end
        end

        # unbounded problem with unique ray:
        # min -x-y
        # s.t. x-y == 0
        # x,y >= 0

        @testset "LP4 unbounded with unique ray" begin
            @testset "with buildlp and solvelp" begin
                m = buildlp([-1,-1],[1 -1],'=',0,solver)
                sol = solvelp(m)
                @test sol.status == :Unbounded
                @test sol.attrs[:unboundedray][1] > 1e-7
                @test isapprox(sol.attrs[:unboundedray][1], sol.attrs[:unboundedray][2])
            end

            @testset "with linprog" begin
                sol = linprog([-1,-1],[1 -1],'=',0,solver)
                @test sol.status == :Unbounded
                @test sol.attrs[:unboundedray][1] > 1e-7
                @test isapprox(sol.attrs[:unboundedray][1], sol.attrs[:unboundedray][2])
            end
        end
    end
end


text:
@testset "tape: record" begin

    tape = Tape()
    xv = rand(3,4)
    yv = rand(4,3)
    x = record!(tape, Input, xv)
    y = record!(tape, Input, yv)
    u = record!(tape, Call, *, (x, y))
    c = record!(tape, Constant, 2.0)
    w = record!(tape, Bcast, *, (u, c))
    z = record!(tape, Call, sum, (u,); kwargs=Dict{Symbol,Any}(:dims=>1))

    uv = xv * yv
    wv = uv .* 2
    zv = sum(uv, dims=1)
    @test getvalue(z) == zv

end

@testset "tape: ops" begin

    tape = Tape()
    xv = rand(3,4)
    yv = rand(4,3)
    x = record!(tape, Input, xv)
    y = record!(tape, Input, yv)
    z = sum(x * y .* 2; dims=1)

    zv = sum(xv * yv .* 2; dims=1)
    @test getvalue(z) == zv

end


text:
const BACKEND = CUDABackend()

mutable struct Decoder
    worddict::Dict
    posdict::Dict
    nn
end

struct Sample
    w::Var
    c::Var
    batchdims_w
    batchdims_c
    p::Var
    heads
end

function create_batch(samples::Vector{Sample}, batchsize::Int)
    batches = Sample[]
    for i = 1:batchsize:length(samples)
        range = i:min(i+batchsize-1,length(samples))
        s = samples[range]
        sort!(s, by = x -> x.batchdims_w[1], rev=true)
        w = cat(2, map(x -> x.w.data, s)...)
        w = Var(BACKEND(w))
        c = cat(2, map(x -> x.c.data, s)...)
        c = Var(BACKEND(c))
        batchdims_w = cat(1, map(x -> x.batchdims_w, s)...)
        batchdims_c = cat(1, map(x -> x.batchdims_c, s)...)
        p = cat(2, map(x -> x.p.data, s)...)
        p = Var(BACKEND(p))
        heads = map(x -> Var(BACKEND(x.heads.data)), s)
        push!(batches, Sample(w,c,batchdims_w,batchdims_c,p,heads))
    end
    batches
end

function Decoder(config::Dict)
    chardict = Dict("UNKNOWN" => 1)
    words = h5read(config["wordvec_file"], "words")
    wordembeds = h5read(config["wordvec_file"], "vectors")
    worddict = Dict(words[i] => i for i=1:length(words))
    posdict = Dict{String,Int}()
    traindata = readconll(config["train_file"], worddict, chardict, posdict)[1:10000]
    testdata = readconll(config["test_file"], worddict, chardict, posdict)
    posembeds = Uniform(-0.01,0.01)(Float32, 50, length(posdict))
    nn = NN(wordembeds, posembeds)

    info("#Training examples:\t$(length(traindata))")
    info("#Testing examples:\t$(length(testdata))")
    info("#Words:\t$(length(worddict))")
    info("#Chars:\t$(length(chardict))")
    info("#POS-tags:\t$(length(posdict))")
    testdata = create_batch(testdata, 100)

    opt = SGD()
    batchsize = config["batchsize"]
    for epoch = 1:config["nepochs"]
        println("Epoch:\t$epoch")
        #opt.rate = learnrate / batchsize
        opt.rate = config["learning_rate"] / batchsize * sqrt(batchsize) / (1 + 0.05*(epoch-1))
        println("Learning rate: $(opt.rate)")

        shuffle!(traindata)
        batches = create_batch(traindata, batchsize)
        prog = Progress(length(batches))
        loss = 0.0
        for i in 1:length(batches)
            s = batches[i]
            y = nn(s, true)
            loss += sum(Array(y.data))
            params = gradient!(y)
            foreach(opt, params)
            ProgressMeter.next!(prog)
        end
        loss /= length(batches)
        println("Loss:\t$loss")

        # test
        println("Testing...")
        preds = Int[]
        golds = Int[]
        for s in testdata
            y = nn(s, false)
            append!(preds, y)
            for h in s.heads
                append!(golds, Array(h.data))
            end
        end
        length(preds) == length(golds) || throw("Length mismatch: $(length(preds)), $(length(golds))")
        count = 0
        for k = 1:length(preds)
            preds[k] == golds[k] && (count += 1)
        end
        acc = round(count/length(preds), 5)
        println("Acc:\t$acc")
        println()
    end
    Decoder(worddict, posdict, nn)
end

function initvocab(path::String)
    chardict = Dict{String,Int}()
    #tagdict = Dict{String,Int}()
    lines = open(readlines, path)
    for line in lines
        isempty(line) && continue
        items = Vector{String}(split(line,"\t"))
        word = strip(items[2])
        chars = Vector{Char}(word)
        for c in chars
            c = string(c)
            if haskey(chardict, c)
                chardict[c] += 1
            else
                chardict[c] = 1
            end
        end
        #tag = strip(items[2])
        #haskey(tagdict,tag) || (tagdict[tag] = length(tagdict)+1)
    end

    chars = String[]
    for (k,v) in chardict
        v >= 3 && push!(chars,k)
    end
    chardict = Dict(chars[i] => i for i=1:length(chars))
    chardict["UNKNOWN"] = length(chardict) + 1
    chardict
end

function readconll(path::String, worddict::Dict, chardict::Dict, posdict::Dict)
    samples = Sample[]
    words = String[]
    postags = String[]
    heads = Int[]
    unkword = worddict["UNKNOWN"]
    unkchar = chardict["UNKNOWN"]

    lines = open(readlines, path)
    push!(lines, "")
    for i = 1:length(lines)
        line = lines[i]
        if isempty(line)
            isempty(words) && continue
            wordids = Int[]
            charids = Int[]
            batchdims_w = [length(words)]
            batchdims_c = Int[]
            for w in words
                w0 = replace(lowercase(w), r"[0-9]", '0')
                id = get(worddict, w0, unkword)
                push!(wordids, id)

                chars = Vector{Char}(w)
                ids = map(chars) do c
                    get(chardict, string(c), unkchar)
                end
                append!(charids, ids)
                push!(batchdims_c, length(ids))
            end
            posids = Int[]
            for p in postags
                id = get!(posdict, p, length(posdict)+1)
                push!(posids, id)
            end

            w = Var(reshape(wordids,1,length(wordids)))
            c = Var(reshape(charids,1,length(charids)))
            p = Var(reshape(posids,1,length(posids)))
            h = Var(heads)
            push!(samples, Sample(w,c,batchdims_w,batchdims_c,p,h))
            empty!(words)
            empty!(postags)
            heads = Int[]
        else
            items = Vector{String}(split(line,"\t"))
            word = strip(items[2])
            @assert !isempty(word)
            push!(words, word)
            push!(postags, strip(items[6]))
            push!(heads, parse(Int,items[9]))
        end
    end
    samples
end


text:
push!(LOAD_PATH, "../src/")

using Documenter, DecisionMakingProblems

makedocs(
    modules = [DecisionMakingProblems],
    format = Documenter.HTML(),
    sitename = "DecisionMakingProblems.jl",
    pages = [
        ##############################################
        ## MAKE SURE TO SYNC WITH docs/src/index.md ##
        ##############################################
        "Basics" => [
            "index.md",
            # "get_started.md",
            # "concepts.md"
           ],
        # "Usage" => [
        #     "mdp.md",
        #     # "pomdp.md",
        #     # "simplegame.md",
        #     # "mg.md",
        #     # "pomg.md",
        #     # "decpomdp.md"
        #    ],

        "MDP Models" => [
            "mdp.md",
            "hexworld.md",
            "2048.md",
            "cart_pole.md",
            "mountain_car.md",
            "simple_lqr.md",
            "collision_avoidance.md"
           ],

        "POMDP Models" => [
            "pomdp.md",
            "crying_baby.md",
            "machine_replacement.md",
            "catch.md"
           ],

        "Simple Games" => [
            "simplegame.md",
            "prisoners_dilemma.md",
            "rock_paper_scissors.md",
            "travelers.md"
           ],

        "POMG Models" => [
            "pomg.md",
            "multicaregiver.md",
           ],

        "Markov Game" => [
            "mg.md",
            "predator_prey.md",
           ],

        "Dec-POMDP" => [
            "decpomdp.md",
            "collab_predator_prey.md",
           ],
        # "Concepts" => [
        #     "concepts.md"
        #    ],
    ]
)

# deploydocs(
#     repo = "github.com/sisl/Alg4DMProblems.jl.git",
# )


text:
export CuDiag

const CuDiag{ElT,VecT} = Diag{ElT,VecT} where {VecT<:CuArray{ElT}}
const NonuniformCuDiagTensor{ElT,N,StoreT,IndsT} =
  Tensor{ElT,N,StoreT,IndsT} where {StoreT<:CuDiag}

CuArray(D::NonuniformCuDiagTensor) = CuArray(dense(D))

function NDTensors.dense(T::NonuniformCuDiagTensor{ElT}) where {ElT}
  R_data = CUDA.zeros(ElT, dim(inds(T)))
  diag_inds = diagind(reshape(R_data, dims(inds(T))), 0)
  R_data[diag_inds] .= data(store(T))
  return Tensor(Dense(R_data), inds(T))
end

function Base.promote_rule(
  ::Type{Diag{ElT2,VecT2}}, ::Type{CuDense}
) where {ElT2,VecT2<:Number}
  return promote_type(DenseT1, ElT2)
end

function Base.promote_rule(
  ::Type{<:Tensor{ElT1,N1,StoreT1}}, ::Type{<:Tensor{ElT2,N2,StoreT2}}
) where {ElT1,ElT2,N1,N2,StoreT1<:CuDense,StoreT2<:UniformDiag}
  ElT3 = promote_type(ElT1, ElT2)
  return Tensor{promote_type(ElT1, ElT2),N3,CuDense{ElT3,CuVector{ElT3}}} where {N3}
end
function Base.promote_rule(
  ::Type{<:Tensor{ElT1,N1,StoreT1}}, ::Type{<:Tensor{ElT2,N2,StoreT2}}
) where {ElT1,ElT2,N1,N2,StoreT1<:UniformDiag,StoreT2<:CuDense}
  ElT3 = promote_type(ElT1, ElT2)
  return Tensor{promote_type(ElT1, ElT2),N3,CuDense{ElT3,CuVector{ElT3}}} where {N3}
end

function Base.promote_rule(
  ::Type{<:UniformDiag{ElT1}}, ::Type{<:NonuniformDiag{ElT2,VecT2}}
) where {ElT1,ElT2,VecT2<:CuArray}
  ElT3 = promote_type(ElT1, ElT2)
  return NonuniformDiag{ElT3,CuVector{ElT3}}
end
function Base.promote_rule(
  ::Type{<:NonuniformDiag{ElT2,VecT2}}, ::Type{<:UniformDiag{ElT1}}
) where {ElT1,ElT2,VecT2<:CuArray}
  ElT3 = promote_type(ElT1, ElT2)
  return NonuniformDiag{ElT3,CuVector{ElT3}}
end
function Base.promote_rule(
  ::Type{DenseT1}, ::Type{Diag{ElT2,VecT2}}
) where {DenseT1<:CuDense,ElT2,VecT2<:Number}
  return promote_type(DenseT1, ElT2)
end

function contraction_output_type(
  TensorT1::Type{<:NonuniformCuDiagTensor}, TensorT2::Type{<:CuDenseTensor}, IndsR
)
  return similar_type(promote_type(TensorT1, TensorT2), IndsR)
end
function contraction_output_type(
  TensorT1::Type{<:CuDenseTensor}, TensorT2::Type{<:NonuniformCuDiagTensor}, IndsR
)
  return contraction_output_type(TensorT2, TensorT1, IndsR)
end

function contraction_output_type(
  TensorT1::Type{<:DiagTensor{<:Number,<:CuDiag}},
  TensorT2::Type{<:DiagTensor{<:Number,<:CuDiag}},
  IndsR::Type,
)
  return similar_type(promote_type(TensorT1, TensorT2), IndsR)
end

function contraction_output_type(
  TensorT1::Type{<:UniformDiagTensor},
  TensorT2::Type{<:DiagTensor{<:Number,<:CuDiag}},
  IndsR::Type,
)
  return similar_type(promote_type(TensorT1, TensorT2), IndsR)
end
function contraction_output_type(
  TensorT1::Type{<:DiagTensor{<:Number,<:CuDiag}},
  TensorT2::Type{<:UniformDiagTensor},
  IndsR::Type,
)
  return contraction_output_type(TensorT2, TensorT1, IndsR)
end

function zero_contraction_output(
  T1::UniformDiagTensor{ElT1,N1}, T2::CuDenseTensor{ElT2,N2}, indsR::IndsR
) where {ElT1,N1,ElT2,N2,IndsR}
  ElT3 = promote_type(ElT1, ElT2)
  dat = CUDA.zeros(ElT3, dim(indsR))
  return Tensor(Dense(dat), indsR)
end
function zero_contraction_output(
  T2::CuDenseTensor{ElT2,N2}, T1::UniformDiagTensor{ElT1,N1}, indsR::IndsR
) where {ElT1,N1,ElT2,N2,IndsR}
  return zero_contraction_output(T1, T2, indsR)
end

function zero_contraction_output(
  T1::TensorT1, T2::TensorT2, indsR
) where {TensorT1<:NonuniformDiagTensor,TensorT2<:NonuniformCuDiagTensor}
  ElT3 = promote_type(eltype(TensorT1), eltype(TensorT2))
  dat = CUDA.zeros(ElT3, length(data(store(T2))))
  return Tensor(Diag(dat), indsR)
end

function zero_contraction_output(
  T1::TensorT1, T2::TensorT2, indsR
) where {TensorT1<:UniformDiagTensor,TensorT2<:NonuniformCuDiagTensor}
  ElT3 = promote_type(eltype(TensorT1), eltype(TensorT2))
  dat = CUDA.zeros(ElT3, length(data(store(T2))))
  return Tensor(Diag(dat), indsR)
end
function zero_contraction_output(
  T1::TensorT1, T2::TensorT2, indsR
) where {TensorT2<:UniformDiagTensor,TensorT1<:NonuniformCuDiagTensor}
  return zero_contraction_output(T2, T1, indsR)
end

function zero_contraction_output(
  T1::TensorT1, T2::TensorT2, indsR
) where {TensorT1<:DiagTensor,TensorT2<:CuDenseTensor}
  ElT3 = promote_type(eltype(TensorT1), eltype(TensorT2))
  dat = CUDA.zeros(ElT3, length(data(store(T2))))
  return Tensor(Dense(dat), indsR)
end
function zero_contraction_output(
  T1::TensorT1, T2::TensorT2, indsR
) where {TensorT2<:DiagTensor,TensorT1<:CuDenseTensor}
  return zero_contraction_output(T2, T1, indsR)
end

function contraction_output(
  T1::UniformDiagTensor, T2::DiagTensor{Elt2,<:CuDiag}, indsR
) where {Elt2}
  return zero_contraction_output(T1, T2, indsR)
end
function contraction_output(
  T1::DiagTensor{Elt1,<:CuDiag}, T2::UniformDiagTensor, indsR
) where {Elt1}
  return contraction_output(T2, T1, indsR)
end

function contract!(
  C::CuDenseTensor{<:Number,NC},
  Clabels,
  A::UniformDiagTensor{<:Number,NA},
  Alabels,
  B::CuDenseTensor{<:Number,NB},
  Blabels,
) where {NC,NA,NB}
  Bstore = data(store(B))
  Astore = data(store(A))
  Cstore = data(store(C))
  return copyto!(Cstore, Astore .* Bstore)
end

function contract!(
  C::CuDenseTensor{<:Number,NC},
  Clabels,
  A::CuDenseTensor{<:Number,NA},
  Alabels,
  B::UniformDiagTensor{<:Number,NB},
  Blabels,
) where {NC,NA,NB}
  return contract!(C, Clabels, B, Blabels, A, Alabels)
end

function contract!(
  C::NonuniformCuDiagTensor{EltC,NC,IndsC},
  Clabels,
  A::UniformDiagTensor{EltA,NA,IndsA},
  Alabels,
  B::NonuniformCuDiagTensor{EltB,NB,IndsB},
  Blabels,
) where {EltC<:Number,EltB<:Number,EltA<:Number,NC,NB,NA,IndsA,IndsB,IndsC}
  Bstore = data(store(B))
  Astore = data(store(A))
  Cstore = data(store(C))
  return copyto!(Cstore, Astore .* Bstore)
end

function contract!(
  C::NonuniformCuDiagTensor{EltC,NC,IndsC},
  Clabels,
  B::NonuniformCuDiagTensor{EltB,NB,IndsB},
  Blabels,
  A::UniformDiagTensor{EltA,NA,IndsA},
  Alabels,
) where {EltC<:Number,EltB<:Number,EltA<:Number,NC,NB,NA,IndsA,IndsB,IndsC}
  Bstore = data(store(B))
  Astore = data(store(A))
  Cstore = data(store(C))
  return copyto!(Cstore, Astore .* Bstore)
end

function contract!(
  C::NonuniformCuDiagTensor{EltC,NC,IndsC},
  Clabels,
  B::NonuniformCuDiagTensor{EltB,NB,IndsB},
  Blabels,
  A::NonuniformCuDiagTensor{EltA,NA,IndsA},
  Alabels,
) where {EltC<:Number,EltB<:Number,EltA<:Number,NC,NB,NA,IndsA,IndsB,IndsC}
  Bstore = data(store(B))
  Astore = data(store(A))
  Cstore = data(store(C))
  return copyto!(Cstore, Astore .* Bstore)
end

# Dense * NonuniformCuDiag
function contract!(
  C::CuDenseTensor, Clabels, A::NonuniformCuDiagTensor, Alabels, B::CuDenseTensor, Blabels
)
  Astore = data(store(A))
  newAstore = CUDA.zeros(eltype(A), dims(inds(A))[1], dims(inds(A))[2])
  adi = diagind(newAstore, 0)
  newAstore[adi] .= Astore[:]
  newA = Tensor(Dense(vec(newAstore)), inds(A))
  return contract!(C, Clabels, newA, Alabels, B, Blabels)
end

function contract!(
  C::CuDenseTensor, Clabels, A::CuDenseTensor, Alabels, B::NonuniformCuDiagTensor, Blabels
)
  return contract!(C, Clabels, B, Blabels, A, Alabels)
end


text:
### Fetch required packages and reaction networks ###
using Catalyst, OrdinaryDiffEq, Random, Test
using ModelingToolkit: get_states, get_ps
include("test_networks.jl")

using StableRNGs
rng = StableRNG(12345)

### Checks the solutions of specific problems ###

# Exponential decay, should be identical to the (known) analytical solution.
exponential_decay = @reaction_network begin
    d, X → ∅
end d

for factor in [1e-2, 1e-1, 1e0, 1e1, 1e2]
    u0 = factor*rand(rng, length(get_states(exponential_decay)))
    p = factor*rand(rng, length(get_ps(exponential_decay)))
    prob = ODEProblem(exponential_decay,u0,(0.,100/factor),p)
    sol = solve(prob,Rosenbrock23(),saveat=range(0.,100/factor,length=101))
    analytic_sol = map(t -> u0[1]*exp(-p[1]*t),range(0.,100/factor,length=101))
    @test all(abs.(first.(sol.u) .- analytic_sol) .< 0.1)
end

# Networks with know equilibrium
known_equilibrium = @reaction_network begin
    (k1,k2), X1 ↔ X2
    (k3,k4), X3+ X4 ↔ X5
    (k5,k6), 2X6 ↔ 3X7
    (k7,k8), ∅ ↔ X8
end k1 k2 k3 k4 k5 k6 k7 k8

for factor in [1e-1, 1e0, 1e1, 1e2, 1e3]
    u0 = factor*rand(rng,length(get_states(known_equilibrium)))
    p = 0.01 .+ factor*rand(rng,length(get_ps(known_equilibrium)))
    prob = ODEProblem(known_equilibrium,u0,(0.,1000000.),p)
    sol = solve(prob,Rosenbrock23())
    @test abs.(sol.u[end][1]/sol.u[end][2] - p[2]/p[1]) < 10000*eps()
    @test abs.(sol.u[end][3]*sol.u[end][4]/sol.u[end][5] - p[4]/p[3]) < 10000*eps()
    @test abs.((sol.u[end][6]^2/factorial(2))/(sol.u[end][7]^3/factorial(3))- p[6]/p[5]) < 1e-8
    @test abs.(sol.u[end][8]  - p[7]/p[8]) < 10000*eps()
end

### Compares to the manually calcualted function ###
identical_networks_1 = Vector{Pair}()

function real_functions_1(du,u,p,t)
    X1,X2,X3 = u
    p1,p2,p3,k1,k2,k3,k4,d1,d2,d3 = p
    du[1] = p1 + k1*X2 - k2*X1*X3^2/factorial(2) - k3*X1 + k4*X3 - d1*X1
    du[2] = p2 - k1*X2 + k2*X1*X3^2/factorial(2) - d2*X2
    du[3] = p3 + 2*k1*X2 - 2*k2*X1*X3^2/factorial(2) + k3*X1 - k4*X3 - d3*X3
end
push!(identical_networks_1, reaction_networks_standard[1] => real_functions_1)

function real_functions_2(du,u,p,t)
    X1,X2 = u
    v1,K1,v2,K2,d = p
    du[1] = v1*K1/(K1+X2) - d*X1*X2
    du[2] = v2*X1/(K2+X1) - d*X1*X2
end
push!(identical_networks_1, reaction_networks_standard[2] => real_functions_2)

function real_functions_3(du,u,p,t)
    X1,X2,X3 = u
    v1,K1,n1,v2,K2,n2,v3,K3,n3,d1,d2,d3 = p
    du[1] = v1*K1^n1/(K1^n1+X3^n1) - d1*X1
    du[2] = v2*K2^n2/(K2^n2+X1^n2) - d2*X2
    du[3] = v3*K3^n3/(K3^n3+X2^n3) - d3*X3
end
push!(identical_networks_1, reaction_networks_hill[2] => real_functions_3)

function real_functions_4(du,u,p,t)
    X1,X2,X3 = u
    k1,k2,k3,k4,k5,k6 = p
    du[1] = -k1*X1 + k2*X2 + k5*X3 - k6*X1
    du[2] = -k3*X2 + k4*X3 + k1*X1 - k2*X2
    du[3] = -k5*X3 + k6*X1 + k3*X2 - k4*X3
end
push!(identical_networks_1, reaction_networks_constraint[1] => real_functions_4)

function real_functions_5(du,u,p,t)
    X,Y,Z = u
    k1,k2,k3,k4 = p
    du[1] = k1 - k2*log(12+X)*X
    du[2] = k2*log(12+X)*X - k3*log(3+Y)*Y
    du[3] = k3*log(3+Y)*Y - log(5,6+k4)*Z
end
push!(identical_networks_1, reaction_networks_weird[2] => real_functions_5)

for (i,networks) in enumerate(identical_networks_1)
    for factor in [1e-2, 1e-1, 1e0, 1e1]
        u0 = factor*rand(rng,length(get_states(networks[1])))
        p = factor*rand(rng,length(get_ps(networks[1])))
        (i==3) && (p = min.(round.(p).+1,10))                      #If parameter in exponent, want to avoid possibility of (-small u)^(decimal). Also avoid large exponents.
        prob1 = ODEProblem(networks[1],u0,(0.,10000.),p)
        sol1 = solve(prob1,Rosenbrock23(),saveat=1.)
        prob2 = ODEProblem(networks[2],u0,(0.,10000.),p)
        sol2 = solve(prob2,Rosenbrock23(),saveat=1.)
        @test all(abs.(hcat((sol1.u .- sol2.u)...)) .< 1e-7)
    end
end

### Tries solving a large number of problem, ensuring there are no errors. ###
for (i,network) in enumerate(reaction_networks_all)
    (i%5 == 0) && println("Iteration "*string(i)*" at line 104 in file solve_ODEs.jl")
    for factor in [1e-1, 1e0, 1e1]
        u0 = factor*rand(rng,length(get_states(network)))
        p = factor*rand(rng,length(get_ps(network)))
        in(i,[[11:20...]...,34,37,42]) && (p = min.(round.(p).+1,10))  #If parameter in exponent, want to avoid possibility of (-small u)^(decimal). Also avoid large exponents.
        prob = ODEProblem(network,u0,(0.,1.),p)
        @test solve(prob,Rosenbrock23()).retcode == :Success
    end
end

### No parameter test ###

no_param_network = @reaction_network begin
    (1.5,2), ∅ ↔ X
end
for factor in [1e0, 1e1, 1e2]
    u0 = factor*rand(rng,length(get_states(no_param_network)))
    prob = ODEProblem(no_param_network,u0,(0.,1000.))
    sol = solve(prob,Rosenbrock23())
    @test abs.(sol.u[end][1] - 1.5/2) < 1e-8
end


text:
using MambaModels
using Test

@test 1 == 1


text:
export DiGraph, plot
using LightGraphs
using TikzGraphs

import LightGraphs: DiGraph
function DiGraph(vtree::Vtree)
    N = num_nodes(vtree)
    g = DiGraph(N)
    nv = num_variables(vtree)
    id = nv
    dict = Dict{PlainVtreeInnerNode, Int}()
    if id < N
        dict[vtree] = (id += 1)
    end

    convert_rec(v::PlainVtreeLeafNode) = ()
    convert_rec(v::PlainVtreeInnerNode) = begin
        foreach(children(v)) do c
            c_id = isleaf(c) ? variable(c) : dict[c] = (id += 1)
            add_edge!(g, dict[v], c_id)
            convert_rec(c)
        end
    end

    convert_rec(vtree)
    @assert id == N
    label = [["$i" for i in 1 : nv];fill(".", nv - 1)]
    g, label
end

function DiGraph(lc::LogicCircuit; on_edge=noop, on_var=noop)
    nv = num_variables(lc)
    nn = num_nodes(lc)
    g = DiGraph(nn)
    id = 0
    dict = Dict{LogicCircuit, Int}()
    label = Vector{String}(undef, nn)
    dict[lc] = (id += 1)

    add_label!(g, dict, c::LogicCircuit) = begin
        label[dict[c]] =
        if isliteralgate(c) "$(literal(c))"
        elseif istrue(c)  "T"
        elseif isfalse(c)  "F"
        elseif is⋁gate(c) "⋁"
        else "⋀"
        end
    end

    if on_var == noop
        on_var = add_label!
    end

    foreach_down(lc) do n
        on_var(g, dict, n)
        if has_children(n)
            foreach(children(n)) do c
                c_id = haskey(dict, c) ? dict[c] : dict[c] = (id += 1)
                add_edge!(g, dict[n], c_id)
                on_edge(g, dict, n, c, dict[n], c_id)
            end
        end
    end
    g, label
end

import TikzGraphs: plot
plot(vtree::Vtree) = TikzGraphs.plot(DiGraph(vtree)...)
plot(lc::LogicCircuit) = TikzGraphs.plot(DiGraph(lc)...)


text:
# License for this file: MIT (expat)
# Copyright 2020-2021, DLR Institute of System Dynamics and Control
#
# This file is part of module ModiaResult (is included by the ModiaPlot_xxxMakie packages).

"""
    plot(result, names;
         heading = "", grid = true, xAxis = nothing,
         figure = 1, prefix = "", reuse = false, maxLegend = 10,
         minXaxisTickLabels = false,
         MonteCarloAsArea = true)

Generate **line plots** of the signals of the `result` data structure that are
identified with the `names` keys using the plot package defined with
`ModiaResult.@usePlotPackage(xxx)`. Possible values for `xxx`:
`"GLMakie", "WGLMakie", "CairoMakie", "PyPlot", "NoPlot", "SilentNoPlot"`).

`result` is any data structures that supports the abstract interface of `ModiaResult`.
Currently these are:

- `result::ModiaResult.ResultDict`
- `result::AbstractDict{String,T}`
- `result::DataFrames.DataFrame`
- `result` supports the `Tables.jl` interface with `Tables.columnaccess(result) = true`.

Argument `names` defines the diagrams to be drawn and the result data to be included in the respective diagram:

- If `names` is a **String**, generate one diagram with one time series of the variable with key `names`.

- If `names` is a **Tuple** of Strings, generate one diagram with the time series of the variables
  with the keys given in the tuple.

- If names is a **Vector** or a **Matrix** of **Strings** and/or **Tuples**,
  generate a vector or matrix of diagrams.

Note, the names (and their units, if available in the result) are automatically used as legends in the
respective diagram.

A signal variable identified by a `String` key can be a scalar of type `<:Number`
or an array of element type `<:Number`. A signal is defined by a vector of time values,
a corresponding vector of signal values, and the signal type (continuous or clocked).

Note, before passing data to the plot package,
it is converted to Float64. This allows to, for example, also plot rational numbers,
even if not supported by the plot package. `Measurements.Measurement{xxx}`
and `MonteCarloMeasurements` is specially handled.

# Optional Arguments
- `heading::AbstractString`: Optional heading above the diagram.

- `grid::Bool`: = true, to display a grid.

- `xAxis::Union{AbstractString,Nothing}`: Name of x-axis. If `xAxis=nothing`, the independent variable
  of the result (usually `"time"` is used as x-axis.

- `figure::Int`: Integer identifier of the window in which the diagrams shall be drawn.

- `prefix::AbstractString`: String that is appended in front of every legend label
  (useful especially if `reuse=true`).

- `reuse::Bool`: If figure already exists and reuse=false, clear the figure before adding the plot.
   Otherwise, include the plot in the existing figure without removing the curves present in the figure.
   `reuse = true` is ignored for `"WGLMakie"` (because not supported).

- `maxLegend::Int`: If the number of legend entries in one plot command `> maxLegend`,
  the legend is suppressed.
  All curves have still their names as labels. In PyPlot, the curves can be inspected by
  their names by clicking in the toolbar of the plot on button `Edit axis, curve ..`
  and then on `Curves`.

- `minXaxisTickLabels::Bool`: = true, if xaxis tick labels shall be
  removed in a vector or array of plots, if not the last row
  (useful when including plots in a document).
  = false, x axis tick labels are always shown (useful when interactively zooming into a plot).

- `MonteCarloAsArea::Bool`: = true, if MonteCarloMeasurements values are shown with the mean value
  and the area between the minimum and the maximum value of all particles.
  = false, if all particles of MonteCarloMeasurements values are shown (e.g. if a value has 2000 particles,
  then 2000 curves are shown in the diagram).

# Examples
```julia
import ModiaResult
using  Unitful

# Generate "using xxx" statement
# (where "xxx" is from a previous ModiaResult.usePlotPackage("xxx"))
ModiaResult.@usingModiaPlot

# Construct result data
t = range(0.0, stop=10.0, length=100);
result = Dict{String,Any}();
result["time"] = t*u"s";
result["phi"]  = sin.(t)*u"rad";
result["w"]    = cos.(t)*u"rad/s";
result["a"]    = 1.2*sin.(t)*u"rad/s^2";
result["r"]    = hcat(0.4 * cos.(t), 0.5 * sin.(t), 0.3*cos.(t))*u"m";

# 1 signal in one diagram (legend = "phi [rad]")
plot(result, "phi")

# 3 signals in one diagram
plot(result, ("phi", "w", "a"), figure=2)

# 3 diagrams in form of a vector (every diagram has one signal)
plot(result, ["phi", "w", "r"], figure=3)

# 4 diagrams in form of a matrix (every diagram has one signal)
plot(result, ["phi" "w";
              "a"   "r[2]" ], figure=4)

# 2 diagrams in form of a vector
plot(result, [ ("phi", "w"), ("a") ], figure=5)

# 4 diagrams in form of a matrix
plot(result, [ ("phi",)           ("phi", "w");
               ("phi", "w", "a")  ("r[2:3]",)     ],figure=6)

# Plot w=f(phi) in one diagram
plot(result, "w", xAxis="phi", figure=7)

# Append signal of the next simulation run to figure=1
# (legend = "Sim 2: phi [rad]")
result["phi"] = 0.5*result["phi"];
plot(result, "phi", prefix="Sim 2: ", reuse=true)
```

Example of a matrix of plots:

![Matrix of plots](../resources/images/matrix-of-plots.png)
"""
plot(result, names::AbstractString; kwargs...) = plot(result, [names]        ; kwargs...)
plot(result, names::Symbol        ; kwargs...) = plot(result, [string(names)]; kwargs...)
plot(result, names::Tuple         ; kwargs...) = plot(result, [names]        ; kwargs...)
plot(result, names::AbstractVector; kwargs...) = plot(result, reshape(names, length(names), 1); kwargs...)

"""
    showFigure(figure)

| Plot package |  Effect                             |
|:-------------|:------------------------------------|
| GLMakie      | Show `figure` in the single window. |
| WGLMakie     | Show `figure` in the single window. |
| CairoMakie   | Call is ignored                     |
| PyPlot       | Call is ignored                     |
| NoPlot       | Call is ignored                     |

# Example

```julia
import ModiaResult
ModiaResult.@usingModiaPlot
...
plot(..., figure=1)
plot(..., figure=2)
plot(..., figure=3)

showFigure(2)
showFigure(1)
```
"""
function showFigure end

"""
    saveFigure(figure, file; kwargs...)

Save figure on file. The file extension defines the image format
(for example `*.png`).

| Plot package |  Supported file extensions          |
|:-------------|:------------------------------------|
| GLMakie      | png, jpg, bmp                       |
| WGLMakie     | png                                 |
| CairoMakie   | png, pdf, svg, eps                  |
| PyPlot       | depends on [backend](https://matplotlib.org/stable/tutorials/introductory/usage.html) (png, pdf, jpg, tiff, svg, ps, eps) |
| NoPlot       | Call is ignored                     |

# Keyword arguments

- resolution: (width::Int, height::Int) of the scene in dimensionless
  units (equivalent to px for GLMakie and WGLMakie).

# Example

```julia
import ModiaResult
ModiaResult.@usingModiaPlot
...

plot(..., figure=1)
plot(..., figure=2)

saveFigure(1, "plot.png")   # save in png-format
saveFigure(2, "plot.svg")   # save in svg-format
```
"""
function saveFigure end

"""
    closeFigure(figure)

Close `figure`.
"""
function closeFigure end

max_stars_repo_name: ModiaSim/ModiaResult

max_issues_repo_name: ModiaSim/ModiaResult

max_forks_repo_name: ModiaSim/ModiaResult


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

# Compute X <- a X + b Y.
function matfun_axpby!(X,a,b,Y)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:m
        @inbounds for j=1:n
            if (b isa ValueOne)
                X[i,j]+=Y[i,j]
            else
                X[i,j]+=b*Y[i,j]
            end
        end
    end
end

@inline function exp_sid_m10(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_sid_m10!(A_copy)
end

@inline function exp_sid_m10!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=9
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    memslots9 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    value_one=ValueOne()
    # Computation order: Bb2 Ba2 B2 Ba3 B3 Ba4 B4 Ba5 B5 Bb6 B6 Ba7 Bb7 B7 Ba8 Bb8_7 B8 Ba9 Bb9 B9 B10 B11 T2k13
    # Computing Bb2 = x*I+x*A
    coeff1=0.0
    coeff2=0.125
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=0.0
    coeff2=0.125
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*A+x*B2
    coeff1=0.125
    coeff2=0.0
    memslots2 .= coeff1.*memslots1 .+ coeff2.*memslots4
    # Computing B3 with operation: mult
    mul!(memslots3,memslots2,memslots4)
    # Deallocating Ba3 in slot 2
    # Computing Ba4 = x*A+x*B3
    coeff1=0.125
    coeff2=0.0
    memslots2 .= coeff1.*memslots1 .+ coeff2.*memslots3
    # Computing B4 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Computing Ba5 = x*A+x*B4
    coeff1=0.125
    coeff2=0.0
    memslots2 .= coeff1.*memslots1 .+ coeff2.*memslots5
    # Computing B5 with operation: mult
    mul!(memslots6,memslots2,memslots5)
    # Deallocating Ba5 in slot 2
    # Computing Bb6 = x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=1.6421017771235788e-7
    coeff2=6.204734935438909e-8
    coeff3=2.957106114715868e-9
    coeff4=1.556371639324141e-10
    coeff5=1.556371639324141e-11
    memslots2 .= coeff1.*memslots1 .+ coeff2.*memslots4 .+ coeff3.*memslots3 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    # Computing B6 with operation: mult
    mul!(memslots7,memslots6,memslots2)
    # Deallocating Bb6 in slot 2
    # Computing Ba7 = x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=0.7439481579382581
    coeff2=0.4155284057336423
    coeff3=0.02479095151834799
    coeff4=0.001283057135586989
    coeff5=3.501669195497238e-5
    coeff6=1.0
    memslots2 .= coeff1.*memslots1 .+ coeff2.*memslots4 .+ coeff3.*memslots3 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    # Computing Bb7 = x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=0.03306559506631931
    coeff2=0.002630043177655382
    coeff3=0.0002100333647757715
    coeff4=3.7537107416419e-5
    coeff5=1.0
    memslots8 .= coeff1.*memslots4 .+ coeff2.*memslots3 .+ coeff3.*memslots5 .+ coeff4.*memslots6 .+ coeff5.*memslots7
    # Computing B7 with operation: mult
    mul!(memslots9,memslots2,memslots8)
    # Deallocating Ba7 in slot 2
    # Deallocating Bb7 in slot 8
    # Computing Ba8 = x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=1.3883617476073524
    coeff2=2.991654767354374
    coeff3=0.2857950268422422
    coeff4=0.03005135891320298
    coeff5=0.002742336655922557
    coeff6=61.75954247606858
    coeff7=1.0
    memslots2 .= coeff1.*memslots1 .+ coeff2.*memslots4 .+ coeff3.*memslots3 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots9
    # Computing Bb8_7 = x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=-3.2977952779222e-5
    coeff2=0.008139086096860678
    coeff3=0.001121744731945438
    coeff4=9.027588625491207e-5
    coeff5=8.572383602707347e-6
    coeff6=1.0
    # Smart lincomb recycle B6
    memslots7 .= coeff1.*memslots1 .+ coeff2.*memslots4 .+ coeff3.*memslots3 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    # Computing B8 with operation: mult
    mul!(memslots8,memslots2,memslots7)
    # Deallocating Ba8 in slot 2
    # Deallocating Bb8_7 in slot 7
    # Computing Ba9 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B8
    coeff1=1.0
    coeff2=0.125
    coeff3=0.5029302610017967
    coeff4=0.07705596948494946
    coeff5=0.004985549176118462
    coeff6=6.263526066651383e-5
    coeff7=1.0
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots3 .+ coeff5.*memslots5 .+ coeff6.*memslots6 .+ coeff7.*memslots8
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb9 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B8
    coeff1=1.0
    coeff2=0.125
    coeff3=0.5029302610017967
    coeff4=0.07705596948494946
    coeff5=0.004985549176118462
    coeff6=6.263526066651383e-5
    coeff7=1.0
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots3 .+ coeff5.*memslots5 .+ coeff6.*memslots6 .+ coeff7.*memslots8
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 3
    # Deallocating B4 in slot 5
    # Deallocating B5 in slot 6
    # Deallocating B8 in slot 8
    # Computing B9 with operation: mult
    mul!(memslots3,memslots2,memslots1)
    # Deallocating Ba9 in slot 2
    # Deallocating Bb9 in slot 1
    # Computing B10 with operation: mult
    mul!(memslots1,memslots3,memslots3)
    # Deallocating B9 in slot 3
    # Computing B11 with operation: mult
    mul!(memslots2,memslots1,memslots1)
    # Deallocating B10 in slot 1
    # Computing T2k13 = x*B7+x*B11
    coeff1=0.0
    coeff2=1.0
    # Smart lincomb recycle B7
    memslots9 .= coeff1.*memslots9 .+ coeff2.*memslots2
    # Deallocating B11 in slot 2
    return memslots9 # Returning T2k13
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
module SimulatedAnnealing

using Distributed
using Statistics

export SA, run, prun, geometric_schedule_fn, constant_Ebsf_done_fn, constant_E_done_fn,
    n_steps_done_fn, n_moves, n_accepted, n_max_moves!, n_max_moves, n_max_accepted,
    n_max_accepted!, set_current_solution!, current_solution, mean_E_total, mean_E_accepted,
    var_E_total, var_E_accepted, Ebsf, bsf, E, T

mutable struct SA{ET<:Real , ST<:Any}
    s::Vector{ST}     # current solution history
    s_best::Vector{ST} # best solution so far history
    T::Vector{AbstractFloat}      # Temperature history
    E::Vector{ET}      # Energy of current solution history
    Ebsf::Vector{ET} # Energy of best solution so far history
    n_max_moves::Int  # maximum number of moves before new T step
    n_max_accepted::Int # maximum number of accepted solutions before new T step
    step::Int         # current step number
    n_moves::Vector{Int}      # number of moves attempted per step
    n_accepted::Vector{Int}   # max. number of accepted states per step
    mean_E_total::Vector{Float64}  # average E of all tried neighbor candidates
    mean_E_accepted::Vector{Float64} # average E of all accepted neighbors
    var_E_total::Vector{Float64}  # average E of all tried neighbor candidates
    var_E_accepted::Vector{Float64} # average E of all accepted neighbors
    schedule::Function
    done::Function
    guess_T0::Bool
    p0::Float64
    fused_neighbor_energy::Bool
    energy_statistics::Bool
    solution_history::Bool
    report_every::Int
    current_solution::ST
    bsf::ST

    function SA{ET,ST}(;schedule, done, n_max_moves, n_max_accepted,
                      guess_T0 = true,
                      p0 = 0.5,
                      fused_neighbor_energy = false,
                      energy_statistics = true,
                      solution_history = false,
                      report_every = 0) where ET<:Real where ST

        new{ET,ST}(ST[], #s
                  ST[],   #s_best
                  Vector{Float64}(),    #T
                  Vector{ET}(),    #E
                  Vector{ET}(),    #Ebsf
                  n_max_moves,        #n_max_moves
                  n_max_accepted,     #n_max_accepted
                  0, #step
                  Int[], # n_moves
                  Int[],  #n_accepted
                  Vector{Float64}(), #mean_E_total
                  Vector{Float64}(), #mean_E_accepted
                  Vector{Float64}(), #var_E_total
                  Vector{Float64}(), #var_E_accepted
                  schedule,
                  done,
                  guess_T0,
                  p0,
                  fused_neighbor_energy,
                  energy_statistics,
                  solution_history,
                  report_every)
    end
end

""" returns total number of moves trialed in last annealing step"""
n_moves(sa::SA{ET,ST} where ET where ST) = sa.n_moves[end]

""" returns number of accepted moves in last annealing step"""
n_accepted(sa::SA{ET,ST} where ET where ST) = sa.n_accepted[end]

""" set max number of trialed moves in an annealing step"""
n_max_moves!(sa::SA{ET,ST} where ET where ST,n) = sa.n_max_moves = n

""" returns max number of trialed moves in an annealing step"""
n_max_moves(sa::SA{ET,ST} where ET where ST) = sa.n_max_moves

""" returns max number of accepted moves in an annealing step"""
n_max_accepted(sa::SA{ET,ST} where ET where ST) = sa.n_max_accepted

""" set max number of accepted moves in an annealing step"""
n_max_accepted!(sa::SA{ET,ST} where ET where ST,n) = sa.n_max_accepted = n

""" set the current solution"""
set_current_solution!(sa::SA{ET,ST} where ET where ST, solution) = sa.current_solution = solution

""" returns the current candidate solution"""
current_solution(sa::SA{ET,ST} where ET where ST) = sa.current_solution

""" returns the mean energy of all trialed moves in last annealing step"""
mean_E_total(sa::SA{ET,ST} where ET where ST) = sa.mean_E_total[end]

""" returns the mean energy of all accepted moves in last annealing step"""
mean_E_accepted(sa::SA{ET,ST} where ET where ST) = sa.mean_E_accepted[end]

""" returns the energy variance of all trialed moves in last annealing step"""
var_E_total(sa::SA{ET,ST} where ET where ST) = sa.var_E_total[end]

""" returns the energy variance of all accepted moves in last annealing step"""
var_E_accepted(sa::SA{ET,ST} where ET where ST) = sa.var_E_accepted[end]

""" returns the energy of the Best So Far solution over the whole annealing process"""
Ebsf(sa::SA{ET,ST} where ET where ST) = sa.Ebsf[end]

""" returns the Best So Far solution over the whole annealing process"""
bsf(sa::SA{ET,ST} where ET where ST) = sa.bsf

""" returns the final energy of the previous annealing step"""
E(sa::SA{ET,ST} where ET where ST) = sa.E[end]

""" returns the temperature of the previous annealing step"""
T(sa::SA{ET,ST} where ET where ST) = sa.T[end]

""" returns the history (vector of length = number of annealing steps so far) of
some quantity in annealing process sa. Symbol s determines the quantity:
:n_moves     - number of moves trialed per step
:n_accepted  - number of moves accepted per step
:mean_E_total - average energy of all moves trialed per step
:mean_E_accepted - average energy of all moves accepted per step
:var_E_total     - energy variance of all moves trialed per step
:var_E_accepted  - energy variance of all moves accepted per step
:Ebsf            - energy of Best solution So Far per step
:E               - final energy per annealing step
:T               - temperature per annealing step
"""
history(sa::SA{ET,ST} where ET where ST, s::Symbol) = getfield(sa, s)

function history(sa::SA{ET,ST} where ET where ST, s::Symbol, n::Int)
    v = getfield(sa,s)
    N = length(v)
    if n > 0
        @view v[n:N]
    else
        @view v[N + n : N]
    end
end

geometric_schedule_fn(alpha) = sa::SA -> sa.T[end] * alpha

function constant_E_done_fn(n_steps)
    function (sa::SA)
        if sa.step < n_steps
            false
        else
            N = sa.step
            E_now = sa.E[N]
            all(x->x ≈ E_now, sa.E[N-(n_steps-1):N-1])
        end
    end
end

function constant_Ebsf_done_fn(n_steps)
    function (sa::SA)
        if sa.step < n_steps
            false
        else
            N = sa.step
            E_now = sa.Ebsf[N]
            all(x->x ≈ E_now, sa.Ebsf[N-(n_steps-1):N-1])
        end
    end
end

function n_steps_done_fn(n_steps)
    ndone(sa::SA) = sa.step >= n_steps
end

function arithmetic_schedule_fn(a,b)
    function asched(sa::SA)
        a*sa.T[end] + b
    end
end

function initial_temperature(s0, energy::Function, neighbor::Function, p0)
    ΔE_avg = 0.0
    E = energy(s0)
    for i in 1:1000
        ΔE_avg += abs(E - (energy ∘ neighbor)(copy(s0)))
    end
    ΔE_avg = ΔE_avg / 1000

    return -(ΔE_avg / log(p0))
end

function report(sa::SA)
    println("- step $(sa.step), T = $(sa.T[end]), E = $(sa.E[end]), Ebsf = $(sa.Ebsf[end])")
end

function prun(;sas::Vector{SA{I,J}} where I where J,
              s0s,
              neighbor::Function,
              energy::Function,
              nprocs::Int,
              T0::Float64 = -1.0,
              revert::Union{Function, Nothing} = nothing)

    #awkward indexed loop, but @threads doesn't support i,j iteration over zips
    Threads.@threads for i in 1:length(sas)
        run(sa = sas[i], s0 = s0s[i], neighbor = neighbor, energy = energy,
            T0 = T0, revert = revert)
    end
    return sas

end

function run(;sa::SA{ET,ST} where ET where ST,
             s0,
             neighbor::Function,
             energy::Function,
             T0::Float64 = -1.0,
             revert::Union{Function, Nothing} = nothing)

    # either guess initial temperature, or user-given T0
    # error-out if T0 not given and no guess
    if sa.guess_T0
        neighbor_fn = sa.fused_neighbor_energy ? x -> neighbor(x)[1] : neighbor
        T = initial_temperature(s0, energy, neighbor_fn, sa.p0)
    elseif T0 < 0.0
        error("if guess_initial_temperature = false, a T0 > 0.0 must be set")
    else
        T = T0
    end

    # initialize best, current solutions, energies and number of moves
    E = energy(s0)
    zero_E = zero(E)
    s = s0
    sa.step = 0
    s_best = s
    Ebsf = E

    # allocate energy vectors once and reuse for every annealing step
    if sa.energy_statistics
        E_candidates = zeros(typeof(E), sa.n_max_moves)
        E_accepted  = zeros(typeof(E), sa.n_max_accepted)
    end

    # initialize copy of current solution
    scopy = copy(s0)

    while !sa.done(sa)
        n_accepted = 0
        n_moves = 0

        for i in 1:sa.n_max_moves

            if sa.fused_neighbor_energy
                trial,ΔE = neighbor(scopy)
                E_trial = E + ΔE
            else
                trial = neighbor(scopy)
                E_trial = energy(trial)
                ΔE = E_trial - E
            end

            n_moves += 1

            if sa.energy_statistics
                E_candidates[n_moves] = E_trial
            end

            # accept a neighbor with lower energy or with P = e^(-ΔE/T)
            if ΔE <= zero_E ||  @fastmath rand() < exp(-ΔE / T)
                s = trial
                E = E_trial
                n_accepted += 1

                if sa.energy_statistics
                    E_accepted[n_accepted] = E_trial
                end

                scopy = copy(s)
            else
                scopy = revert != nothing ? revert(scopy, s) : copy(s)
            end

            # update best solution/energy
            if E_trial < Ebsf
                Ebsf = E_trial
                s_best = trial
            end

            # max number of accepted solutions reached? -> next annealing step
            if n_accepted == sa.n_max_accepted
                break
            end
        end

        if sa.energy_statistics
            push!(sa.mean_E_total, Statistics.mean(E_candidates[1:n_moves]))
            push!(sa.mean_E_accepted, Statistics.mean(E_accepted[1:n_accepted]))
            push!(sa.var_E_total, Statistics.var(E_candidates[1:n_moves]))
            push!(sa.var_E_accepted, Statistics.var(E_accepted[1:n_accepted]))
        end

        if sa.solution_history
            push!(sa.s, s)
            push!(sa.s_best, s_best)
        end

        push!(sa.n_moves, n_moves)
        push!(sa.n_accepted, n_accepted)
        push!(sa.E, E)
        push!(sa.Ebsf, Ebsf)

        sa.step += 1
        push!(sa.T,T)
        sa.current_solution = s
        sa.bsf = s_best

        if sa.report_every > 0 && sa.step % sa.report_every == 0
            report(sa)
        end

        T = sa.schedule(sa)
    end

    if sa.report_every >= 0
        report(sa)
    end

    return sa
end

end


text:
"""

Attention:
    EN COURS DE MISE AU POINT : NE PAS UTILISER : voir Pkg AnsiColor
    <http://ascii-table.com/ansi-escape-sequences.php>
    <http://misc.flogisoft.com/bash/tip_colors_and_formatting>
    Voir aussi ma commande bash : diam_color.sh

Example
  println(to_sc("== TEST to_sc pour BLUE_STABILO ==", :BLUE_STABILO))

Sur les versions récentes de julia on peut faire par exemple :
printstyled("Solution invalide...", color=:red)

Cependant :

- julia doit être lancé avec l'option --color=yes
- cette commande ne permet pas de créer un chaine pour l'afficher plus tard.

"""
function to_sc(str::AbstractString, color::Symbol=:STABILO)
    red="9";    red_l="210";  red_d="1"
    grn="2" ;   grn_l="46";   gre_d="22"
    blu="12" ;  grn_l="14" ;  blu_d="4"
    yel="226";  yel_l="228";  yel_d="3"
    ora="208";  ora_l="215";  ora_d="130"
    colors = Dict(
        :STABILO         => "\e1\e[38;5;$(red)m\e[48;5;$(yel)m",
        :RED_STABILO     => "\e1\e[38;5;$(red)m\e[48;5;$(yel)m",
        :GREEN_STABILO   => "\e1\e[38;5;$(grn)m\e[48;5;$(yel)m",
        :BLUE_STABILO    => "\e1\e[38;5;$(blu)m\e[48;5;$(yel)m",
        :ORANGE_STABILO  => "\e1\e[38;5;$(ora_d)m\e[48;5;$(yel)m",
        :RED             => "\e1\e[38;5;$(red)m",
        :GREEN           => "\e1\e[38;5;$(grn)m",
        :BLUE            => "\e1\e[38;5;$(blu)m",
        :YELLOW          => "\e1\e[38;5;$(yel)m",
        :ORANGE          => "\e1\e[38;5;$(ora_d)m",
    )
    if haskey(colors, color)
        col = colors[color]
    else
        error("to_s_color: couleur inconnue : $(color).")
    end
    def = "\e[0m"
    buf = IOBuffer()
    print(buf, "$(col)", str, "$(def)")
    # return takebuf_string(buf)
    return String(take!(buf))
end

function prt_c(str::AbstractString, color::Symbol=:STABILO)
    print(to_sc(str, color))
end
function prn_c(str::AbstractString, color::Symbol=:STABILO)
    println(to_sc(str, color))
end

"""
Renvoie une chaine composée du caractère unicode spécifié par un symbol
(perso et non standard)
"""
function get_unicode(sym::Symbol)

    SYMS = Dict(
        # ligne brisée jougnant le nord et l'est passant par le centre
        :NE   => "\U2517", # Nord East
        :NS   => "\U2503",
        :NW   => "\U251B",
        :ES   => "\U250F",
        :EW   => "\U2501",
        :SW   => "\U2513", # South West

        :PP   => "\U254B", # croix totale style plus "+" bold

        :ne   => "\U2514",
        :ns   => "\U2502", # barre verticale light
        :nw   => "\U2518",
        :es   => "\U250C",
        :ew   => "\U2500", # barre horizontale light
        :sw   => "\U2510",
        :pp   => "\U253C", # croix totale style plus "+" light

        :tn   => "\U2534", # Tee vers N
        :te   => "\U251C", # Tee vers E
        :ts   => "\U252C", # Tee vers S
        :tw   => "\U2524", # Tee vers W

        :cp   => "\U00B7", # point centré
        :xx   => "\U2573", # croix totale X
    )
    if haskey(SYMS, sym)
        return SYMS[sym]
    else
        error("get_unicode: symbol UNICODE inconnu : $(sym).")
    end
end