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

using Markdown
using InteractiveUtils

# ╔═╡ 1eb34f10-b4b3-11ec-19f5-c7f26108d374
begin
	using LinearAlgebra
	using StatsPlots
md"""
# HW 11

Eric Nguyen
"""
end

# ╔═╡ 6c4b3bf7-26e9-4d6f-9305-65fbc851818f
let
md"""
### Problem 1

Use the properties of determinants to prove that ``\underline{\underline{A}} \,\underline{\underline{B}} = \underline{\underline{0}}`` implies that either ``|\underline{\underline{A}}| = 0`` or ``|\underline{\underline{B}}| = 0``, or both, where ``\underline{\underline{0}}`` is the matrix of all zeros.
Demonstrate this with the matrices

$\underline{\underline{A}} = \begin{bmatrix} 1 & 1 \\ 2 & 2 \end{bmatrix} \qquad\text{and}\qquad \underline{\underline{B}} = \begin{bmatrix} a & b \\ -a & -b \end{bmatrix}$

where ``a`` and ``b`` can be any complex numbers.
Which matrix has zero determinant?
"""
let
	A = [1 1; 2 2]
	a = rand(Complex{Float64})
	b = rand(Complex{Float64})

	B = [a b; -a -b]
	det(A), det(B), eigvals(A), eigvals(B)
end
md"""
### Problem 1

Given that ``|\underline{\underline{A}} \,\underline{\underline{B}}| = 0`` and using the sixth property of determinants as listed in the course notes, ``|\underline{\underline{A}} \,\underline{\underline{B}}| = |\underline{\underline{A}}| |\underline{\underline{B}}|``, we get

$|\underline{\underline{A}} \,\underline{\underline{B}}| = |\underline{\underline{A}}| |\underline{\underline{B}}| = 0$

which implies that at least ``|\underline{\underline{A}}| = 0`` or ``|\underline{\underline{B}}| = 0``.
The determinant ``|\underline{\underline{A}}|`` is

$|\underline{\underline{A}}| = (1)(2) - (1)(2) = 2 - 2 = 0$

so ``\underline{\underline{A}}`` is the matrix with the zero determinant.
"""
end

# ╔═╡ 2743185f-b0d7-4195-8e8e-a1bfcf132934
let
md"""
### Problem 2

Follow the procedure we went through in class to find the symmetry axes of the conic section ``6x^2 + 12xy + y^2 = 16``, and find the angle they make with the ``x,y`` axes.
What kind of curve is this?
A plot would be helpful.
You can do this with MATHEMATICA if you want, but the necessary algebra is rather simple.
"""
md"""
### Problem 2

$\underline{\underline{A}} = \begin{bmatrix} 6 & 6 \\ 6 & 1 \end{bmatrix} \implies \underline{\underline{R}} = \frac{1}{\sqrt{13}} \begin{bmatrix} -2 & 3 \\ 3 & 2 \end{bmatrix} \implies \underline{\underline{R}}^T \underline{\underline{A}} \,\underline{\underline{R}} = \begin{bmatrix} -3 & 0 \\ 0 & 10 \end{bmatrix}$
"""
end

# ╔═╡ 982fab25-1414-4789-bc16-343ab29b4514
let
	A = [6 6; 6 1]
	v1 = eigvecs(A)[:,1]
	v2 = eigvecs(A)[:,2]
	R = [v1 v2]
	R' * A * R
end

# ╔═╡ 7e53c83c-ead7-484f-8762-4a0956ba7210
md"""
This produces the curve

$-3x^2 + 10y^2 = 16$

which is a hyperbola
"""

# ╔═╡ c8977146-1409-4eb1-bede-13436fd7214b
let
	plot(x -> sqrt(16 + 3x^2) / 10, color=:black, legend=false)
	plot!(x -> -sqrt(16 + 3x^2) / 10, color=:black, legend=false)
end

# ╔═╡ 201fb421-39b1-47db-8dec-609706461dca
let
md"""
### Problem 3

Find the eigenvalues and eigenvectors for the matrix

$\underline{\underline{\sigma}}_y = \begin{bmatrix} 0 & -1 \\ i & 0 \end{bmatrix}$

one of the three **Pauli matrices**.
Do this by hand, not with MATHEMATICA.
Normalize the eigenvectors and show that they are orthogonal.
"""
md"""
### Problem 3

First, we find the characteristic polynomial.

$\det(\underline{\underline{\sigma}}_y - \lambda \,\underline{\underline{I}}) = \begin{vmatrix} -\lambda & -i \\ i & -\lambda \end{vmatrix} = \lambda^2 + i^2 = 0 \implies \lambda = ±1$

The eigenvalues are

$\lambda^{(1)} = -1 \qquad \lambda^{(2)} = 1$

For ``\lambda^{(1)} = -1`` we have

$\begin{align*}
-v_1 - i v_2 &= 0 \\
i v_1 - v_2 &= 0
\end{align*}$

$\implies v_1 = 1, \, v_2 = i$

giving the normalized eigenvector

$\underline{v} = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ i \end{bmatrix}$

For ``\lambda^{(2)} = 1`` we have

$\begin{align*}
u_1 - i u_2 &= 0 \\
i u_1 - u_2 &= 0
\end{align*}$

$\implies u_1 = 1, \, u_2 = i$

giving the normalized eigenvector

$\underline{u} = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ -i \end{bmatrix}$

We can show that the eigenvectors are orthogonal:

$⟨v \vert u⟩ = v_i u_i^* = \left(\frac{1}{\sqrt{2}}\right)^2 + \left(\frac{i}{\sqrt{2}}\right)^2 = \frac{1}{\sqrt{2}} - \frac{1}{\sqrt{2}} = 0$
"""
end

# ╔═╡ e12e34c7-90db-4f31-8792-ce2a39acc435
let
md"""
### Problem 4

Find the eigenvalues, two of which equal each other, of the real symmetric matrix

$\underline{\underline{A}} = \begin{bmatrix} 13 & 4 & -2 \\ 4 & 13 & -2 \\ -2 & -2 & 10 \end{bmatrix}$

Construct the three eigenvectors by hand, not with MATHEMATICA.
You will have more freedom than you would have thought.
Do you see how to use this freedom to make all three eigenvectors mutually orthogonal?
"""
md"""
### Problem 4

We first calculate the characteristic equation

$\det(\underline{\underline{A}} - \lambda \, \underline{\underline{I}}) = \begin{vmatrix} 13 - \lambda & 4 & -2 \\ 4 & 13 - \lambda & -2 \\ -2 & -2 & 10 - \lambda \end{vmatrix} = -(\lambda - 18)(\lambda - 9)^2$

So the eigenvalues are

$\lambda^{(1)} = 9 \qquad \lambda^{(2)} = 18$

For ``\lambda^{(1)} = 9`` we have

$\begin{align*}
4v_1 + 4v_2 - 2v_3 &= 0 \\
4v_1 + 4v_2 - 2v_3 &= 0 \\
-2v_1 - 2v_2 + v_3 &= 0
\end{align*}$

$\implies v_1 = 1, v_2 = 0, v_3 = 2 \qquad \text{or} \qquad v_1 = -1, v_2 = 1, v_3 = 0$

which gives the normalized eigenvectors

$\underline{v} = \frac{1}{\sqrt{5}} \begin{bmatrix} 1 \\ 0 \\ 2 \end{bmatrix} \qquad \underline{u} = \frac{1}{\sqrt{2}} \begin{bmatrix} -1 \\ 1 \\ 0 \end{bmatrix}$

For ``\lambda^{(2)} = 18`` we have

$\begin{align*}
-5w_1 + 4w_2 - 2w_3 &= 0 \\
4w_1 - 5w_2 - 2w_3 &= 0 \\
-2w_1 - 2w_2 - 8w_3 &= 0
\end{align*}$

$\implies w_1 = -2, w_2 = -2, w_3 = 1$

which gives the normalized eigenvector

$\underline{w} = \frac{1}{3} \begin{bmatrix} -2 \\ -2 \\ 1 \end{bmatrix}$
"""
end

# ╔═╡ 1aaa8364-b72e-4128-b26a-82d4d4601634
let
md"""
### Problem 5

Find the eigenfrequencies and eigenmodes for the mechanical system

`|~~~k~~~[ m ]~~~k~~~[ m ]~~~k~~~[ m ]~~~k~~~|`

Make a plot that shows the motions of each of the three masses, for the three sets of initial conditions where the masses start at rest with position given by each of the three eigenvectors.
Briefly describe the motions of the three masses, for each of the eigenmodes.
"""
md"""
### Problem 5

The coupled differential equations can be written as

$\begin{align*}
mẍ_1 &= -kx_1 + k(x_2 - x_1) = -2kx_1 + kx_2 \\
mẍ_2 &= -k(x_2 - x_1) + k(x_3 - x_2) = kx_1 - 2kx_2 + kx_3 \\
mẍ_3 &= -k(x_3 - x_2) + k(-x_3) = kx_2 - 2kx_3
\end{align*}$

We can rewrite them using ``\omega_0 = \sqrt{k / m}``,

$\begin{align*}
ẍ_1 &= -2{\omega_0}^2 x_1 + {\omega_0}^2 x_2 \\
ẍ_2 &= {\omega_0}^2 x_1 - 2{\omega_0}^2 x_2 + {\omega_0}^2 x_3 \\
ẍ_3 &= {\omega_0}^2 x_2 - 2{\omega_0}^2 x_3
\end{align*}$

In terms of vectors and matrices, we can write this as

$\underline{ẍ}(t) = -{\omega_0}^2 \underline{\underline{Ω}} \,\underline{x}(t)$

where

$\underline{x} = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} \qquad Ω = \begin{bmatrix} -2 & 1 & 0 \\ 1 & -2 & 1 \\ 0 & 1 & -2 \end{bmatrix}$

with the help of MATHEMATICA, we find the eigenfrequencies to be

$\omega_1 = \sqrt{(2 - \sqrt{2}) \omega_0}, \, \omega_2 = \sqrt{2\omega_0}, \, \omega_3 = \sqrt{(2 + \sqrt{2}) \omega_0}$

and the normalized eigenmodes to be

$\underline{a}^{(1)} = \frac{1}{2} \begin{bmatrix} 1 \\ -\sqrt{2} \\ 1 \end{bmatrix} \qquad \underline{a}^{(2)} = \frac{1}{\sqrt{2}} \begin{bmatrix} -1 \\ 0 \\ 1 \end{bmatrix} \qquad \underline{a}^{(3)} = \frac{1}{2} \begin{bmatrix} 1 \\ \sqrt{2} \\ 1 \end{bmatrix}$
"""
end

# ╔═╡ 7387c353-51b7-402b-a500-adec0a34b059
let
	Ω = [-2 1 0; 1 -2 1; 0 1 -2]
	eigen(Ω)
end

# ╔═╡ b40838d0-8391-4f1e-a3fc-a4be028b7bc7
(2 - sqrt(2))

# ╔═╡ ec26fefb-18d1-4ac7-b0c1-0d3c71ae1305
begin
	ω₁ = sqrt(2 - sqrt(2))
	ω₂ = sqrt(2)
	ω₃ = sqrt(2 + sqrt(2))

	a1 = (1/2) .* [1; -sqrt(2); 1]
	a2 = (1/sqrt(2)) .* [-1; 0; 1]
	a3 = (1/2) .* [1; sqrt(2); 1]

	s1(t, A, ϕ) = A .* a1 .* cos(ω₁ * t + ϕ)
	s2(t, A, ϕ) = A .* a2 .* cos(ω₂ * t + ϕ)
	s3(t, A, ϕ) = A .* a3 .* cos(ω₃ * t + ϕ)
end;

# ╔═╡ 355540ec-db9b-4ce5-a5d2-0c10497da934
let
	x1(t) = s1(t, 1, 0)[1]
	x2(t) = s1(t, 1, 0)[2]
	x3(t) = s1(t, 1, 0)[3]
	plot(x1, label="x₁")
	plot!(x2, label="x₂")
	plot!(x3, label="x₃")
	title!("Eigenmode 1")
end

# ╔═╡ de44b0c7-eee4-49ed-9264-40a5339e7d6c
md"""
Eigenmode 1: x1 and x3 oscillate with each other, opposite to x2
"""

# ╔═╡ 6f073f29-39b5-496e-b644-45051c81ab54
let
	x1(t) = s2(t, 1, 0)[1]
	x2(t) = s2(t, 1, 0)[2]
	x3(t) = s2(t, 1, 0)[3]
	plot(x1, label="x₁")
	plot!(x2, label="x₂")
	plot!(x3, label="x₃")
	title!("Eigenmode 2")
end

# ╔═╡ 5500a012-48fd-4227-a70f-515425c5122a
md"""
Eigenmode 2: x1 and x3 oscillate against each other while x2 remains in place
"""

# ╔═╡ e0a8d8ea-e5bd-48ca-b44f-4b6625181a17
let
	x1(t) = s3(t, 1, 0)[1]
	x2(t) = s3(t, 1, 0)[2]
	x3(t) = s3(t, 1, 0)[3]
	plot(x1, label="x₁")
	plot!(x2, label="x₂")
	plot!(x3, label="x₃")
	title!("Eigenmode 3")
end

# ╔═╡ cc32c2a5-d377-4a77-8477-b3acc63cca3b
md"""
Eigenmode 3: x1, x2, and x3 all oscillate with each other, however x2 oscillation is amplified
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"

[compat]
StatsPlots = "~0.14.33"
"""

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

julia_version = "1.8.0-beta1"
manifest_format = "2.0"
project_hash = "94b1ff0197986fb39c41ec4fbe1e3c63bdbc66ff"

[[deps.AbstractFFTs]]
deps = ["ChainRulesCore", "LinearAlgebra"]
git-tree-sha1 = "6f1d9bc1c08f9f4a8fa92e3ea3cb50153a1b40d4"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.1.0"

[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.3"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
version = "1.1.1"

[[deps.Arpack]]
deps = ["Arpack_jll", "Libdl", "LinearAlgebra", "Logging"]
git-tree-sha1 = "91ca22c4b8437da89b030f08d71db55a379ce958"
uuid = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
version = "0.5.3"

[[deps.Arpack_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "OpenBLAS_jll", "Pkg"]
git-tree-sha1 = "5ba6c757e8feccf03a1554dfaf3e26b3cfc7fd5e"
uuid = "68821587-b530-5797-8361-c406ea357684"
version = "3.5.1+1"

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

[[deps.AxisAlgorithms]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"]
git-tree-sha1 = "66771c8d21c8ff5e3a93379480a2307ac36863f7"
uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950"
version = "1.0.1"

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

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

[[deps.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 = "4b859a208b2397a7a623a03449e4636bdb17bcf2"
uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
version = "1.16.1+1"

[[deps.Calculus]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad"
uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
version = "0.5.1"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.14.0"

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

[[deps.Clustering]]
deps = ["Distances", "LinearAlgebra", "NearestNeighbors", "Printf", "SparseArrays", "Statistics", "StatsBase"]
git-tree-sha1 = "75479b7df4167267d75294d14b58244695beb2ac"
uuid = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5"
version = "0.14.2"

[[deps.ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "12fc73e5e0af68ad3137b886e3f7c1eacfca2640"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.17.1"

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

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

[[deps.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 = "96b0bc6c52df76506efc8a441c6cf1adcb1babc4"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.42.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
version = "0.5.0+0"

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

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

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

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

[[deps.DataValues]]
deps = ["DataValueInterfaces", "Dates"]
git-tree-sha1 = "d88a19299eba280a6d062e135a43f00323ae70bf"
uuid = "e7dc6d0d-1eca-5fa6-8ad6-5aecde8b7ea5"
version = "0.4.13"

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

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

[[deps.DensityInterface]]
deps = ["InverseFunctions", "Test"]
git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
version = "0.4.0"

[[deps.Distances]]
deps = ["LinearAlgebra", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "3258d0659f812acde79e8a74b11f17ac06d0ca04"
uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
version = "0.10.7"

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

[[deps.Distributions]]
deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
git-tree-sha1 = "5a4168170ede913a2cd679e53c2123cb4b889795"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.53"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

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

[[deps.DualNumbers]]
deps = ["Calculus", "NaNMath", "SpecialFunctions"]
git-tree-sha1 = "5837a837389fccf076445fce071c8ddaea35a566"
uuid = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
version = "0.6.8"

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

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

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

[[deps.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"

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

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

[[deps.FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[deps.FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "246621d23d1f43e3b9c368bf3b72b2331a27c286"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.13.2"

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

[[deps.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"

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

[[deps.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"

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

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

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

[[deps.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 = "cd6efcf9dc746b06709df14e462f0a3fe0786b1e"
uuid = "d2c73de3-f751-5644-a686-071e5b155ba9"
version = "0.64.2+0"

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

[[deps.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"

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

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

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

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

[[deps.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"

[[deps.HypergeometricFunctions]]
deps = ["DualNumbers", "LinearAlgebra", "SpecialFunctions", "Test"]
git-tree-sha1 = "65e4589030ef3c44d3b90bdc5aac462b4bb05567"
uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
version = "0.3.8"

[[deps.IniFile]]
git-tree-sha1 = "f550e6e32074c939295eb5ea6de31849ac2c9625"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.1"

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

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

[[deps.Interpolations]]
deps = ["AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"]
git-tree-sha1 = "b15fc0a95c564ca2e0a7ae12c1f095ca848ceb31"
uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
version = "0.13.5"

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "91b5dcf362c5add98049e6c29ee756910b03051d"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.3"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.IterTools]]
git-tree-sha1 = "fa6287a4469f5e048d763df38279ee729fbd44e5"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.4.0"

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

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

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

[[deps.KernelDensity]]
deps = ["Distributions", "DocStringExtensions", "FFTW", "Interpolations", "StatsBase"]
git-tree-sha1 = "591e8dc09ad18386189610acafb970032c519707"
uuid = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
version = "0.6.3"

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

[[deps.LERC_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "bf36f528eec6634efc60d7ec062008f171071434"
uuid = "88015f11-f218-50d7-93a8-a6af411a945d"
version = "3.0.0+1"

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

[[deps.LaTeXStrings]]
git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.3.0"

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

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

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

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

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

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

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

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

[[deps.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"

[[deps.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"

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

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

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

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

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

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

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "58f25e56b706f95125dcb796f39e1fb01d913a71"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.10"

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

[[deps.MKL_jll]]
deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
git-tree-sha1 = "e595b205efd49508358f7dc670a940c790204629"
uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
version = "2022.0.0+0"

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

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

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

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
version = "2.28.0+0"

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

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

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

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
version = "2022.2.1"

[[deps.MultivariateStats]]
deps = ["Arpack", "LinearAlgebra", "SparseArrays", "Statistics", "StatsAPI", "StatsBase"]
git-tree-sha1 = "7008a3412d823e29d370ddc77411d593bd8a3d03"
uuid = "6f286f6a-111f-5878-ab1e-185364afe411"
version = "0.9.1"

[[deps.NaNMath]]
git-tree-sha1 = "737a5957f387b17e74d4ad2f440eb330b39a62c5"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "1.0.0"

[[deps.NearestNeighbors]]
deps = ["Distances", "StaticArrays"]
git-tree-sha1 = "ded92de95031d4a8c61dfb6ba9adb6f1d8016ddd"
uuid = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
version = "0.4.10"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
version = "1.2.0"

[[deps.Observables]]
git-tree-sha1 = "fe29afdef3d0c4a8286128d4e45cc50621b1e43d"
uuid = "510215fc-4207-5dde-b226-833fc4488ee2"
version = "0.4.0"

[[deps.OffsetArrays]]
deps = ["Adapt"]
git-tree-sha1 = "043017e0bdeff61cfbb7afeb558ab29536bbb5ed"
uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
version = "1.10.8"

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

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
version = "0.3.17+2"

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
version = "0.8.1+0"

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

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

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

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

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

[[deps.PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "e8185b83b9fc56eb6456200e873ce598ebc7f262"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.7"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "621f4f3b4977325b9128d5fae7a8b4829a0c2222"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.4"

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

[[deps.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"
version = "1.8.0"

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

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

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

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "d3538e7f8a790dc8903519090857ef8e1283eecd"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.5"

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

[[deps.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"

[[deps.QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.2"

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

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

[[deps.Ratios]]
deps = ["Requires"]
git-tree-sha1 = "dc84268fe0e3335a62e315a3a7cf2afa7178a734"
uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439"
version = "0.4.3"

[[deps.RecipesBase]]
git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.1"

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

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

[[deps.RelocatableFolders]]
deps = ["SHA", "Scratch"]
git-tree-sha1 = "cdbd3b1338c72ce29d9584fdbe9e9b70eeb5adca"
uuid = "05181044-ff0b-4ac5-8273-598c1e38db00"
version = "0.1.3"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[deps.Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
version = "0.7.0"

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

[[deps.SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "6a2f7d70512d205ca8c7ee31bfa9f142fe74310c"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.12"

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

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

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

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

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

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

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.1.4"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "4f6ec5d99a28e1a749559ef7dd518663c5eca3d5"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.4.3"

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

[[deps.StatsAPI]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "c3d8ba7f3fa0625b062b82853a7d5229cb728b6b"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.2.1"

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

[[deps.StatsFuns]]
deps = ["ChainRulesCore", "HypergeometricFunctions", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "72e6abd6fc9ef0fa62a159713c83b7637a14b2b8"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.17"

[[deps.StatsPlots]]
deps = ["AbstractFFTs", "Clustering", "DataStructures", "DataValues", "Distributions", "Interpolations", "KernelDensity", "LinearAlgebra", "MultivariateStats", "Observables", "Plots", "RecipesBase", "RecipesPipeline", "Reexport", "StatsBase", "TableOperations", "Tables", "Widgets"]
git-tree-sha1 = "4d9c69d65f1b270ad092de0abe13e859b8c55cad"
uuid = "f3b207a7-027a-5e70-b257-86293d7955fd"
version = "0.14.33"

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

[[deps.SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

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

[[deps.TableOperations]]
deps = ["SentinelArrays", "Tables", "Test"]
git-tree-sha1 = "e383c87cf2a1dc41fa30c093b2a19877c83e1bc1"
uuid = "ab02a1b2-a7df-11e8-156e-fb1833f50b87"
version = "1.2.0"

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

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"]
git-tree-sha1 = "5ce79ce186cc678bbb5c5681ca3379d1ddae11a1"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.7.0"

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

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

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

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

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

[[deps.UnicodeFun]]
deps = ["REPL"]
git-tree-sha1 = "53915e50200959667e78a92a418594b428dffddf"
uuid = "1cfade01-22cf-5700-b092-accc4b62d6e1"
version = "0.4.1"

[[deps.Unzip]]
git-tree-sha1 = "34db80951901073501137bdbc3d5a8e7bbd06670"
uuid = "41fe7b60-77ed-43a1-b4f0-825fd5a5650d"
version = "0.1.2"

[[deps.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"

[[deps.Wayland_protocols_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4528479aa01ee1b3b4cd0e6faef0e04cf16466da"
uuid = "2381bf8a-dfd0-557d-9999-79630e7b1b91"
version = "1.25.0+0"

[[deps.Widgets]]
deps = ["Colors", "Dates", "Observables", "OrderedCollections"]
git-tree-sha1 = "505c31f585405fc375d99d02588f6ceaba791241"
uuid = "cc8bc4a8-27d6-5769-a93b-9d913e69aa62"
version = "0.6.5"

[[deps.WoodburyMatrices]]
deps = ["LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "de67fa59e33ad156a590055375a30b23c40299d3"
uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6"
version = "0.5.5"

[[deps.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"

[[deps.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"

[[deps.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"

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

[[deps.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"

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

[[deps.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"

[[deps.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"

[[deps.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"

[[deps.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"

[[deps.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"

[[deps.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"

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

[[deps.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"

[[deps.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"

[[deps.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"

[[deps.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"

[[deps.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"

[[deps.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"

[[deps.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"

[[deps.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"

[[deps.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"

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

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

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

[[deps.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"

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

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

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

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

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
version = "1.41.0+1"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
version = "16.2.1+1"

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

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

[[deps.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:
# ╟─1eb34f10-b4b3-11ec-19f5-c7f26108d374
# ╟─6c4b3bf7-26e9-4d6f-9305-65fbc851818f
# ╟─2743185f-b0d7-4195-8e8e-a1bfcf132934
# ╠═982fab25-1414-4789-bc16-343ab29b4514
# ╟─7e53c83c-ead7-484f-8762-4a0956ba7210
# ╟─c8977146-1409-4eb1-bede-13436fd7214b
# ╟─201fb421-39b1-47db-8dec-609706461dca
# ╟─e12e34c7-90db-4f31-8792-ce2a39acc435
# ╟─1aaa8364-b72e-4128-b26a-82d4d4601634
# ╠═7387c353-51b7-402b-a500-adec0a34b059
# ╠═b40838d0-8391-4f1e-a3fc-a4be028b7bc7
# ╟─ec26fefb-18d1-4ac7-b0c1-0d3c71ae1305
# ╟─355540ec-db9b-4ce5-a5d2-0c10497da934
# ╟─de44b0c7-eee4-49ed-9264-40a5339e7d6c
# ╟─6f073f29-39b5-496e-b644-45051c81ab54
# ╟─5500a012-48fd-4227-a70f-515425c5122a
# ╟─e0a8d8ea-e5bd-48ca-b44f-4b6625181a17
# ╟─cc32c2a5-d377-4a77-8477-b3acc63cca3b
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
module Mods

import Base: isequal, (==), (+), (-), (*), (inv), (/), (//), (^), hash, show
import Base: zero, one, rand, conj

export Mod, modulus, value, AbstractMod
export isequal, ==, +, -, *, is_invertible, inv, /, ^
export hash, CRT

abstract type AbstractMod <: Number end

"""
`Mod{m}(v)` creates a modular number in mod `m` with value `mod(v,m)`.
`Mod{m}()` is equivalent to `Mod{m}(0)`.
"""
struct Mod{N} <: AbstractMod
    val::Int
    function Mod(x::Int,N::Int)
        @assert N>1 "modulus must be at least 2"
        new{N::Int}(mod(x,N))
    end
end

Mod{N}(x::Int=0) where N = Mod(x,N)

function Mod{N}(x::T) where {N,T<:Integer}
    xx = mod(x,N)
    return Mod{N}(Int(xx))
end
Mod(x::T, N::Int) where T<:Integer = Mod{N}(x)

"""
`modulus(a::Mod)` returns the modulus of this `Mod` number.
```
julia> a = Mod{13}(11);

julia> modulus(a)
13
```
"""
modulus(a::Mod{N}) where N = N

"""
`value(a::Mod)` returns the value of this `Mod` number.
```
julia> a = Mod{13}(11);

julia> value(a)
11
```
"""
value(a::AbstractMod) = a.val

zero(::Mod{N}) where N = Mod{N}()
zero(::Type{Mod{N}}) where N = Mod{N}()

one(::Mod{N}) where N = Mod{N}(1)
one(::Type{Mod{N}}) where N = Mod{N}(1)

conj(x::Mod) = x   # so matrix transpose will work

function hash(x::AbstractMod, h::UInt64= UInt64(0))
    v = value(x)
    m = modulus(x)
    return hash(v,hash(m,h))
end

# Test for equality
isequal(x::Mod, y::Mod) = value(x)==value(y) && modulus(x)==modulus(y)

==(x::Mod,y::Mod) = isequal(x,y)

"""
`modcheck(x::Mod,y::Mod)` checks if `x` and `y` have the same modulus.
If so, return than modulus. If not, throw an error.
"""
function modcheck(x::Mod, y::Mod)::Int
    mx = modulus(x)
    my = modulus(y)
    if mx != my
        error("Cannot operate on two Mod objects with different moduli")
    end
    return mx
end

# Easy arithmetic
function +(x::Mod, y::Mod)::Mod
    m = modcheck(x,y)
    s,flag = Base.add_with_overflow(x.val,y.val)
    if !flag
        return Mod(x.val+y.val, m)
    end
    s = widen(x.val) + widen(y.val)    # add with added precision
    s = Int(mod(s,m))                  # reduce by modulus
    return Mod(s,m)
end

function -(x::Mod{M}) where M
    return Mod(-x.val, M)
end

-(x::Mod,y::Mod) = x + (-y)

function *(x::Mod, y::Mod)
    m = modcheck(x,y)
    p,flag = Base.mul_with_overflow(x.val,y.val)
    if !flag
        return Mod(x.val*y.val, m)
    end
    p = widemul(x.val, y.val)         # multipy with added precision
    p = Int(mod(p,m))                  # reduce by the modulus
    return Mod(p,m) # return with proper type
end

# Division stuff
"""
`is_invertible(x::Mod)` determines if `x` is invertible.
"""
function is_invertible(x::Mod{M})::Bool where M
    return gcd(value(x),M) == 1
end

"""
`inv(x::Mod)` gives the multiplicative inverse of `x`.
This may be abbreviated by `x'`.
"""
function inv(x::Mod{M}) where M
    (g, v, ignore) = gcdx(x.val, M)
    if g != 1
        error("$x is not invertible")
    end
    return Mod(v, M)
end

function /(x::Mod, y::Mod)
    modcheck(x,y)
    return x * inv(y)
end

(//)(x::Mod,y::Mod) = x/y

# Operations with Integers

(+)(x::Mod{M}, k::Integer) where M = Mod(k,M)+x
(+)(k::Integer, x::Mod) = x+k

(-)(x::Mod, k::Integer) = x + (-k)
(-)(k::Integer, x::Mod) = (-x) + k

(*)(x::Mod{M}, k::Integer) where M = Mod(k,M) * x
(*)(k::Integer, x::Mod) = x*k

(/)(x::Mod{M}, k::Integer) where M = x / Mod(k, M)
(/)(k::Integer, x::Mod{M}) where M = Mod(k, M) / x

(//)(x::Mod{M}, k::Integer) where M = x / Mod(k, M)
(//)(k::Integer, x::Mod{M}) where M = Mod(k, M) / x

# Operations with rational numbers

Mod{N}(k::Rational) where N = Mod{N}(numerator(k))/Mod{N}(denominator(k))

function +(x::Mod{N}, k::Rational) where N
    return x + Mod{N}(k)
end
(+)(k::Rational,x::Mod) = x+k

(-)(x::Mod,k::Rational) = x + (-k)
(-)(k::Rational,x::Mod) = k + (-x)

function (*)(x::Mod{N},k::Rational) where N
    return x * Mod{N}(k)
end
(*)(k::Rational,x::Mod) = x*k

function (/)(x::Mod,k::Rational)
    return x * (1/k)
end
(/)(k::Rational,x::Mod) = k * inv(x)

(//)(x::Mod,k::Rational) = x/k
(//)(k::Rational,x::Mod) = k/x

# Comparison with Integers

isequal(x::Mod{M}, k::Integer) where M = mod(k,M) == x.val
isequal(k::Integer, x::Mod) = isequal(x,k)
(==)(x::Mod, k::Integer) = isequal(x,k)
(==)(k::Integer, x::Mod) = isequal(x,k)

# Comparisons with Rationals
function isequal(x::Mod{N}, k::Rational) where N
    return x == Mod{N}(k)
end
isequal(k::Rational,x::Mod) = isequal(x,k)
(==)(x::Mod, k::Rational) = isequal(x,k)
(==)(k::Rational, x::Mod) = isequal(x,k)

# Random

rand(::Type{Mod{N}}) where N = Mod{N}(rand(Int))
rand(::Type{Mod{N}},dims::Integer...) where N = Mod{N}.(rand(Int,dims...))

# Chinese remainder theorem functions

# private helper function
function CRT_work(x::Mod{n}, y::Mod{m}) where {n,m}
    # n = x.mod
    # m = y.mod
    if gcd(n,m) != 1
        error("Moduli must be pairwise relatively prime")
    end

    a = x.val
    b = y.val

    k = inv(Mod(n,m)) * (b-a)

    z = a + k.val*n

    return Mod(z, n*m)
end

# public interface
"""
`CRT(m1,m2,...)`: Chinese Remainder Theorem
```
julia> CRT( Mod(4,11), Mod(8,14) )
Mods.Mod(92,154)

julia> 92%11
4

julia> 92%14
8
```
"""
function CRT(mtuple::Mod...)
    n = length(mtuple)
    if n == 0
        return 1
    end

    result::Mod = mtuple[1]

    for k=2:n
        result = CRT_work(result,mtuple[k])
    end

    return result
end

include("GaussMods.jl")

end # end of module Mods

ext: jl


text:
# Unit test for Dygraphs file reading
@testset "Dygraphs file reading" begin
	data = readdygraphs(joinpath(dirname(@__DIR__),"test","input","dygraphs_data.csv"),datestring="yyyymmdd");
	@test size(data) == (10,3)
	@test names(data)==[:datetime,:High,:Low]
	@test data[:datetime][1] == DateTime(2007,01,01,00,00,00)
	@test sum(data[:High]) == 600.
	@test data[:Low][end] == 37.
end

@testset "Dygraphs file write" begin
	data = DataFrame(temp=[10.,11.,12.,14.],grav=[9.8123,9.9,NaN,9.7],
	       datetime=[DateTime(2010,1,1,0),DateTime(2010,1,1,1),
	           DateTime(2010,1,1,2),DateTime(2010,1,1,4)]);
	writedygraphs(data,joinpath(dirname(@__DIR__),"test","output","dygraphs_data.csv"),decimal=[1,3]);
	data_read = readdlm(joinpath(dirname(@__DIR__),"test","output","dygraphs_data.csv"),',',skipstart=1);
	@test size(data_read) == (4,3)
	@test data_read[1,1] == "2010/01/01 00:00:00"
	@test data_read[2,end] == 9.9
	@test data_read[end,2] == 14
end


text:
gaussjacobi(n::Number, a::Number, b::Number) =
    gaussjacobi(Int(n),Float64(a),Float64(b))

function gaussjacobi(n::Integer, a::Float64, b::Float64)
    #GAUSS-JACOBI QUADRATURE NODES AND WEIGHTS
    if n < 0
        error("gaussjacobi($n,$a,$b) not defined: n must be positive.")
    elseif a == 0. && b == 0.
        gausslegendre(n)
    elseif a == -0.5 && b == -0.5
        gausschebyshev(n, 1)
    elseif a == 0.5 && b == 0.5
        gausschebyshev(n, 2)
    elseif a == -0.5 && b == 0.5
        gausschebyshev(n, 3)
    elseif a == 0.5 && b == -0.5
        gausschebyshev(n, 4)
    elseif n == 0
        Float64[], Float64[]
    elseif n == 1
        [(b - a) / (a + b + 2)], [2^(a + b + 1) * beta(a + 1, b + 1)]
    elseif min(a,b) <= -1.
        error("The Jacobi parameters correspond to a nonintegrable weight function")
    elseif n <= 100 && max(a,b) < 5.
        JacobiRec(n, a, b)
    elseif n > 100 && max(a,b) < 5.
        JacobiAsy(n, a, b)
    elseif n <= 4000 && max(a,b)>=5.
        JacobiGW(n, a, b)
    else
        error("gaussjacobi($n,$a,$b) is not implemented: n must be ≤ 4000 for max(a,b)≥5.")
    end
end

# Convenience function: convert any kind of numbers a and b to a joint floating point type
JacobiRec(n::Integer, a::Number, b::Number) = JacobiRec(n, promote(float(a), float(b))...)

function JacobiRec(n::Integer, a::T, b::T) where {T <: AbstractFloat}
    #Compute nodes and weights using recurrrence relation.
    x11, x12 = HalfRec(n, a, b, 1)
    x21, x22 = HalfRec(n, b, a, 0)

    x = Array{T}(undef,n)
    w = Array{T}(undef,n)
    m1 = length(x11)
    m2 = length(x21)
    sum_w = zero(T)
    @inbounds for i in 1:m2
        idx = m2 + 1 - i
        xi = -x21[i]
        der = x22[i]
        wi = 1 / ((1 - xi^2) * der^2)
        w[idx] = wi
        x[idx] = xi
        sum_w += wi
    end
    @inbounds for i in 1:m1
        idx = m2 + i
        xi = x11[i]
        der = x12[i]
        wi = 1 / ((1 - xi^2) * der^2)
        w[idx] = wi
        x[idx] = xi
        sum_w += wi
    end
    c = (2^(a + b + 1) * gamma(2 + a) *
         gamma(2 + b) / (gamma(2 + a + b) * (a + 1) * (b + 1)))
    rmul!(w, c / sum_w)
    x, w
end

function HalfRec(n::Integer, a::T, b::T, flag) where {T <: AbstractFloat}
    # HALFREC  Jacobi polynomial recurrence relation.
    # Asymptotic formula - only valid for positive x.
    r = (flag == 1) ? (ceil(n / 2):-1:1) : (floor(n / 2):-1:1)
    m = length(r)
    c1 = 1 / (2 * n + a + b + 1)
    a1 = one(T)/4 - a^2
    b1 = one(T)/4 - b^2
    c1² = c1^2
    x = Array{T}(undef,m)
    @inbounds for i in 1:m
        C = muladd(2*one(T), r[i], a - one(T)/2) * (T(π) * c1)
        C_2 = C/2
        x[i] = cos(muladd(c1², muladd(-b1, tan(C_2), a1 * cot(C_2)), C))
    end

    P1 = Array{T}(undef,m)
    P2 = Array{T}(undef,m)
    # Loop until convergence:
    for _ in 1:10
        innerJacobiRec!(n, x, a, b, P1, P2)
        dx2 = 0.0
        @inbounds for i in 1:m
            dx = P1[i] / P2[i]
            _dx2 = abs2(dx)
            dx2 = ifelse(_dx2 > dx2, _dx2, dx2)
            x[i] = x[i] - dx
        end
        dx2 > eps(T) / 1e6 || break
    end
    # Once more for derivatives:
    innerJacobiRec!(n, x, a, b, P1, P2)
    x, P2
end

function innerJacobiRec!(n, x, a::T, b::T, P, PP) where {T <: AbstractFloat}
    # EVALUATE JACOBI POLYNOMIALS AND ITS DERIVATIVE USING THREE-TERM RECURRENCE.
    N = length(x)
    @inbounds for j = 1:N
        xj = x[j]

        Pj = (a - b + (a + b + 2) * xj)/2
        Pm1 = one(T)
        PPj = (a + b + 2)/2
        PPm1 = zero(T)
        for k = 1:n-1
            k0 = muladd(2*one(T), k, a + b)
            k1 = k0 + 1
            k2 = k0 + 2
            A = 2 * (k + 1) * (k + (a + b + 1)) * k0
            B = k1 * (a^2 - b^2)
            C = k0 * k1 * k2
            D = 2 * (k + a) * (k + b) * k2
            c1 = muladd(C, xj, B)
            Pm1, Pj = Pj, muladd(-D, Pm1, c1 * Pj) / A
            PPm1, PPj = PPj, muladd(c1, PPj, muladd(-D, PPm1, C * Pm1)) / A
        end
        P[j] = Pj
        PP[j] = PPj
    end
    nothing
end

function innerJacobiRec(n, x, a::T, b::T) where {T <: AbstractFloat}
    # EVALUATE JACOBI POLYNOMIALS AND ITS DERIVATIVE USING THREE-TERM RECURRENCE.
    N = length(x)
    P = Array{T}(undef,N)
    PP = Array{T}(undef,N)
    innerJacobiRec!(n, x, a, b, P, PP)
    P, PP
end

function weightsConstant(n, a, b)
    # Compute the constant for weights:
    M = min(20, n - 1)
    C = 1.0
    p = -a * b / n
    for m = 1:M
        C += p
        p *= -(m + a) * (m + b) / (m + 1) / (n - m)
        abs(p / C) < eps(Float64) / 100 && break
    end
    2^(a + b + 1) * C
end

function JacobiAsy(n, a, b)
    # ASY  Compute nodes and weights using asymptotic formulae.

    # Determine switch between interior and boundary regions:
    nbdy = 10
    bdyidx1 = n - (nbdy - 1):n
    bdyidx2 = nbdy:-1:1

    # Interior formula:
    x, w = asy1(n, a, b, nbdy)

    # Boundary formula (right):
    xbdy = boundary(n, a, b, nbdy)
    x[bdyidx1], w[bdyidx1] = xbdy

    # Boundary formula (left):
    if a != b
        xbdy = boundary(n, b, a, nbdy)
    end
    x[bdyidx2] = -xbdy[1]
    w[bdyidx2] = xbdy[2]

    rmul!(w, weightsConstant(n, a, b))
    x, w
end

function asy1(n::Integer, a::Float64, b::Float64, nbdy)
    # Algorithm for computing nodes and weights in the interior.

    # Approximate roots via asymptotic formula: (Gatteschi and Pittaluga, 1985)
    K = (2*(n:-1:1).+a.-0.5).*pi./(2*n+a+b+1)
    tt = K .+ (1/(2*n+a+b+1)^2).*((0.25-a^2).*cot.(0.5 .* K).-(0.25-b^2).*tan.(0.5*K))

    # First half (x > 0):
    t = tt[tt .<= pi/2]'
    mint = t[end-nbdy+1]
    idx = 1:max(findfirst(vec(t) .< mint)-1, 1)

    dt = 1.0
    counter = 0
    # Newton iteration
    while ( norm(dt,Inf) > sqrt(eps(Float64))/100 && counter < 10 )
        vals = feval_asy1(n, a, b, t, idx)         # Evaluate
        dt = vals[1]./vals[2]                                 # Newton update
        t += dt                                     # Next iterate
        counter += 1
        dt = dt[idx]
    end
    vals = feval_asy1(n, a, b, t, idx)      # Once more for luck
    t += vals[1]./vals[2]                                 # Newton update.

    # Store:
    x = cos.(t)
    w = 1 ./ vals[2].^2

    # Second half (x < 0):
    tmp = a; a = b; b = tmp
    t = pi .- tt[1:(n-length(x))]'
    mint = t[nbdy]
    idx = max(findfirst(vec(t) .> mint), 1):length(t)

    dt = 1.0; counter = 0;
    # Newton iteration
    while ( norm(dt,Inf) > sqrt(eps(Float64))/100 && counter < 10 )
        vals = feval_asy1(n, a, b, t, idx)  # Evaluate.
        dt = vals[1]./vals[2]                                # Newton update.
        t += dt                                     # Next iterate.
        counter += 1
        dt = dt[idx]
    end
    vals = feval_asy1(n, a, b, t, idx)     # Once more for luck.
    t += vals[1]./vals[2]                                 # Newton update.

    # Store:
    [(-).(cos.(vec(t)));vec(x)],[1 ./vec(vals[2]).^2;vec(w)]
end

function feval_asy1(n::Integer, a::Float64, b::Float64, t, idx)
    # Evaluate the interior asymptotic formula at x = cos(t).

    # Number of terms in the expansion:
    M = 20

    # Some often used vectors/matrices:
    onesT = fill(1.0,length(t))'; onesM = fill(1.0,M); MM = collect(0:M-1);

    # The sine and cosine terms:
    alpha = @. (0.5*(2*n+a+b+1+MM))*onesT * (onesM*t) - 0.5*(a+0.5)*pi
    cosA = cos.(alpha); sinA = sin.(alpha)

    sinT = onesM*sin.(t)
    cosT = onesM*cos.(t)
    cosA2 = cosA.*cosT .+ sinA.*sinT
    sinA2 = sinA.*cosT .- cosA.*sinT

    one = fill!(similar(t),1.0)
    sinT = vcat( one , cumprod(onesM[2:end].*(0.5.*csc.(0.5.*t))))
    cosT = 0.5.*sec.(0.5.*t)

    j = transpose(0:M-2)
    vec = (0.5.+a.+j).*(0.5.-a.+j)./(j.+1)./(2*n.+a.+b.+j.+2)
    P1 = [1 cumprod(vec,2)]
    P1[3:4:end] = -P1[3:4:end]
    P1[4:4:end] = -P1[4:4:end]
    P2 = Matrix(1.0I, M, M)
    for l = 0:M-1
        j = transpose(0:(M-l-2))
        vec = @. (0.5+b+j)*(0.5-b+j)/(j+1)/(2*n+a+b+j+l+2)
        P2[(l+1).+(1:length(j)),l+1] = cumprod(vec,2)
    end
    PHI = repeat(P1,M,1).*P2

    j = transpose(0:M-2)
    vec = @. (0.5+a+j)*(0.5-a+j)/(j+1)/(2*(n-1)+a+b+j+2)
    P1 = [1 cumprod(vec,2)]
    P1[3:4:end] = -P1[3:4:end]
    P1[4:4:end] = -P1[4:4:end]
    P2 = Matrix(1.0I, M, M)
    for l = 0:M-1
        j = transpose(0:(M-l-2))
        vec = @. (0.5+b+j)*(0.5-b+j)/(j+1)/(2*(n-1)+a+b+j+l+2)
        P2[(l+1).+(1:length(j)),l+1] = cumprod(vec,2)
    end
    PHI2 = repeat(P1,M,1).*P2

    S = 0; S2 = 0;
    SC = sinT
    for m = 1:M
        l = 1:2:m
        phi = PHI[m, l]'
        dS1 = phi * SC[l, :] .* cosA[m, :]'
        phi2 = PHI2[m, l]'
        dS12 = phi2*SC[l, :] .* cosA2[m, :]'
        l = 2:2:m
        phi = PHI[m, l]'
        dS2 = phi * SC[l, :] .* sinA[m, :]'
        phi2 = PHI2[m, l]'
        dS22 = phi2 * SC[l, :] .* sinA2[m, :]'
        if m - 1 > 10 && norm(dS1[idx] + dS2[idx], Inf) < eps(Float64) / 100
            break
        end
        S = S .+ dS1 .+ dS2
        S2 = S2 .+ dS12 .+ dS22
        SC[1:m,:] = SC[1:m,:].*cosT
    end

    # Constant out the front:
    dsa = 0.5*(a^2)/n; dsb = 0.5*(b^2)/n; dsab = 0.25*(a+b)^2/n
    ds = dsa + dsb - dsab; s = ds; j = 1
    dsold = ds # to fix a = -b bug.
    while ( (abs(ds/s) + dsold) > eps(Float64)/10 )
        dsold = abs(ds/s)
        j += 1
        tmp = -(j-1)/(j+1)/n
        dsa = tmp*dsa*a
        dsb = tmp*dsb*b
        dsab = .5*tmp*dsab*(a+b)
        ds = dsa + dsb - dsab
        s = s + ds
    end
    p2 = exp(s)*sqrt(2*pi)*sqrt((n+a)*(n+b)/(2*n+a+b))/(2*n+a+b+1)
    g = [1, 1/12, 1/288, -139/51840, -571/2488320, 163879/209018880,
         5246819/75246796800, -534703531/902961561600,
         -4483131259/86684309913600, 432261921612371/514904800886784000]
    f(g,z) = sum(g.*[1;cumprod(fill(1.0,9)./z)])
    C = p2*(f(g,n+a)*f(g,n+b)/f(g,2*n+a+b))*2/pi
    C2 = C*(a+b+2*n).*(a+b+1+2*n)./(4*(a+n).*(b+n))

    vals = C*S
    S2 = C2*S2

    # Use relation for derivative:
    ders = (n.*(a.-b.-(2*n+a+b).*cos.(t)).*vals .+ 2 .*(n+a).*(n+b).*S2)./(2*n+a+b)./sin.(t)
    t .= abs.( t )
    denom = 1 ./ real.(sin.(t/2).^(a+0.5).*cos.(t ./ 2).^(b+0.5))
    vals = vals.*denom
    ders = ders.*denom

    return (vals, ders)
end

function boundary(n::Integer, a::Float64, b::Float64, npts)
# Algorithm for computing nodes and weights near the boundary.

    # Use Newton iterations to find the first few Bessel roots:
    smallK = min(30, npts)
    jk = besselroots(a, min(npts, smallK))

    # Approximate roots via asymptotic formula: (see Olver 1974)
    phik = jk/(n + .5*(a + b + 1))
    x = cos.( phik .+ ((a^2-0.25).*(1 .-phik.*cot.(phik))./(8*phik) .- 0.25.*(a^2-b^2).*tan.(0.5.*phik))./(n + 0.5*(a + b + 1))^2 )

    dx = 1.0; counter = 0;
    # Newton iteration:
    while ( norm(dx,Inf) > sqrt(eps(Float64))/200 && counter < 10)
        vals = innerJacobiRec(n, x, a, b)   # Evaluate via asymptotic formula.
        dx = -vals[1]./vals[2]                   # Newton update.
        x += dx                             # Next iterate.
        counter += 1
    end
    vals = innerJacobiRec(n, x, a, b);     # Evaluate via asymptotic formula.
    dx = -vals[1]./vals[2]                        # Newton update
    x += dx

    # flip:
    x = x[npts:-1:1]
    ders = vals[2]; vals = vals[1]
    ders = ders[npts:-1:1]

    # Revert to x-space:
    w = @. 1 /((1-x^2)*ders^2)
    return x, w
end

function JacobiGW( n::Integer, a::Float64, b::Float64 )
    # Golub-Welsh for Gauss--Jacobi quadrature. This is used when max(a,b)>5.
    ab = a + b;
    ii = 2:n-1;
    abi = 2*ii .+ ab;
    aa = Float64[(b - a)/(2 + ab);
          (b^2 - a^2)./((abi .- 2).*abi);
          (b^2 - a^2)./((2*n - 2+ab).*(2*n+ab))] ::Vector{Float64}
    bb = Float64[2*sqrt( (1 + a)*(1 + b)/(ab + 3))/(ab + 2) ;
          2 .*sqrt.(ii.*(ii .+ a).*(ii .+ b).*(ii .+ ab)./(abi.^2 .- 1))./abi] ::Vector{Float64}
    TT = SymTridiagonal(aa, bb)  # Jacobi matrix.
    x, V = eigen( TT )                       # Eigenvalue decomposition.
    # Quadrature weights:
    w = V[1,:].^2 .*( 2^(ab+1)*gamma(a+1)*gamma(b+1)/gamma(2+ab) );
    x, vec(w)
end


text:
"""
    struct Polygon{N,T} <: AbstractPlanarSurfacePrimitive{T}

Surface primitive describing a planar polygon, e.g. the base of a [`RegularPrism`](@ref).

## Parametric types
* `T`: Precision type.
* `N`: Number of vertices of the `Polygon`.

## Fields
* `points::SVector{N, CartesianPoint{T}}`: Vertices of the polygon in the order with which they are connected.
"""
struct Polygon{N,T} <: AbstractPlanarSurfacePrimitive{T}
    points::SVector{N, CartesianPoint{T}}
end
const Triangle{T} = Polygon{3, T}
const Quadrangle{T} = Polygon{4, T}

normal(p::Polygon) = normalize((p.points[2] - p.points[1]) × (p.points[3] - p.points[1]))

vertices(p::Polygon) = p.points
vertices(p::Polygon, n::Int64) = vertices(p)

function _sample_excluding_border(t::Triangle{T}, spacing::T)::Vector{CartesianPoint{T}} where {T}
    push = spacing/2.5#best value for not showing triangle edges
    u = t.points[2] - t.points[1]
    v = t.points[3] - t.points[1]
    w = t.points[3] - t.points[2]
    nu = norm(u)
    nv = norm(v)
    nw = norm(w)
    sθ = norm(u × v)/(nu*nv)
    sφ = norm(u × w)/(nu*nw)
    p1 = t.points[1] + (push/sθ)*(u/nu + v/nv)
    p2 = t.points[2] + (push/sφ)*(-u/nu + w/nw)
    if (p2-p1) ⋅ u > 0
        su = norm(p2-p1)/nu
        [(su*a*u + su*b*v) + p1 for a in range(0, 1, length = max(2, 1 + Int(ceil(su*nu*sθ/spacing)))) for b in range(0, (1 - a), length = max(2, 1 + Int(ceil(su*nv*(1 - a)/spacing))))]
    else
        []
    end
end

triangles(p::Polygon{N,T}) where {N,T} = [Triangle{T}([p.points[1], p.points[i], p.points[i+1]]) for i in 2:N-1]
triangles(p::Polygon, n_arc::Int64) = triangles(p)

function sample(p::Polygon{N,T}, spacing::T)::Vector{CartesianPoint{T}} where {N,T}
    v = [s for e in edges(p) for s in sample(e,n=max(2,Int(ceil(norm(e.b-e.a)/spacing))))]
    for t in triangles(p)
        append!(v, _sample_excluding_border(t, spacing))
    end
    v
end

function extremum(p::Polygon{N,T})::T where {N,T}
    c = sum(p.points)/N
    m = maximum([norm(point - c) for point in p.points])
end

connections(p::Polygon{N}) where {N} = [collect(1:N)]
connections(p::Polygon, ::Int64) = connections(p)
connections(p::Polygon{N}, ::Int64, ::Int64) where {N} = [[i%N+1, (i+1)%N+1] for i in 0:N-1]
get_label_name(::Polygon) = "Polygon"

extreme_points(p::Polygon) = p.points

function edges(p::Triangle{T}) where {T}
    vs = vertices(p)
    return (
        Edge(vs[1], vs[2]),
        Edge(vs[2], vs[3]),
        Edge(vs[3], vs[1])
    )
end
function edges(p::Quadrangle{T}) where {T}
    vs = vertices(p)
    return (
        Edge(vs[1], vs[2]),
        Edge(vs[2], vs[3]),
        Edge(vs[3], vs[4]),
        Edge(vs[4], vs[1])
    )
end

function edges(p::Polygon{N,T})::SVector{N, Edge{T}} where {N,T}
    vs = vertices(p)
    SVector{N, Edge{T}}([Edge(vs[i],vs[i%N+1]) for i in 1:N])
end

flip(p::Polygon) = Polygon(reverse(p.points))

lines(p::Polygon) = edges(p)

Plane(p::Polygon{N, T}) where {N, T} = Plane{T}(p.points[1], (p.points[2] - p.points[1]) × (p.points[3] - p.points[1]))

function _get_rot_for_rotation_on_xy_plane(p::Polygon{<:Any, T}) where {T}
    n = normal(p)
    rot_angle_rad = CartesianVector{T}(0,0,1) ⋅ n
    return if abs(rot_angle_rad) == 1
        RotationVec(zero(T), zero(T), zero(T))
    else
        rot_angle = -acos(rot_angle_rad)
        rot_axis = rot_angle * normalize(CartesianVector{T}(0,0,1) × n)
        RotationVec(rot_axis[1], rot_axis[2], rot_axis[3])
    end
end

function _rotate_on_xy_plane(p::Polygon{<:Any, T}) where {T}
    rot = _get_rot_for_rotation_on_xy_plane(p)
    map(pt -> rot * pt, p.points)
end

function in(pt::CartesianPoint{T}, p::Polygon{N,T}) where {N,T}
    plane = Plane(p)
    if (pt - origin(plane)) ⋅ normal(plane) == 0
        rot = _get_rot_for_rotation_on_xy_plane(p)
        vs = vertices(p)
        pts2d = SVector{N+1, SVector{2, T}}(
            view(rot * vs[i%N + 1], 1:2) for i in 0:N
        )
        # PolygonOps.inpolygon -> in = 1, on = -1, out = 0)
        b = PolygonOps.inpolygon(view(rot * pt, 1:2), pts2d) != 0
    end
    return b
end

function _above_or_below_polygon(pt::AbstractCoordinatePoint, p::Quadrangle{T}) where {T}
    rot = _get_rot_for_rotation_on_xy_plane(p)
    vs = vertices(p)
    pts2d = SVector{5, SVector{2, T}}(
        view(rot * vs[1], 1:2),
        view(rot * vs[2], 1:2),
        view(rot * vs[3], 1:2),
        view(rot * vs[4], 1:2),
        view(rot * vs[1], 1:2),
    )
    # PolygonOps.inpolygon -> in = 1, on = -1, out = 0)
    return PolygonOps.inpolygon(view(rot * pt, 1:2), pts2d) != 0
end

function _transform_into_global_coordinate_system(poly::Polygon{N, T}, p::AbstractPrimitive{T}) where {N, T}
    Polygon{N,T}(broadcast(pt -> _transform_into_global_coordinate_system(pt, p), poly.points))
end

function distance(pt::CartesianPoint, p::Polygon)
    return if _above_or_below_polygon(pt, p)
        distance(pt, Plane(p))
    else
        es = edges(p)
        ds = map(e -> distance(pt, e), es)
        minimum(ds)
    end
end


text:
module TravelingSalesmanExact

using JuMP, UnicodePlots, Logging, LinearAlgebra
import MathOptInterface
const MOI = MathOptInterface
export get_optimal_tour,
       plot_cities,
       simple_parse_tsp,
       with_optimizer,
       set_default_optimizer!

# added in Julia 1.2
if VERSION < v"1.2.0-"
    import Base.>
    >(x) = Base.Fix2(>, x)
end

# added in Julia 1.1
if VERSION < v"1.1.0-"
    isnothing(::Any) = false
    isnothing(::Nothing) = true
end

const default_optimizer = Ref{Union{OptimizerFactory,Nothing}}(nothing)

"""
    set_default_optimizer(O::OptimizerFactory)

Sets the default optimizer. For example,

    using GLPK
    set_default_optimizer(with_optimizer(GLPK.Optimizer))
"""
set_default_optimizer!(O::OptimizerFactory) = default_optimizer[] = O

"""
    get_default_optimizer()

Gets the default optimizer, which is set by `set_default_optimizer`.
"""
get_default_optimizer() = default_optimizer[]

reset_default_optimizer!() = default_optimizer[] = nothing

"""
    plot_cities(cities)

Uses `UnicodePlots`'s `lineplot` to make a plot of the tour of the cities in
`cities`, in order (including going from the last city back to the first).
"""
function plot_cities(cities)
    n = length(cities)
    inc(a) = a == n ? one(a) : a + 1
    lineplot([cities[inc(j)][1] for j = 0:n], [cities[inc(j)][2] for j = 0:n])
end

"""
    find_cycle(perm_matrix, starting_ind)

Returns the cycle in the permutation described by `perm_matrix` which includes
`starting_ind`.
"""
function find_cycle(perm_matrix, starting_ind = 1)
    cycle = [starting_ind]
    prev_ind = ind = starting_ind
    while true
        # the comparisons `x > (0.5)` should mean `x == 1`. Due to floating point results returned
        # by the solvers, instead we sometimes have `x ≈ 1.0` instead. Since these are binary
        # values, we might as well just compare to 1/2.
        next_ind = findfirst(>(0.5), @views(perm_matrix[ind, 1:prev_ind-1]))
        if isnothing(next_ind)
            next_ind = findfirst(>(0.5), @views(perm_matrix[ind, prev_ind+1:end])) +
                       prev_ind
        end
        next_ind == starting_ind && break
        push!(cycle, next_ind)
        prev_ind, ind = ind, next_ind
    end
    cycle
end

"""
    get_cycles(perm_matrix)

Returns a list of cycles from the permutation described by `perm_matrix`.
"""
function get_cycles(perm_matrix)
    N = size(perm_matrix, 1)
    remaining_inds = Set(1:N)
    cycles = Vector{Int}[]
    while length(remaining_inds) > 0
        cycle = find_cycle(perm_matrix, first(remaining_inds))
        push!(cycles, cycle)
        setdiff!(remaining_inds, cycle)
    end
    cycles
end

"""
    show_tour(cities, perm_matrix)

Show a plot of the tour described by `perm_matrix` of the cities in the vector `cities`.
"""
function plot_tour(cities, perm_matrix)
    cycles = get_cycles(perm_matrix)
    tour = reduce(vcat, cycles)
    plot_cities(cities[tour])
end

"""
    remove_cycles!(model, tour_matrix)

Find the (non-maximal-length) cycles in the current solution `tour_matrix`
and add constraints to the JuMP model to disallow them. Returns the
number of cycles found.
"""
function remove_cycles!(model, tour_matrix; symmetric)
    tour_matrix_val = value.(tour_matrix)
    cycles = get_cycles(tour_matrix_val)
    length(cycles) == 1 && return 1
    for cycle in cycles
        constr = symmetric ? 2 * length(cycle) - 2 : length(cycle) - 1
        @constraint(model, sum(tour_matrix[cycle, cycle]) <= constr)
    end
    return length(cycles)
end

"""
    euclidean_distance(city1, city2)

The usual Euclidean distance measure.
"""
euclidean_distance(city1, city2) = sqrt((city1[1] - city2[1])^2 + (city1[2] - city2[2])^2)

"""
    ATT(city1, city2)

The `ATT` distance measure as specified in TSPLIB:
<https://www.iwr.uni-heidelberg.de/groups/comopt/software/TSPLIB95/tsp95.pdf>.
"""
function ATT(city1, city2)
    xd = city1[1] - city2[1]
    yd = city1[2] - city2[2]
    r = sqrt((xd^2 + yd^2) / 10.0)
    t = round(Int, r)
    if t < r
        d = t + 1
    else
        d = t
    end
    return d
end

"""
    get_optimal_tour(
        cities::AbstractVector,
        with_optimizer = get_default_optimizer();
        verbose = false,
        distance = euclidean_distance,
        symmetric = true,
    )

Solves the travelling salesman problem for a list of cities using JuMP by
formulating a MILP using the Dantzig-Fulkerson-Johnson formulation and
adaptively adding constraints to disallow non-maximal cycles. Returns an optimal
tour and the cost of the optimal path. Optionally specify a distance metric.

The second argument is mandatory if a default optimizer has not been set (via
`set_default_optimizer`). This argument should be the result of a call to
`JuMP.with_optimizer`, e.g.

    get_optimal_tour(cities, with_optimizer(GLPK.Optimizer))
"""
function get_optimal_tour(
    cities::AbstractVector,
    with_optimizer = get_default_optimizer();
    verbose = false,
    distance = euclidean_distance,
    symmetric = true,
)
    isnothing(with_optimizer) && throw(ArgumentError("An optimizer is required if a default optimizer has not been set."))
    N = length(cities)
    cost = [distance(cities[i], cities[j]) for i = 1:N, j = 1:N]
    return _get_optimal_tour(cost, with_optimizer, symmetric, verbose, cities)
end

"""
    get_optimal_tour(
        cost::AbstractMatrix,
        with_optimizer = get_default_optimizer();
        verbose = false,
        symmetric = issymmetric(cost),
    )

Solves the travelling salesman problem for a square cost matrix using JuMP by
formulating a MILP using the Dantzig-Fulkerson-Johnson formulation and
adaptively adding constraints to disallow non-maximal cycles. Returns an optimal
tour and the cost of the optimal path.

The second argument is mandatory if a default optimizer has not been set (via
`set_default_optimizer`). This argument should be the result of a call to
`JuMP.with_optimizer`, e.g.

    get_optimal_tour(cities, with_optimizer(GLPK.Optimizer))
"""
function get_optimal_tour(
    cost::AbstractMatrix,
    with_optimizer = get_default_optimizer();
    verbose = false,
    symmetric = issymmetric(cost),
)
    size(cost, 1) == size(cost, 2) || throw(ArgumentError("First argument must be a square matrix"))
    isnothing(with_optimizer) && throw(ArgumentError("An optimizer is required if a default optimizer has not been set."))
    return _get_optimal_tour(cost, with_optimizer, symmetric, verbose)
end

function _get_optimal_tour(
    cost::AbstractMatrix,
    with_optimizer,
    symmetric,
    verbose,
    cities = nothing,
)
    N = size(cost, 1)
    has_cities = !isnothing(cities)

    model = Model(with_optimizer)
    if symmetric
         # `tour_matrix` has tour_matrix[i,j] = 1 iff cities i and j should be connected
        @variable(model, tour_matrix[1:N, 1:N], Symmetric, binary = true)

        # cost of the tour
        @objective(model, Min, sum(tour_matrix[i, j] * cost[i, j] for i = 1:N, j = 1:i))
        for i = 1:N
            @constraint(model, sum(tour_matrix[i, :]) == 2) # degree of each city is 2
            @constraint(model, tour_matrix[i, i] == 0) # rule out cycles of length 1
        end
    else
        # `tour_matrix` will be a permutation matrix
        @variable(model, tour_matrix[1:N, 1:N], binary = true)
        @objective(model, Min, sum(tour_matrix[i, j] * cost[i, j] for i = 1:N, j = 1:N))
        for i = 1:N
            @constraint(model, sum(tour_matrix[i, :]) == 1) # row-sum is 1
            @constraint(model, sum(tour_matrix[:, i]) == 1) # col-sum is 1
            @constraint(model, tour_matrix[i, i] == 0) # rule out cycles of length 1
            for j = 1:N
                @constraint(model, tour_matrix[i, j] + tour_matrix[j, i] <= 1) # rule out cycles of length 2
            end
        end
    end

    if has_cities && verbose
        @info "Starting optimization." plot_cities(cities)
    elseif verbose
        @info "Starting optimization."
    end

    iter = 0 # count for logging
    tot_cycles = 0 # count for logging
    num_cycles = 2 # just something > 1
    while num_cycles > 1
        t = @elapsed optimize!(model)
        status = termination_status(model)
        status == MOI.OPTIMAL || throw(ErrorException("Error: problem status $status"))
        num_cycles = remove_cycles!(model, tour_matrix; symmetric = symmetric)
        tot_cycles += num_cycles
        iter += 1
        if has_cities && verbose
            @info "Iteration $iter took $(round(t, digits=3))s, disallowed $num_cycles cycles." plot_tour(
                cities,
                value.(tour_matrix),
            )
        elseif verbose
            @info "Iteration $iter took $(round(t, digits=3))s, disallowed $num_cycles cycles."
        end
    end
    tot_cycles -= 1 # remove the true cycle

    status = termination_status(model)
    status == MOI.OPTIMAL || @warn(status)

    if verbose
        @info "Optimization finished; adaptively disallowed $tot_cycles cycles."
        @info "Final path has length $(objective_value(model))."
        @info "Final problem has $(num_constraints(model, VariableRef, MOI.ZeroOne)) binary variables,
            $(num_constraints(model,
            GenericAffExpr{Float64,VariableRef}, MOI.LessThan{Float64})) inequality constraints, and
            $(num_constraints(model, GenericAffExpr{Float64,VariableRef}, MOI.EqualTo{Float64})) equality constraints."
    end
    return find_cycle(value.(tour_matrix)), objective_value(model)
end

"""
    simple_parse_tsp(filename; verbose = true)

Try to parse the ".tsp" file given by `filename`. Very simple implementation
just to be able to test the optimization; may break on other files. Returns a
list of cities for use in `get_optimal_tour`.
"""
function simple_parse_tsp(filename; verbose = true)
    cities = Vector{Int}[]
    for line in readlines(filename)
        if startswith(line, '1':'9')
            nums = split(line, " ")
            @assert length(nums) == 3
            x = parse(Int, nums[2])
            y = parse(Int, nums[3])
            push!(cities, [x, y])
        elseif verbose
            println(line)
        end
    end
    return cities
end

end # module


text:
using Test

using ChipSort
using SIMD

include("testutils.jl")

@testset "Bitonic Merge" begin

function test_merge_types_sizes(T, N)
    va = Vec(tuple(sort(rand(T, N))...))
    vb = Vec(tuple(sort(rand(T, N))...))
    a, b = bitonic_merge(va, vb)
    ab = [[a[n] for n in 1:N]; [b[n] for n in 1:N]]
    abref = sort([[va[n] for n in 1:N]; [vb[n] for n in 1:N]])
    @test ab == abref
    @test all(ab[2:end] .>= ab[1:end-1])
end

@testset "Basic bitonic merge" begin
    for T in [Int8, Int16, Int32, Int64, Float32, Float64]
        for N in 2 .^ (1:3)
            test_merge_types_sizes(T, N)
        end
    end
end

function test_merge_heaviside(N, Na, Nb)
    va = Vec(ntuple(n->if n<=Na 0x0 else 0x1 end, N))
    vb = Vec(ntuple(n->if n<=Nb 0x0 else 0x1 end, N))
    a, b = bitonic_merge(va, vb)
    ab = [[a[n] for n in 1:N]; [b[n] for n in 1:N]]
    @test all(ab[1:Na+Nb] .== 0)
    @test all(ab[Na+Nb+1:end] .== 1)
end

@testset "Exaustive binary arrays" begin
    for N in 2 .^ (1:8)
        for Na in 0:N
            for Nb in 0:N
                test_merge_heaviside(N, Na, Nb)
            end
        end
    end
end

function test_merge_vecs(T, N, L)
    T=Int64
    N=8
    aa = sort(rand(T, N, L); dims=1)
    va = ntuple(k->Vec(ntuple(j->aa[j,k],N)), L)
    srtref = sort(aa[:])

    srthat = merge_vecs(va...)
    srthat_arr = [srthat[n] for n in 1:N*L]
    @test all(srthat_arr .== srtref)
end

@testset "Multiple vectors" begin
    for T in [Int8, Int16, Int32, Int64, Float32, Float64]
        for N in [4, 16]
            for L in [4, 16]
                test_merge_vecs(T, N, L)
            end
        end
    end
end

function test_merge_interleaved(T, V, J, K)
    aa = randa(T, V*J*K)
    for k in 1:K
        sort!(@view aa[1+(k-1)*J*V:k*J*V])
    end

    srt = chipsort_merge_medium(aa, Val(V), Val(J), Val(K))
    @test srt == sort(aa)
end

@testset "Basic bitonic merge interleaved" begin
    for T in [Int8, Int32, Float64]
        for V in 2 .^ (2:3)
            for J in [2,8]
                for K in [2,4]
                    test_merge_interleaved(T,V,J,K)
                end
            end
        end
    end
end

end


text:
function sparse_from_coo(::Type{<:SparseMatrixCSC},I,J,V,m,n)
  sparse(I,J,V,m,n)
end

@inline function is_entry_stored(::Type{<:SparseMatrixCSC},i,j)
  true
end

function finalize_coo!(::Type{<:SparseMatrixCSC},I,J,V,m,n)
  nothing
end

function nz_index(A::SparseArrays.AbstractSparseMatrixCSC,i0,i1)
    if !(1 <= i0 <= size(A, 1) && 1 <= i1 <= size(A, 2)); throw(BoundsError()); end
    ptrs = SparseArrays.getcolptr(A)
    r1 = Int(ptrs[i1])
    r2 = Int(ptrs[i1+1]-1)
    (r1 > r2) && return -1
    r1 = searchsortedfirst(rowvals(A), i0, r1, r2, Base.Order.Forward)
    ((r1 > r2) || (rowvals(A)[r1] != i0)) ? -1 : r1
end

@inline function push_coo!(::Type{<:SparseMatrixCSC},I,J,V,i,j,v)
 push!(I,i)
 push!(J,j)
 push!(V,v)
 nothing
end

#@inline function add_entries!(
#  combine::Function,
#  A::SparseMatrixCSC,
#  vs::AbstractMatrix{<:Number},
#  is,js)
#
#  if issorted(is)
#    nz = A.nzval
#    ptrs = A.colptr
#    rows = A.rowval
#    for (lj,j) in enumerate(js)
#      if j>0
#        pini = ptrs[j]
#        pend = ptrs[j+1]-1
#        li = 1
#        for p in pini:pend
#          _i = rows[p]
#          if _i == is[li]
#            vij = vs[li,lj]
#            Aij = nz[p]
#            nz[p] = combine(Aij,vij)
#            li += 1
#          end
#        end
#      end
#    end
#  else
#    for (lj,j) in enumerate(js)
#      if j>0
#        for (li,i) in enumerate(is)
#          if i>0
#            vij = vs[li,lj]
#            add_entry!(combine,A,vij,i,j)
#          end
#        end
#      end
#    end
#  end
#  A
#end

struct CounterCSC{Tv,Ti,L}
  tv::Type{Tv}
  nrows::Int
  ncols::Int
  colnnzmax::Vector{Ti}
  loop_style::L
end

LoopStyle(::Type{CounterCSC{Tv,Ti,L}}) where {Tv,Ti,L} = L()

@inline function add_entry!(::typeof(+),a::CounterCSC{Tv,Ti,Loop},v,i,j) where {Tv,Ti}
  a.colnnzmax[j] += Ti(1)
  nothing
end

@inline function add_entry!(::typeof(+),a::CounterCSC{Tv,Ti,DoNotLoop},v,i,j) where {Tv,Ti}
  nothing
end

struct InserterCSC{Tv,Ti}
  nrows::Int
  ncols::Int
  colptr::Vector{Ti}
  colnnz::Vector{Ti}
  rowval::Vector{Ti}
  nzval::Vector{Tv}
end

LoopStyle(::Type{<:InserterCSC}) = Loop()

@inline function add_entry!(::typeof(+),a::InserterCSC{Tv,Ti},v::Nothing,i,j) where {Tv,Ti}
  pini = Int(a.colptr[j])
  pend = pini + Int(a.colnnz[j]) - 1
  p = searchsortedfirst(a.rowval,i,pini,pend,Base.Order.Forward)
  if (p>pend)
    # add new entry
    a.colnnz[j] += 1
    a.rowval[p] = i
  elseif a.rowval[p] != i
    # shift one forward from p to pend
    @check  pend+1 < Int(a.colptr[j+1])
    for k in pend:-1:p
      o = k + 1
      a.rowval[o] = a.rowval[k]
    end
    # add new entry
    a.colnnz[j] += 1
    a.rowval[p] = i
  end
  nothing
end

@noinline function add_entry!(::typeof(+),a::InserterCSC{Tv,Ti},v,i,j) where {Tv,Ti}
  pini = Int(a.colptr[j])
  pend = pini + Int(a.colnnz[j]) - 1
  p = searchsortedfirst(a.rowval,i,pini,pend,Base.Order.Forward)
  if (p>pend)
    # add new entry
    a.colnnz[j] += 1
    a.rowval[p] = i
    a.nzval[p] = v
  elseif a.rowval[p] != i
    # shift one forward from p to pend
    @check  pend+1 < Int(a.colptr[j+1])
    for k in pend:-1:p
      o = k + 1
      a.rowval[o] = a.rowval[k]
      a.nzval[o] = a.nzval[k]
    end
    # add new entry
    a.colnnz[j] += 1
    a.rowval[p] = i
    a.nzval[p] = v
  else
    # update existing entry
    a.nzval[p] += v
  end
  nothing
end

## index of the first value of vector a that is greater than or equal to x;
## returns length(v)+1 if x is greater than all values in v.
#function _searchsortedfirst(v, x, lo::T, hi::T) where T<:Integer
#  if x <= v[lo]
#    return lo
#  end
#  u = one(T)
#  if x > v[hi]
#    return hi + u
#  end
#  d = T(10) - u
#  @inbounds while lo < hi - d
#    m = Base.Sort.midpoint(lo, hi)
#    if v[m] < x
#      lo = m
#    else
#      hi = m
#    end
#  end
#  i = lo
#  if x<=v[i]
#    return i
#  end
#  i += u
#  if x<=v[i]
#    return i
#  end
#  i += u
#  if x<=v[i]
#    return i
#  end
#  i += u
#  if x<=v[i]
#    return i
#  end
#  i += u
#  if x<=v[i]
#    return i
#  end
#  i += u
#  if x<=v[i]
#    return i
#  end
#  i += u
#  if x<=v[i]
#    return i
#  end
#  i += u
#  if x<=v[i]
#    return i
#  end
#  i += u
#  if x<=v[i]
#    return i
#  end
#  i += u
#  if x<=v[i]
#    return i
#  end
#  return hi + u
#end

function nz_counter(
  builder::SparseMatrixBuilder{SparseMatrixCSC{Tv,Ti},<:MinMemory},
  axes) where {Tv,Ti}
  nrows = length(axes[1])
  ncols = length(axes[2])
  maxnnz = builder.approach.maxnnz
  if isa(maxnnz,Nothing)
    colnnzmax = zeros(Ti,ncols)
    CounterCSC(Tv,nrows,ncols,colnnzmax,Loop())
  elseif isa(maxnnz,Integer)
    colnnzmax = fill(Ti(maxnnz),ncols)
    CounterCSC(Tv,nrows,ncols,colnnzmax,DoNotLoop())
  elseif isa(maxnnz,Vector{<:Integer})
    colnnzmax = maxnnz
    @assert length(colnnzmax) == ncols
    CounterCSC(Tv,nrows,ncols,colnnzmax,DoNotLoop())
  else
    @notimplemented
  end
end

function nz_counter(
  builder::SparseMatrixBuilder{SparseMatrixCSC{Tv,Ti},<:MinCPU},
  axes) where {Tv,Ti}

  nrows = length(axes[1])
  ncols = length(axes[2])
  rowptrs = zeros(Ti,nrows+1)
  CounterCSRR(Tv,nrows,ncols,rowptrs)
end

#function nz_counter(::Type{SparseMatrixCSC{Tv,Ti}},axes) where {Tv,Ti}
#  builder = SparseMatrixBuilder(SparseMatrixCSC{Tv,Ti},MinMemory())
#  nz_counter(builder,axes)
#end

function nz_allocation(a::CounterCSC{Tv,Ti}) where {Tv,Ti}
  colptr = Vector{Ti}(undef,a.ncols+1)
  @inbounds for i in 1:a.ncols
    colptr[i+1] = a.colnnzmax[i]
  end
  length_to_ptrs!(colptr)
  ndata = colptr[end] - one(Ti)
  rowval = Vector{Ti}(undef,ndata)
  nzval = zeros(Tv,ndata)
  colnnz = a.colnnzmax
  fill!(colnnz,zero(Ti))
  InserterCSC(a.nrows,a.ncols,colptr,colnnz,rowval,nzval)
end

function create_from_nz(a::InserterCSC{Tv,Ti}) where {Tv,Ti}
  k = 1
  for j in 1:a.ncols
    pini = Int(a.colptr[j])
    pend = pini + Int(a.colnnz[j]) - 1
    for p in pini:pend
      a.nzval[k] = a.nzval[p]
      a.rowval[k] = a.rowval[p]
      k += 1
    end
  end
  @inbounds for j in 1:a.ncols
    a.colptr[j+1] = a.colnnz[j]
  end
  length_to_ptrs!(a.colptr)
  nnz = a.colptr[end]-1
  resize!(a.rowval,nnz)
  resize!(a.nzval,nnz)
  SparseMatrixCSC(a.nrows,a.ncols,a.colptr,a.rowval,a.nzval)
end

struct CounterCSRR{Tv,Ti}
  tv::Type{Tv}
  nrows::Int
  ncols::Int
  rowptrs::Vector{Ti}
end

LoopStyle(::Type{<:CounterCSRR}) = Loop()

@inline function add_entry!(::typeof(+),a::CounterCSRR{Tv,Ti},v,i,j) where {Tv,Ti}
  a.rowptrs[i+1] += Ti(1)
  nothing
end

struct CSRR{Tv,Ti}
  nrows::Int
  ncols::Int
  rowptrs::Vector{Ti}
  colvals::Vector{Ti}
  nzvals::Vector{Tv}
  work::Vector{Ti}
end

LoopStyle(::Type{<:CSRR}) = Loop()

@inline function add_entry!(::typeof(+),a::CSRR{Tv,Ti},v::Nothing,i,j) where {Tv,Ti}
  p = a.rowptrs[i]
  a.colvals[p] = j
  a.rowptrs[i] = p+Ti(1)
  nothing
end

@inline function add_entry!(::typeof(+),a::CSRR{Tv,Ti},v,i,j) where {Tv,Ti}
  p = a.rowptrs[i]
  a.colvals[p] = j
  a.nzvals[p] = v
  a.rowptrs[i] = p+Ti(1)
  nothing
end

function nz_allocation(a::CounterCSRR{Tv,Ti}) where {Tv,Ti}
  rowptrs = a.rowptrs
  length_to_ptrs!(rowptrs)
  ndata = rowptrs[end]-1
  colvals = Vector{Ti}(undef,ndata)
  nzvals = zeros(Tv,ndata)
  work = Vector{Ti}(undef,a.ncols)
  CSRR(a.nrows,a.ncols,rowptrs,colvals,nzvals,work)
end

function create_from_nz(a::CSRR{Tv,Ti}) where {Tv,Ti}
  rewind_ptrs!(a.rowptrs)
  colptrs = Vector{Ti}(undef,a.ncols+1)
  cscnnz = _csrr_to_csc_count!(colptrs,a.rowptrs,a.colvals,a.nzvals,a.work)
  rowvals = Vector{Ti}(undef,cscnnz)
  nzvalscsc = Vector{Tv}(undef,cscnnz)
  _csrr_to_csc_fill!(colptrs,rowvals,nzvalscsc,a.rowptrs,a.colvals,a.nzvals)
  SparseMatrixCSC(a.nrows,a.ncols,colptrs,rowvals,nzvalscsc)
end

# Notation
# csrr: csr with repeated and unsorted columns
# csru: csr witu unsorted columns
# csc: csc with sorted columns

# Adapted form SparseArrays
function _csrr_to_csc_count!(
  colptrs::Vector{Ti},
  rowptrs::Vector{Tj},
  colvals::Vector{Tj},
  nzvalscsr::Vector{Tv},
  work::Vector{Tj}) where {Ti,Tj,Tv}

  nrows = length(rowptrs)-1
  ncols = length(colptrs)-1
  if nrows == 0 || ncols == 0
    fill!(colptrs, Ti(1))
    return Tj(0)
  end

  # Convert csrr to csru by identifying repeated cols with array work.
  # At the same time, count number of unique rows in colptrs shifted by one.
  fill!(colptrs, Ti(0))
  fill!(work, Tj(0))
  writek = Tj(1)
  newcsrrowptri = Ti(1)
  origcsrrowptri = Tj(1)
  origcsrrowptrip1 = rowptrs[2]
  @inbounds for i in 1:nrows
    for readk in origcsrrowptri:(origcsrrowptrip1-Tj(1))
      j = colvals[readk]
      if work[j] < newcsrrowptri
        work[j] = writek
        if writek != readk
          colvals[writek] = j
          nzvalscsr[writek] = nzvalscsr[readk]
        end
        writek += Tj(1)
        colptrs[j+1] += Ti(1)
      else
        klt = work[j]
        nzvalscsr[klt] = +(nzvalscsr[klt], nzvalscsr[readk])
      end
    end
    newcsrrowptri = writek
    origcsrrowptri = origcsrrowptrip1
    origcsrrowptrip1 != writek && (rowptrs[i+1] = writek)
    i < nrows && (origcsrrowptrip1 = rowptrs[i+2])
  end

  # Convert colptrs from counts to ptrs shifted by one
  # (ptrs will be corrected below)
  countsum = Tj(1)
  colptrs[1] = Ti(1)
  @inbounds for j in 2:(ncols+1)
    overwritten = colptrs[j]
    colptrs[j] = countsum
    countsum += overwritten
    @check Base.hastypemax(Ti) && (countsum <= typemax(Ti))
  end

  cscnnz = countsum - Tj(1)
  cscnnz
end

function _csrr_to_csc_fill!(
  colptrs::Vector{Ti},rowvals::Vector{Ti},nzvalscsc::Vector{Tv},
  rowptrs::Vector{Tj},colvals::Vector{Tj},nzvalscsr::Vector{Tv}) where {Ti,Tj,Tv}

  nrows = length(rowptrs)-1
  ncols = length(colptrs)-1
  if nrows == 0 || ncols == 0
    return nothing
  end

  # From csru to csc
  # Tracking write positions in colptrs corrects
  # the column pointers to the final value.
  @inbounds for i in 1:nrows
    for csrk in rowptrs[i]:(rowptrs[i+1]-Tj(1))
      j = colvals[csrk]
      x = nzvalscsr[csrk]
      csck = colptrs[j+1]
      colptrs[j+1] = csck + Ti(1)
      rowvals[csck] = i
      nzvalscsc[csck] = x
    end
  end

  nothing
end


text:
using SafeTestsets

@safetestset "Beam" begin include("BeamTestset.jl") end
@safetestset "Column" begin include("ColumnTestset.jl") end
@safetestset "InternalForces" begin include("InternalForcesTestset.jl") end
@safetestset "AISIS10016" begin include("AISIS10016Testset.jl") end
# @safetestset "AISIS10024" begin include("AISIS10024Testset.jl") end


text:
using Plots, Statistics, MLDatasets, LinearAlgebra, Clustering, Printf, Dates, Random, Debugger
path1=realpath(dirname(@__FILE__)*"/../..")
include(string(path1,"/src/ANNBN.jl"))

i_train=100
i_test=i_train-1
vars=1
xx_train=range(0.05,0.95,length=i_train)
yy_train=[zeros(50);ones(50)]
tol1=0.0000000001
ANNBN.isigm1(tol1)
ANNBN.isigm1(1.0-tol1)
mi=minimum(yy_train);yy_train.-=mi;ma=maximum(yy_train);yy_train./=ma;yy_train.*=1.0-2tol1;yy_train.+=tol1
scatter(xx_train,yy_train,label="input yy_train")
xx_test=(xx_train[2:end].+xx_train[1:end-1])./2
yy_test=(yy_train[2:end].+yy_train[1:end-1])./2

# neurons=Int64(floor(i_train/(vars+1)))
# neurons=4
# this automatically generates number of neurons which equals to the number of clusters (§3.1)
# inds_all,n_per_part,items_per_neuron=ANNBN.___clustering(neurons,xx_train,200)

inds_all=(1:i_train)[:,1]
# inds_all=randperm(i_train)
div1=20
neurons=Int64(i_train/div1)
items_per_neuron=div1*ones(Int64,neurons)   # (Int64(floor(i_train/(neurons))))
n_per_part=[0;cumsum(items_per_neuron)]

@time a_all,a_layer1,layer1=ANNBN.train_layer_1_sigmoid_fast(neurons,vars,i_train,n_per_part,inds_all,xx_train,yy_train)
# i1=reduce(hcat, a_all).==[0;0]
# i2=i1[1,:].+i1[2,:]
# sum(i2.==2)
#
# @time ll1=layer1'*layer1
# @time det(ll1)
# @time ll1i=inv(ll1)
# @time a_layer1=ll1i*layer1'*y
# plot(a_all)
# det(layer1'*layer1)
predl1,layer1_train=ANNBN.predict_new(a_all,a_layer1,xx_train,i_train,neurons)
predl1_test,layer1_test=ANNBN.predict_new(a_all,a_layer1,xx_test,i_test,neurons)
maetr=mean(abs.(yy_train-predl1))
maete=mean(abs.(yy_test-predl1_test))
scatter!(xx_train,predl1,label=@sprintf("ANNBN-%1d-neurons",neurons))
scatter!(xx_test,predl1_test,label=@sprintf("ANNBN-%1d-neurons",neurons))
i1=predl1.<=0.5;predl1[i1].=tol1;i1=predl1.>0.5;predl1[i1].=1.0-tol1;i2=abs.(predl1.-yy_train).>0.01;maetr=sum(i2)/length(i2)
i1=predl1_test.<=0.5;predl1_test[i1].=tol1;i1=predl1_test.>0.5;predl1_test[i1].=1.0-tol1;i2=abs.(predl1_test.-yy_test).<0.01;
maete=1.0-sum(i2)/length(i2)

# plot(ANNBN.isigm1.(yy_train))

@time bb=[ones(i_train) xx_train]\ANNBN.isigm1.(yy_train)
# # pinv([ones(i_train) xx_train])*[ones(i_train) xx_train]
# X=[ones(i_train) xx_train]
# Xt=[ones(i_train) xx_train]'
# XtX=Xt*X
# @time XtX_i=inv(XtX)
# det(XtX)
# @time bb=XtX_i*Xt*ANNBN.isigm1.(yy_train)
# # bb=[-500;1000]
predl1=ANNBN.sigm1.([ones(i_train) xx_train]*bb);
maetr=mean(abs.(yy_train-predl1))
# i1=predl1.<=0.5;predl1[i1].=0.0001;i1=predl1.>0.5;predl1[i1].=0.9998;i2=abs.(predl1.-yy_train).<0.01;maetr=sum(i2)/length(i2)
predl1_test=ANNBN.sigm1.([ones(i_test) xx_test]*bb);
maete=mean(abs.(yy_test-predl1_test))
# i1=predl1_test.<=0.5;predl1_test[i1].=0.0001;i1=predl1_test.>0.5;predl1_test[i1].=0.9998;i2=abs.(predl1_test.-yy_test).<0.01;maete=1.0-sum(i2)/length(i2)
scatter!(xx_train,predl1,label="Logistic",legend=false)
i1=predl1.<=0.5;predl1[i1].=tol1;i1=predl1.>0.5;predl1[i1].=1.0-tol1;i2=abs.(predl1.-yy_train).>0.01;maetr=sum(i2)/length(i2)
i1=predl1_test.<=0.5;predl1_test[i1].=tol1;i1=predl1_test.>0.5;predl1_test[i1].=1.0-tol1;i2=abs.(predl1_test.-yy_test).<0.01;
maete=1.0-sum(i2)/length(i2)

xx_train=collect(xx_train)[:,:];xx_test=collect(xx_test)[:,:]
cc1=100; @time a_all,a_layer1,layer1=ANNBN.train_layer_1_rbf(neurons,vars,i_train,n_per_part,inds_all,xx_train,yy_train,cc1);
predl1=[layer1 ones(i_train)]*a_layer1
maetr=mean(abs.(yy_train-predl1))
maximum(abs.(yy_train-predl1))
i1=predl1.<=0.5;predl1[i1].=0.0001;i1=predl1.>0.5;predl1[i1].=0.9998;i2=abs.(predl1.-yy_train).<0.01;maetr=sum(i2)/length(i2)
@time predl1_test,layer1_test=ANNBN.predict_new_rbf(a_all,a_layer1,xx_test,i_test,neurons,n_per_part,inds_all,xx_train,vars,cc1);
maete=mean(abs.(yy_test-predl1_test))
maximum(abs.(yy_test-predl1_test))
i1=predl1_test.<=0.5;predl1_test[i1].=0.0001;i1=predl1_test.>0.5;predl1_test[i1].=0.9998;i2=abs.(predl1_test.-yy_test).<0.01;maete=1.0-sum(i2)/length(i2)