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

using Markdown
using InteractiveUtils

# ╔═╡ 0cd92d78-1da6-435b-ad8e-60af29502e1d
begin
	# Necessary packages
	using PlutoUI, Random, LinearAlgebra, LightGraphs, GraphPlot
	using SparseArrays, Plots, Distances, Arpack
end

# ╔═╡ 7ef4b0b3-9d0a-4356-bd16-14487f14e0ab
TableOfContents(title="📚 Table of Contents", aside=true)

# ╔═╡ cb61f761-b2ee-4818-9149-f85d95b76a1e
md"""
# Spectral Graph Bipartitioning

Many data clustering problems can be interpreted as clustering of vertices of graphs. __Graph bipartitioning problem__ is to partition vertices into subsets such that the connections within subsets are stronger than the connections between different subsets.

Partition of the vertices into two subsetts is done according to signs of the eigenvectors of the second smallest eigenvalue of the Laplacian matrix.

__Prerequisites__

The reader should be familiar with the basic graph theory, linear algebra and, in particular,  eigenvalues and eigenvectors.

__Competences__

The reader should be able to apply graph spectral bipartitioning and recursive bipartitioning to data clustering problems.

__Credits.__ The notebook was initially derived from M.Sc. Thesis of Ivančica Mirošević.
"""

# ╔═╡ 5926c79a-6038-4dbc-b49a-5b05d02250a0
md"""
# Graphs

For more details, see [W. H. Haemers, Matrices and Graphs, in L. Hogben, Ed., 'Handbook of Linear Algebra', pp. 39.1-39.14, CRC Press, Boca Raton, 2014.](https://www.routledge.com/Handbook-of-Linear-Algebra/Hogben/p/book/9781138199897) and [S. Butler and F. Chung, Spectral Graph Theory, ibid., pp. 47.1-47.6](https://www.routledge.com/Handbook-of-Linear-Algebra/Hogben/p/book/9781138199897) and the references therein.

## Definitions

A __weighted graph__ is an ordered triplet $G=(V,E,\omega)$, where $V=\{1,2,3,...,n\}$ is the set of __vertices__ , $E=\{(i,j)\}$ is a set of __edges__ connecting vertices, and $\omega$ is a set of __weights__ of edges. We assume $G$ is undirected.
"""

# ╔═╡ a66a6e80-1f80-11eb-15f8-71ffaa5050e7
md"
__Adjacency matrix__ of graph $G$ is the matrix $A$ defined as

$$A_{ij}=\begin{cases} 1 \quad \textrm{if}\ (i,j)\in E, \\
0\quad  \textrm{otherwise} \end{cases}.$$

__Weight matrix__ of graph $G$ is the matrix $W$ defined as

$W_{ij}=\begin{cases} \omega(e) \quad \textrm{if}\ e=(i,j)\in E, \\
0\quad  \textrm{otherwise} \end{cases}.$

__Laplacian matrix__ of graph $G$ is the matrix

$L=D-W,$

where
$D=\mathop{\mathrm{diag}}(d_1,d_2,\ldots,d_n)$ with $d_i=\sum_{k=1}^n W_{ik}$ for $i=1,\ldots,n$.

__Normalized Laplacian matrix__ is the matrix

$L_n=D^{-1/2} L D^{-1/2}\equiv D^{-1/2} (D-W) D^{-1/2}$

(__diagonally scaled $L$__).

__Incidence matrix__ of graph $G$ is the $|V|\times |E|$ matrix $I_G$. Each row of $I_G$ corresponds to a vertex of $G$ and each column corresponds to an edge of $G$.
In the column corresponding to en edge $e=(i,j)$, all elements are zero except the ones in the $i$-th and $j$-th row, which are equal to $\sqrt{\omega(e)}$ and $-\sqrt{\omega(e)}$, respectively.
"

# ╔═╡ d4cf12ce-12fb-4562-bbea-e9b85cba94c2
md"""
### Examples

Graph types and algorithms are implemented in the package [LightGraphs.jl](https://github.com/JuliaGraphs/LightGraphs.jl).

Plotting graphs is done by the packages
[GraphPlot.jl](https://github.com/JuliaGraphs/GraphPlot.jl).

As a small inconvenience, we can only plot unweighted graphs and plot weights as edge labels.
"""

# ╔═╡ 87cca1d3-b944-4d74-bf4c-1562f34cf91a
begin
	# Sources, targets and weights
	n=7
	sn=[1,1,1,2,2,3,3,3,5,5,6]
	tn=[2,3,4,4,5,4,6,7,6,7,7]
	wn=[2,3,4,7,1,3,2,1,7,3,5]
	[sn tn wn]
end

# ╔═╡ 3d7b4a01-9a7d-4111-a45d-1d31f1ef6278
begin
	# Create the graph
	G=Graph(n)
	for i=1:length(sn)
	    add_edge!(G,sn[i],tn[i])
	end
	G
end

# ╔═╡ 3b5ed1fa-75b0-43fb-a918-4c2e8f52d039
# What is the optimal bipartition?
gplot(G, nodelabel=1:n, edgelabel=wn)

# ╔═╡ 8f7d1f6c-ee4c-407f-80d9-d2abc2af949e
begin
	# We define some functions
	function WeightMatrix(src::Array,dst::Array,weights::Array)
	    n=nv(G)
	    sparse([src;dst],[dst;src],[weights;weights],n,n)
	end

	Laplacian(W::AbstractMatrix)=spdiagm(0=>vec(sum(W,dims=2)))-W

	function NormalizedLaplacian(L::AbstractMatrix)
	    D=1.0./sqrt.(diag(L))
		Diagonal(D)*L*Diagonal(D)
	end
end

# ╔═╡ b7f127d8-e8cd-4771-8091-c88f9ffe6b08
W=WeightMatrix(sn,tn,wn)

# ╔═╡ 3f5447bb-95cd-4520-b0b9-e85f33b323bd
Matrix(W)

# ╔═╡ ecb2a46e-f9dc-4e57-812f-0a2788afb202
begin
	L=Laplacian(W)
	Matrix(L)
end

# ╔═╡ 86e5f00f-13ee-4867-b0cd-187f2c782a46
Lₙ=NormalizedLaplacian(L)

# ╔═╡ cf26dbe9-3c4f-4c95-88ef-f197f360d759
Matrix(Lₙ)

# ╔═╡ 8fcd2440-4013-4b37-a747-2c0c7f6d4aeb
issymmetric(Lₙ)

# ╔═╡ 0c8d3572-3436-446f-952a-b19d0370ab38
# Let us compute the incidence matrix
function IncidenceMatrix(G::Graph, weights::Array)
    A=zeros(nv(G),ne(G))
    k=1
    for a in edges(G)
        A[a.dst,k]=sqrt.(weights[k])
        A[a.src,k]=-sqrt(weights[k])
        k+=1
    end
    A
end

# ╔═╡ 1437aadc-1c26-4eff-bb22-9cf9c15e1f5f
Iᵧ=IncidenceMatrix(G,wn)

# ╔═╡ 6529324b-49f7-4f01-b5cd-33668979b56b
md"""
## Facts

1.  $L=I_{G}I_{G}^{T}$.

2.  $L$ is symmetric PSD matrix.

3.  $L\mathbf{1}=0$ for $\mathbf{1}=[1,...,1]^{T}$, thus $0$ is an eigenvalue of $L$  and $\mathbf{1}$ is the corresponding eigenvector.

4. If $G$ has $c$ connected components, then $L$ has $c$ eigenvalues equal to $0$.

5. For every $x\in \mathbb{R}^{n}$, it holds
$x^{T}L x=\sum\limits_{i<j}W_{ij}(x_{i}-x_{j})^{2}$.

6. For every $x\in\mathbb{R}^{n}$ and $\alpha,\beta\in\mathbb{R}$, it holds
$(\alpha x+\beta \mathbf{1})^{T} L (\alpha x+\beta \mathbf{1})
=\alpha^{2} x^{T}L x$.

7. Assume that the eigenvalues of $L$ are increasingly ordered. Then,

$$
0=\lambda_1(L)\leq \lambda_2(L)\leq \cdots \leq\lambda_{n}(L)\leq
2\max\limits_{i=1,\cdots ,n}d_{i}.$$

8.  $\sigma(L_n) \subseteq [0,2]$.
"""

# ╔═╡ a9ba64a6-bd17-4b86-acb4-e7c7d8f87dc1
md"""
### Examples
"""

# ╔═╡ 6a6f5db5-7dd4-4a26-ac1a-ef3ca3a5d6be
# Fact 1
norm(L-Iᵧ*Iᵧ')

# ╔═╡ e11bc23d-03e2-4780-9488-93f0adfc2fcb
# Facts 2 and 7
issymmetric(L), eigs(L)[1], 2*maximum(diag(L))

# ╔═╡ 7192f7d0-28b9-11eb-28dd-8db0019504ca
eigen(Matrix(L))

# ╔═╡ e5890dca-368f-437e-bc45-6a395714dae2
# Fact 3
L*ones(n)

# ╔═╡ 61a978a7-b2d2-43e7-b3f7-95bf4c4f54b9
begin
	# Fact 5
	x=rand(n)
	x'*L*x, sum([W[i,j]*(x[i]-x[j])^2 for i=1:n, j=1:n])/2
end

# ╔═╡ da7f85ed-73b0-44b9-a97f-4a8ed5e9a1d1
begin
	# Fact 6
	α,β=rand(),rand()
	(α*x+β*ones(n))'*L*(α*x+β*ones(n)), α^2*x'*L*x
end

# ╔═╡ fed7f592-be03-4baf-8ace-ff9ee78fee71
# Fact 8
eigvals(Matrix(Lₙ))

# ╔═╡ 9f31928e-c933-43bb-9a5c-bf55f29737cd
md"""
# Bipartitioning

## Definitions

Let $\pi=\{V_{1},V_{2}\}$ be a partition of $V$ with $V_1,V_2\neq \emptyset$.

__Cut__ of partition $\pi$ is the sum of weights of all
edges between $V_1$ and $V_2$,

$$\mathop{\mathrm{cut}}(\pi)\equiv \mathop{\mathrm{cut}}(V_1,V_2)=\sum\limits_{{\displaystyle i\in V_{1} \atop \displaystyle j\in V_{2}}}W_{ij}.$$

__Weight__ of vertex $i\in V$ is the sum of the weights of all egdges emanating from $i$,
$\omega(i)=\sum\limits_{j=1}^{n}W_{ij}$.

__Weight__ of a subset $\bar V\subset V$ is the sum of the weights of all vertices in $\bar V$,

$\omega(\bar V)=\sum\limits_{\displaystyle i\in\bar V} \omega(i)$.

__Proportional cut__ of partition $\pi$ is

$$
\mathop{\mathrm{pcut}}(\pi)=\displaystyle\frac{\mathop{\mathrm{cut}}(\pi)}{|V_{1}|}+\frac{\mathop{\mathrm{cut}}(\pi)}{|V_{2}|}.$$

__Normalized cut__ of partition $\pi$ is

$$
\mathop{\mathrm{ncut}}(\pi)=\displaystyle\frac{\mathop{\mathrm{cut}}(\pi)}{\omega(V_{1})}+\frac{\mathop{\mathrm{cut}}(\pi)}{\omega(V_{2})}.$$
"""

# ╔═╡ 1e6d4c72-f25d-4bf5-99d6-ae0446be25dd
md"""
### Example

Consider the following partitions (all edges have unit weights):

![](https://ivanslapnicar.github.io/NumericalMathematics/files/cut2.png)

Left partition is $\pi$, right partition is $\pi'$.

|     Cut \ Partition  | $\pi$            |  $\pi'$     |
| ------- | ---------------- | ------------|
| $\mathop{\mathrm{cut}}$  |  $2$        |     $3$        |
| $\mathop{\mathrm{pcut}}$ | $\frac{2}{1}+\frac{2}{11}=2.18$|    $\frac{3}{6}+\frac{3}{6}=1$      |
| $\mathop{\mathrm{ncut}}$ | $\frac{2}{2}+\frac{2}{50}=1.04$ |  $\frac{3}{27}+\frac{3}{25}=0.23$|
"""

# ╔═╡ 02035eed-6d49-4253-a755-68a999e7e90e
md"""
## Facts

1. The informal description of the bipartitioning problem can be formulated as two problems,

$$\mathop{\textrm{arg min}}\limits_{\pi} \mathop{\mathrm{pcut}}(\pi) \quad \textrm{or} \quad \mathop{\textrm{arg min}}\limits_{\pi} \mathop{\mathrm{ncut}}(\pi).$$

The first problem favors partitions into subsets with similar numbers of vertices, while the second problem favors partitions into subsets with similar weights.

2. Both problems are NP-hard.

3. __Approximate solutions can be computed by suitable relaxations in $O(n^2)$ operations.__

4. The partition $\pi$ is defined by the vector $y$ such that

$$
y_{i}=
\begin{cases}
\frac{1}{2} & \text{for } i\in V_1 \\
-\frac{1}{2} & \text{for } i\in V_2
\end{cases}$$

The proportional cut problem can be formulated as the  __discrete proportional cut__ problem

$$
\underset{\displaystyle \big|\mathbf{y}^{T}\mathbf{1} \big|\leq \beta}
{\min\limits_{\displaystyle y_{i}\in \{-\frac{1}{2},\frac{1}{2}\}}}
\frac{1}{2}\sum_{i,j}(y_{i}-y_{j})^{2}W_{ij}.$$

Parameter $\beta$ controls the number of vertices in each subset.

5. The normalized cut problem can be formulated as the __discrete normalized cut__ problem

$$
\underset{\displaystyle \big|y^{T}D\mathbf{1} \big|\leq \beta}
{\min\limits_{\displaystyle y_{i}\in \{-\frac{1}{2},\frac{1}{2}\}}}
\frac{1}{2}\sum_{i,j}(y_{i}-y_{j})^{2}W_{ij}.$$

Parameter $\beta$ controls the weights of each subset.

6. Using the Fact 5 above, the discrete proportional cut problem can be formulated as the __relaxed proportional cut__ problem

$$
\underset{\displaystyle y^{T}y=1}{\underset{\displaystyle \big| y^{T}\mathbf{1} \big|
\leq 2\frac{\beta}{\sqrt{n}}}
{\min\limits_{\displaystyle y\in \mathbb{R}^{n}}}} y^{T}L y.$$

Similarly, the discrete normalized cut problem can be formulated as the __relaxed normalized cut__ problem

$$
\underset{\displaystyle y^{T}Dy=1}{\underset{\displaystyle \big| y^{T}D\mathbf{1}\big|
\leq \displaystyle \frac{\beta}{\sqrt{\theta n}}}{\min\limits_{\displaystyle y\in
\mathbb{R}^{n}}}}y^{T}L_n y.$$

7. __The Main Theorem.__ Let $A\in \mathbb{R}^{n\times n}$ be a symmetric matrix with eigenvalues $\lambda _{1}<\lambda _{2}<\lambda_{3}\leq \cdots \leq \lambda _{n}$ and let $v^{[1]},v^{[2]},\ldots,v^{[n]}$ be the corresponding eigenvectors. For the fixed $0\leq \alpha <1$, the solution of the problem

$$
\underset{\displaystyle y^{T}y=1}{\underset{\displaystyle \left|y^{T}v^{[1]}\right|\leq \alpha}
{\min\limits_{\displaystyle y\in \mathbb{R}^{n}}}} y^{T}Ay$$

is $y=\pm \alpha v^{[1]}\pm \sqrt{1-\alpha^{2}}v^{[2]}$.

For the proof see [D. J. Higham and M. Kibble, A Unified View of Spectral Clustering, Theorem 3.1, p. 7](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.115.1591&rep=rep1&type=pdf).

8. For $0\leq \beta <\frac{n}{2}$, the solution of the relaxed proportional cut problem is

$$
y=\pm \frac{2\beta}{n\sqrt{n}}\mathbf{1}\pm
\sqrt{1-4\frac{\beta ^{2}}{n^{2}}}v^{[2]},$$

where $v^{[2]}$ is an eigenvector corresponding to $\lambda_2(L)$. $v^{[2]}$ the __Fiedler vector__. Since the first summand carries no information, $V$ is partitioned according to the signs of the components of $v^{[2]}$:

$$
V_{1}=\{i:v^{[2]}_i <0\}, \quad V_{2}=\{i:v^{[2]}_i \geq 0\}.$$

_Notice that the value of $\beta$ is irrelevant for the solution._

9. For $0\leq \beta <\sqrt{\theta n}\left\Vert D^{\frac{1}{2}}\mathbf{1} \right\Vert _{2},$ the solution of the relaxed normalized cut problem is

$$
y=\pm \frac{\beta }{\sqrt{\theta n}\left\Vert
D^{\frac{1}{2}} \mathbf{1}\right\Vert _{2}^{2}}\mathbf{1}\pm
\sqrt{1-\frac{\beta ^{2}}{ \theta n\left\Vert
D^{\frac{1}{2}}\mathbf{1}\right\Vert _{2}^{2}}}D^{-\frac{1
}{2}} v_n^{[2]},$$

where $v_n^{[2]}$ is an eigenvector corresponding to $\lambda_2(L_n)$. $V_n$ is partitioned according to the signs of the components of $v_n^{[2]}$, as above.

10. Neither of the relaxed algorithms is guaranteed to solve exactly the true proportional / normalized cut problem. However, the computed solutions are in the right direction. Whether to use proportional or normalized cut formulation, depends upon the specific problem.
"""

# ╔═╡ c188b7cf-03e7-40b3-bf90-773978c7b410
# Voila!
eigs(L,nev=2,which=:SM, v0=ones(n))

# ╔═╡ f60b51b6-79c3-4a3f-8bfb-8f198cdb92a8
# For the normalized cut
eigs(Lₙ,nev=2,which=:SM, v0=ones(n))

# ╔═╡ bfac7a55-0d82-4549-b999-7db3272eaa20
md"""
### Concentric rings

A __complete graph__ has edges connecting each pair of vertices.

To a set of points $X=\{x_{1},x_{2},\cdots ,x_{m}\}$ , where $x_{i}\in\mathbb{R}^{n}$, we assign a weighted complete graph $G=(V,E)$ with $m$ vertices, where the vertex $j\in V$ corresponds to the point $x_j\in X$.

The main idea is to assign weight of an edge $e=(i,j)$ which reflects the distance between $x_i$ and $x_j$, something like $\omega(e)=\displaystyle\frac{1}{\mathop{\mathrm{dist}}(x_i,x_j)}$.

However, this has to be implemented with care. For example, using simple Euclidean distance yield the same results as the function `kmeans()`. In this example we use Gaussian kernel, that is

$$
\omega(e)=e^{\displaystyle -\|x_i-x_j\|_2^2/\sigma^2},$$

where the choice of $\sigma$ is based on experience.

The computation of various distances is implemented in the package [Distances.jl](https://github.com/JuliaStats/Distances.jl).

We will construct the Laplace matrix directly.
"""

# ╔═╡ 80ed91ad-6f9a-4f12-a1ee-bc8dc6b58168
begin
	# Two concentric circles
	k=2
	# Center
	Random.seed!(541)
	# center=[rand(-5:5),rand(-5:5)]
	center=[0,0]
	# Radii
	radii=randperm(10)[1:k]
	# Number of points in circles
	sizes=rand(1000:2000,k)
	center,radii,sizes
end

# ╔═╡ 2d8689c9-1b9e-4360-865c-23c2c8c5dcb9
begin
	# Generate points
	X=Array{Float64}(undef,2,sum(sizes))
	csizes=cumsum(sizes)
	# Random angles
	ϕ=2*π*rand(sum(sizes))
	for i=1:csizes[1]
		X[:,i]=center+radii[1]*[cos(ϕ[i]);sin(ϕ[i])] + (rand(2).-0.5)/50
	end
	for j=2:k
		for i=csizes[j-1]+1:csizes[j]
			X[:,i]=center+radii[j]*[cos(ϕ[i]);sin(ϕ[i])] + (rand(2).-0.5)/50
		end
	end
	scatter(X[1,:],X[2,:],title="Concentric rings", aspect_ratio=1,label="Points")
end

# ╔═╡ fa4ae713-a0b2-46ce-b5f7-558b26d82019
# Weight matrix
W₁=1 ./pairwise(SqEuclidean(),X)

# ╔═╡ 760e8ca2-596f-435f-befb-140b3d204f7f
begin
	# Laplacian matrix
	m=csizes[end]
	for i=1:m
	    W₁[i,i]=0
	end
	L₁=Diagonal(vec(sum(W₁,dims=2)))-W₁
	# Check Fact 3
	norm(L₁*ones(m))
end

# ╔═╡ 75956488-76ab-4262-a32f-0c53e1b3dc17
# Notice λ₁=0
E=eigs(L₁,nev=2,which=:SM, v0=ones(m))

# ╔═╡ a55260dd-f05d-4a65-9fdf-cb3fd4638857
begin
	# Define clusters
	C=ones(Int64,m)
	C[findall(E[2][:,2].>0)].=2
	C
end

# ╔═╡ ba164ff3-5ef4-4480-8af2-9dd9a99b1e8d
# Yet another plotting function
function plotKpartresult(C::Vector,X::Array)
	scatter(aspect_ratio=1)
    k=maximum(C)
    for j=1:k
        scatter!(X[1,findall(C.==j)],X[2,findall(C.==j)],label="Cluster $j")
    end
	scatter!(aspect_ratio=1)
end

# ╔═╡ 46de28bc-73ed-495c-a2e0-f565c1ce5651
plotKpartresult(C,X)

# ╔═╡ e004e4bf-ca04-4b41-947b-9f1ecd058abb
md"""
This is the same partitioning as obtained by `kmeans()`. Let us try Gaussian kernel. A rule of thumb is: if rings are close, use $\sigma<1$, if rings are apart, use $\sigma>1$.
"""

# ╔═╡ bb9d4d98-3aff-439e-b8a3-8347cf345839
begin
	σ=0.7 # 0.1
	W₂=exp.(-pairwise(SqEuclidean(),X)/σ^2)-I
	L₂=Diagonal(vec(sum(W₂,dims=2)))-W₂
	E₂=eigs(L₂,nev=2,which=:SM, v0=ones(m))
	C₂=ones(Int64,m)
	C₂[findall(E₂[2][:,2].>0)].=2
	plotKpartresult(C₂,X)
end

# ╔═╡ 012fa54a-5171-4980-a0ba-474a22c25981
md"""
# Recursive bipartitioning

## Definitions

Let $G=(V,E)$ be a weighted graph with weights $\omega$.

Let $\pi_k =\{V_{1},V_{2},...,V_{k}\}$ be a $k$-partition of $V$, with $V_i\neq \emptyset$ for $i=1,\ldots,k$.

The previous definition of $cut(\pi)\equiv cut(\pi_2)$ extends naturally to $k$-partition.

A __cut__ of a partition $\pi_k$ is

$$
\mathop{\mathrm{cut}}(\pi_k)=\sum\limits_{\displaystyle i<j} \mathop{\mathrm{cut}}(V_{i},V_{j}),$$

where $\mathop{\mathrm{cut}}(V_{i},V_{j})$ is interpreted as a cut of the bipartition of the subgraph of $G$ with vertices $V_1\cup V_2$.

__Proportional cut__ of a partition $\pi_k$ is

$$
\mathop{\mathrm{pcut}}(\pi_k)=\underset{i<j}{\sum\limits_{i,j=1}^{k}} \left(
\frac{\mathop{\mathrm{cut}}(V_{i},V_{j})}{|V_{i}|}+\frac{\mathop{\mathrm{cut}}(V_{i},V_{j})}{|V_{j}|}\right) =
\sum_{i=1}^{k}\frac{\mathop{\mathrm{cut}}(V_{i},V\backslash V_{i})}{|V_{i}|}.$$

__Normalized cut__ of a partition $\pi_k$ is

$$
\mathop{\mathrm{ncut}}(\pi_k)=\underset{i<j}{\sum\limits_{i,j=1}^{k}} \left(
\frac{\mathop{\mathrm{cut}}(V_{i},V_{j})}{\omega(V_{i})}+\frac{\mathop{\mathrm{cut}}(V_{i},V_{j})}{\omega(V_{j})}\right) =
\sum_{i=1}^{k}\frac{\mathop{\mathrm{cut}}(V_{i},V\backslash V_{i})}{ \omega(V_{i})}.$$

## Facts

If we want to cluster vertices of graph $G=(V,E)$ into $k$ clusters, we can apply the following recursive algorithm:

1. __Initialization.__ Compute the bipartition $\pi=\{V_{1},V_{2}\}$ of $V$. Set the counter $c=2$.

2. __Recursion.__ While $c<k$ repeat:

    1. Compute the bipartition of each subset of $V$.

    2. Among all $(c+1)$-partitions, choose the one with the smallest $\mathop{\mathrm{pcut}}(\pi_{c+1})$ or $\mathop{\mathrm{ncut}}(\pi_{c+1})$, respectively.

    3. Set $c=c+1$.

3. __Stop.__

There is no guarantee for optimality of this algorithm. Clearly, the optimal $k$-partiton may be a subpartition of one of the discarded partitions.
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
GraphPlot = "a2cc645c-3eea-5389-862e-a155d0052231"
LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[compat]
Arpack = "~0.5.3"
Distances = "~0.10.4"
GraphPlot = "~0.4.4"
LightGraphs = "~1.3.5"
Plots = "~1.21.3"
PlutoUI = "~0.7.9"
"""

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

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

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

[[ArnoldiMethod]]
deps = ["LinearAlgebra", "Random", "StaticArrays"]
git-tree-sha1 = "f87e559f87a45bece9c9ed97458d3afe98b1ebb9"
uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
version = "0.1.0"

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

[[Arpack_jll]]
deps = ["Libdl", "OpenBLAS_jll", "Pkg"]
git-tree-sha1 = "e214a9b9bd1b4e1b4f15b22c0994862b66af7ff7"
uuid = "68821587-b530-5797-8361-c406ea357684"
version = "3.5.0+3"

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

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

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

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

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

[[ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "32a2b8af383f11cbb65803883837a149d10dfe8a"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.10.12"

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

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

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

[[Compose]]
deps = ["Base64", "Colors", "DataStructures", "Dates", "IterTools", "JSON", "LinearAlgebra", "Measures", "Printf", "Random", "Requires", "Statistics", "UUIDs"]
git-tree-sha1 = "c6461fc7c35a4bb8d00905df7adafcff1fe3a6bc"
uuid = "a81c6b42-2e10-5240-aca2-a61377ecd94b"
version = "0.9.2"

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

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

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

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

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

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

[[Distances]]
deps = ["LinearAlgebra", "Statistics", "StatsAPI"]
git-tree-sha1 = "9f46deb4d4ee4494ffb5a40a27a2aced67bdd838"
uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
version = "0.10.4"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[GraphPlot]]
deps = ["ArnoldiMethod", "ColorTypes", "Colors", "Compose", "DelimitedFiles", "LightGraphs", "LinearAlgebra", "Random", "SparseArrays"]
git-tree-sha1 = "dd8f15128a91b0079dfe3f4a4a1e190e54ac7164"
uuid = "a2cc645c-3eea-5389-862e-a155d0052231"
version = "0.4.4"

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

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

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

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

[[Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[LightGraphs]]
deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
git-tree-sha1 = "432428df5f360964040ed60418dd5601ecd240b6"
uuid = "093fc24a-ae57-5d10-9952-331d41423f4d"
version = "1.3.5"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "438d35d2d95ae2c5e8780b330592b6de8494e779"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.3"

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

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

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

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

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

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[SimpleTraits]]
deps = ["InteractiveUtils", "MacroTools"]
git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
version = "0.9.4"

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

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

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

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

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

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

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

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

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# ╔═╡ Cell order:
# ╠═0cd92d78-1da6-435b-ad8e-60af29502e1d
# ╠═7ef4b0b3-9d0a-4356-bd16-14487f14e0ab
# ╟─cb61f761-b2ee-4818-9149-f85d95b76a1e
# ╟─5926c79a-6038-4dbc-b49a-5b05d02250a0
# ╟─a66a6e80-1f80-11eb-15f8-71ffaa5050e7
# ╟─d4cf12ce-12fb-4562-bbea-e9b85cba94c2
# ╠═87cca1d3-b944-4d74-bf4c-1562f34cf91a
# ╠═3d7b4a01-9a7d-4111-a45d-1d31f1ef6278
# ╠═3b5ed1fa-75b0-43fb-a918-4c2e8f52d039
# ╠═8f7d1f6c-ee4c-407f-80d9-d2abc2af949e
# ╠═b7f127d8-e8cd-4771-8091-c88f9ffe6b08
# ╠═3f5447bb-95cd-4520-b0b9-e85f33b323bd
# ╠═ecb2a46e-f9dc-4e57-812f-0a2788afb202
# ╠═86e5f00f-13ee-4867-b0cd-187f2c782a46
# ╠═cf26dbe9-3c4f-4c95-88ef-f197f360d759
# ╠═8fcd2440-4013-4b37-a747-2c0c7f6d4aeb
# ╠═0c8d3572-3436-446f-952a-b19d0370ab38
# ╠═1437aadc-1c26-4eff-bb22-9cf9c15e1f5f
# ╟─6529324b-49f7-4f01-b5cd-33668979b56b
# ╟─a9ba64a6-bd17-4b86-acb4-e7c7d8f87dc1
# ╠═6a6f5db5-7dd4-4a26-ac1a-ef3ca3a5d6be
# ╠═e11bc23d-03e2-4780-9488-93f0adfc2fcb
# ╠═7192f7d0-28b9-11eb-28dd-8db0019504ca
# ╠═e5890dca-368f-437e-bc45-6a395714dae2
# ╠═61a978a7-b2d2-43e7-b3f7-95bf4c4f54b9
# ╠═da7f85ed-73b0-44b9-a97f-4a8ed5e9a1d1
# ╠═fed7f592-be03-4baf-8ace-ff9ee78fee71
# ╟─9f31928e-c933-43bb-9a5c-bf55f29737cd
# ╟─1e6d4c72-f25d-4bf5-99d6-ae0446be25dd
# ╟─02035eed-6d49-4253-a755-68a999e7e90e
# ╠═c188b7cf-03e7-40b3-bf90-773978c7b410
# ╠═f60b51b6-79c3-4a3f-8bfb-8f198cdb92a8
# ╟─bfac7a55-0d82-4549-b999-7db3272eaa20
# ╠═80ed91ad-6f9a-4f12-a1ee-bc8dc6b58168
# ╠═2d8689c9-1b9e-4360-865c-23c2c8c5dcb9
# ╠═fa4ae713-a0b2-46ce-b5f7-558b26d82019
# ╠═760e8ca2-596f-435f-befb-140b3d204f7f
# ╠═75956488-76ab-4262-a32f-0c53e1b3dc17
# ╠═a55260dd-f05d-4a65-9fdf-cb3fd4638857
# ╠═ba164ff3-5ef4-4480-8af2-9dd9a99b1e8d
# ╠═46de28bc-73ed-495c-a2e0-f565c1ce5651
# ╟─e004e4bf-ca04-4b41-947b-9f1ecd058abb
# ╠═bb9d4d98-3aff-439e-b8a3-8347cf345839
# ╟─012fa54a-5171-4980-a0ba-474a22c25981
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_path: Lectures/C2 Spectral Graph Bipartitioning.jl


text:
using TSVD
using LinearAlgebra
using ToeplitzMatrices
using ARMA

make_mpr(data, basis) = pinv(basis)*data
make_time_domain(mpr, basis) = basis*mpr
function make_std_residuals(data, basis)
    mpr = make_mpr(data,basis) # model pulse reduced
    td = make_time_domain(mpr,basis)
    std_residuals = std(data-td, dims=1)[1,:]
end

function getall(ljh, maxrecords=typemax(Int))
    records = collect(ljh[1:min(Int(maxrecords),length(ljh))])
    pulses = Array{Float32}(undef, ljh.record_nsamples,length(records))
    for (i,record) in enumerate(records)
        pulses[:,i]=record.data
    end
    pulses
end

function evenly_distributed_inds(inds,n_wanted)
    frac_keep = n_wanted/length(inds)
    keep_each_n = max(1,floor(Int,1/frac_keep))
    inds[1:keep_each_n:length(inds)]
end

function choose_new_train_inds(residuals, train_inds, frac_keep)
    n_keep = round(Int, frac_keep*length(train_inds))
    sort(train_inds, by = i -> residuals[i])[1:n_keep]
end

function train_loop(data, n_pulses_for_train, n_basis, n_presamples, n_loop, frac_keep_per_loop, make_basis)
    last_train_inds = train_inds = evenly_distributed_inds(1:size(data,2),n_pulses_for_train)
    # Use only a 3-Dimensional basis the first time through. This helps reject pileup before it
    # can sneak into the SVD.
    nb = min(n_basis, 3)
    for i=1:n_loop
        basis, singular_values = make_basis(data[:,train_inds],nb)
        residuals = make_std_residuals(data,basis)
        last_train_inds = train_inds
        train_inds = choose_new_train_inds(residuals, train_inds, frac_keep_per_loop)
        if i==n_loop # variables created in loop are not availble outside the loop
            return basis, residuals, last_train_inds, train_inds, singular_values
        end
        nb = n_basis
    end
end

function TSVD_tsvd(data_train::Matrix{<:AbstractFloat},n_basis)
    U, S, V = TSVD.tsvd(data_train, n_basis)
    U, S
end

"""    TSVD_tsvd_mass3(data_train::Matrix{<:AbstractFloat}, n_basis, n_presamples, noise_model::ARMA.ARMAModel, noise_solver::ARMA.ARMASolver)

This function is called as the first step towards creating a basis with `n_basis` elements (minimum 3).
Elements will be 1) constant, 2) derivative-like, 3) pulse-like with first `n_presamples` equal to zero.
Elements 4 and higher are based on SVD of the noise whitened residuals after projecting into a basis of
just the first 3 elements. This ensures that adding additional components will not change the projector
for the first 3 components, and therefore the energy resolution obtained with only the first 3 projectors
will be independing of number of components.

The basis is paired with projectors, which allow you to project into the basis while minimizing
the noise whitened residuals (aka the mahalonobis distance).

This function returns `U,S` much like a truncated singular value decomposition, but such that the
basis caluclated from `U,S` will meet the description above."""
function TSVD_tsvd_mass3(data_train::Matrix{<:AbstractFloat}, n_basis, n_presamples,
        noise_model::ARMA.ARMAModel, noise_solver::ARMA.ARMASolver)
    @assert n_basis>=3 "mean, derivative and average pulse are 3 components, must request at least 3 components"
    # Average pulse is the pretrigger-mean-subtracted average pulse, rescaled to have a maximum value of 1.0
    average_pulse = mean(data_train, dims=2)[:]
    if n_presamples > 0
        average_pulse .-= mean(average_pulse[1:n_presamples])
    end
    average_pulse[1:n_presamples] .= 0.0
    average_pulse /= maximum(abs.(average_pulse))
    # Calculate the derivative like component; assume clean baseline, so low derivative at start
    derivative_like = [0.0;diff(average_pulse)]
    constant_component = ones(size(data_train,1))
    mass3_basis = hcat(constant_component,derivative_like,average_pulse)
    if n_basis == 3
        return mass3_basis, [NaN,NaN,NaN]
    end
    projectors3, _ = computeprojectors(mass3_basis, noise_model)
    mpr = projectors3 * data_train # model pulse reduced
    td = mass3_basis * mpr
    data_residual = data_train .- td
    # Whiten residual before taking the TSVD
    white_residual = ARMA.whiten(noise_solver, data_residual)
    Uwhite, S, _ = TSVD.tsvd(white_residual, n_basis-3)
    # Unwhiten and renormalize U before combining
    Uunnorm = ARMA.unwhiten(noise_solver, Uwhite)
    U = Uunnorm ./ sqrt.(sum(Uunnorm.^2, dims=1))
    U_combined = hcat(mass3_basis,U)
    S_combined = vcat([NaN,NaN,NaN],S) # first 3 elemnts are not from svd, have no meaningful singular value
    U_combined, S_combined
end

"""
    TSVD_tsvd_mass3(data_train::Matrix{<:AbstractFloat}, n_basis, n_presamples, autocorr::Vector)

Create mass3+SVD as the other method does, but use the noise autocorrelation (`autocorr`)
instead of using a noise model. This should mirror the MASS behavior.
"""
function TSVD_tsvd_mass3(data_train::Matrix{<:AbstractFloat}, n_basis, n_presamples, autocorr::Vector)
    @assert n_basis>=3 "mean, derivative and average pulse are 3 components, must request at least 3 components"
    # Average pulse is the pretrigger-mean-subtracted average pulse, rescaled to have a maximum value of 1.0
    average_pulse = mean(data_train, dims=2)[:]
    if n_presamples > 0
        average_pulse .-= mean(average_pulse[1:n_presamples])
    end
    average_pulse[1:n_presamples] .= 0.0
    average_pulse /= maximum(abs.(average_pulse))
    # Calculate the derivative like component; assume clean baseline, so low derivative at start
    derivative_like = [0.0;diff(average_pulse)]
    constant_component = ones(size(data_train,1))
    mass3_basis = hcat(constant_component,derivative_like,average_pulse)
    if n_basis == 3
        return mass3_basis, [NaN,NaN,NaN]
    end
    projectors3, _ = computeprojectors(mass3_basis, autocorr)
    mpr = projectors3 * data_train # model pulse reduced
    td = mass3_basis * mpr
    data_residual = data_train .- td
    # Whiten residual before taking the TSVD
    R = SymmetricToeplitz(autocorr)
    Winv = cholesky(R).L  # The inverse of Winv will whiten a data record.
    white_residual = Winv \ data_residual
    Uwhite, S, _ = TSVD.tsvd(white_residual, n_basis-3)
    # Unwhiten and renormalize U before combining
    Uunnorm = Winv * Uwhite
    U = Uunnorm ./ sqrt.(sum(Uunnorm.^2, dims=1))
    U_combined = hcat(mass3_basis, U)
    S_combined = vcat([NaN,NaN,NaN], S) # first 3 elemnts are not from svd, have no meaningful singular value
    U_combined, S_combined
end

function full_svd(data_train, n_basis)
    U,S,V = svd(data_train)
    U[:,1:n_basis], S[1:n_basis]
end
tsvd_dict = Dict(
    "TSVD" => TSVD_tsvd,
    "full" => full_svd,
    "noisemass3" => TSVD_tsvd_mass3,
    "TSVDmass3" => TSVD_tsvd_mass3)

struct SVDBasis <: AbstractBasisAnalyzer
    basis::Array{Float32,2} # unwhitened time domain representation, tall: basis*mpr = reconsituted pulse
    projectors::Array{Float32,2} # whitened time domain representation, wide: projectors*data = mpr
    projector_covariance::Array{Float32,2}
    noise_result::NoiseAnalysis.NoiseResult
end
(a::SVDBasis)(record) = record2dataproduct(a,record)

struct SVDBasisWithCreationInfo
    svdbasis::SVDBasis
    singular_values::Vector{Float32}
    example_pulses::Array{UInt16,2}
    std_residuals_of_example_pulses::Vector{Float32}
    percentiles_of_sample_pulses::Vector{Float32}
    n_loop::Int
    noise_model_file::String
    pulse_file::String
    std_residuals::Vector{Float32}
    tsvd_method::String
    channel_number::Int
    noise_std_dev::Float64
end

function create_basis_one_channel(ljh, noise_result, frac_keep, n_loop, n_pulses_for_train, n_basis,
    tsvd_method_string)
    data = getall(ljh, ceil(Int,n_pulses_for_train/frac_keep))
    create_basis_one_channel(data, noise_result, frac_keep, n_loop,
        n_pulses_for_train, n_basis, tsvd_method_string,
        LJH.first_rising_sample(ljh[1]),
        LJH.filename(ljh), LJH.channel(ljh))
end

function create_basis_one_channel(data::Matrix{<:AbstractFloat}, noise_result, frac_keep, n_loop, n_pulses_for_train, n_basis,
    tsvd_method_string, n_presamples, datasource_filename, datasource_channel)
    tsvd_func = tsvd_dict[tsvd_method_string]
    if tsvd_method_string == "TSVDmass3"
        noise_solver = ARMA.ARMASolver(noise_result.model, size(data, 1))
        function tsvd_closure(data, n_basis)
            TSVD_tsvd_mass3(data, n_basis, n_presamples, noise_result.model, noise_solver)
        end
        tsvd_func = tsvd_closure
    elseif tsvd_method_string == "noisemass3"
        function tsvd_closure2(data, n_basis)
            TSVD_tsvd_mass3(data, n_basis, n_presamples, noise_result.autocorr)
        end
        tsvd_func = tsvd_closure2
    end

    std_noise = sqrt(noise_result.autocorr[1]) # the first element in the autocorrelation is the variance
    frac_keep_per_loop = exp(log(frac_keep)/n_loop)
    basis, residual_stds, last_train_inds, train_inds, singular_values = train_loop(data,
        n_pulses_for_train, n_basis, n_presamples, n_loop, frac_keep_per_loop, tsvd_func)
    if tsvd_method_string == "noisemass3"
        projectors, pcovar = computeprojectors(basis,noise_result.autocorr)
    else
        projectors, pcovar = computeprojectors(basis,noise_result.model)
    end
    svdbasis = SVDBasis(
        basis,
        projectors,
        pcovar,
        noise_result)
    sortinds = sortperm(residual_stds)
    percentiles = [10:10:90; 91:99]
    percentile_indicies = sortinds[[round(Int,(p/100)*length(residual_stds)) for p in percentiles]]
    example_pulses = round.(UInt16,data[:,percentile_indicies])
    svdbasis_with_info = SVDBasisWithCreationInfo(
        svdbasis, singular_values, example_pulses,
        residual_stds[percentile_indicies], percentiles,
        n_loop, noise_result.datasource, datasource_filename,
        residual_stds, tsvd_method_string, datasource_channel,
        √(svdbasis.noise_result.autocorr[1])
        )
    return svdbasis,svdbasis_with_info
end

function hdf5save(g::HDF5.DataFile, svdbasis::SVDBasis)
    g["basis"]=svdbasis.basis
    g["projectors"]=svdbasis.projectors
    g["projector_covariance"]=svdbasis.projector_covariance
    Pope.NoiseAnalysis.hdf5save(g_create(g,"noise_result"),svdbasis.noise_result)
end

function hdf5save(g::HDF5.DataFile, x::SVDBasisWithCreationInfo)
    hdf5save(g_create(g,"svdbasis"),x.svdbasis)
    g["singular_values"]=x.singular_values
    g["example_pulses"]=x.example_pulses
    g["std_residuals_of_example_pulses"]=x.std_residuals_of_example_pulses
    g["percentiles_of_sample_pulses"]=x.percentiles_of_sample_pulses
    g["n_loop"]=x.n_loop
    g["noise_model_file"]=x.noise_model_file
    g["pulse_file"]=x.pulse_file
    g["std_residuals"]=x.std_residuals
    g["tsvd_method"]=x.tsvd_method
    g["channel_number"]=x.channel_number
    g["noise_std_dev"]=x.noise_std_dev
end

function hdf5load(T::Type{SVDBasis},g::HDF5.DataFile)
    SVDBasis(
    read(g["basis"]),
    read(g["projectors"]),
    read(g["projector_covariance"]),
    NoiseAnalysis.hdf5load(g["noise_result"])
    )
end

function hdf5load(T::Type{SVDBasisWithCreationInfo},g::HDF5.DataFile)
    SVDBasisWithCreationInfo(
    hdf5load(SVDBasis,g["svdbasis"]),
    read(g["singular_values"]),
    read(g["example_pulses"]),
    read(g["std_residuals_of_example_pulses"]),
    read(g["percentiles_of_sample_pulses"]),
    read(g["n_loop"]),
    read(g["noise_model_file"]),
    read(g["pulse_file"]),
    read(g["std_residuals"]),
    read(g["tsvd_method"]),
    read(g["channel_number"]),
    read(g["noise_std_dev"])
    )
end

function make_basis_one_channel(outputh5, ljhname, noise_filename, frac_keep, n_loop,
    n_pulses_for_train, n_basis, tsvd_method)
    ljh = LJH.ljhopen(ljhname)
    @show LJH.channel(ljh)
    @show ljhname
    @show ljh
    noise_result = NoiseAnalysis.hdf5load(noise_filename,LJH.channel(ljh))
    svdbasis, svdbasiswithcreationinfo = create_basis_one_channel(ljh, noise_result,
        frac_keep,
        n_loop,
        n_pulses_for_train,
        n_basis,
        tsvd_method)
    hdf5save(g_create(outputh5,"$(LJH.channel(ljh))"), svdbasiswithcreationinfo)
end

function make_basis_all_channel(outputh5, ljhdict, noise_filename, frac_keep, n_loop,
    n_pulses_for_train, n_basis, tsvd_method)
    noise_h5 = h5open(noise_filename,"r")
    noise_channels = sort!(parse.(Int,names(noise_h5)))
    ljh_channels = keys(ljhdict)
    bothchannels = sort!(intersect(noise_channels,ljh_channels))
    @show noise_channels
    @show ljh_channels
    @show bothchannels
    for channel_number in bothchannels
        println("\nmaking basis for channel $channel_number")
        ljhname = ljhdict[channel_number]
        try
            make_basis_one_channel(outputh5, ljhname, noise_filename, frac_keep, n_loop,
                n_pulses_for_train, n_basis, tsvd_method)
        catch ex
            println("channel $channel_number failed:")
            println(ex)
        end
    end
    println("Channels in noise_file but not in ljh: $(setdiff(noise_channels,ljh_channels))")
    println("Channels in ljh but not in noise_file: $(setdiff(ljh_channels,noise_channels))")
end

ext: jl


text:
# Default values for the integrator

"Default values for the integrator"
module defaults
import Ogre: bar, K, R_earth

# These were mostly used in stage 1 for the interior models only
"Surface pressure"
const P_surf = 1bar
"Surface temperature"
const T_surf = 300K
"Mass fractions of the different layers"
const mass_fractions = [1.]
"Total number of points to use in the integration"
const total_points = 500
"Radius search bracket"
const R_bracket = [0., 10.] .* R_earth
end


text:
module Venti

using Genie, Logging, LoggingExtras

function main()
  Base.eval(Main, :(const UserApp = Venti))

  Genie.genie(; context = @__MODULE__)

  Base.eval(Main, :(const Genie = Venti.Genie))
  Base.eval(Main, :(using Genie))
end

end


text:
"""
    rotation_matrix_zyx(α1, α2, α3)

Calculate the rotation matrix ``R_{ab}`` that transforms a vector in b frame to a frame
(``v_{a} = R_{ab} v_{b}``) given the rotations (α1, α2, α3) (rad).

Frame b is obtained from frame a performing three intrinsic rotations of magnitude
α1, α2 and α3 in ZYX order.
"""
function rotation_matrix_zyx(α1, α2, α3)

    sα1, cα1 = sin(α1), cos(α1)
    sα2, cα2 = sin(α2), cos(α2)
    sα3, cα3 = sin(α3), cos(α3)

    R = @SMatrix [
        cα2 * cα1   (sα3 * sα2 * cα1 - cα3 * sα1)   (cα3 * sα2 * cα1 + sα3 * sα1)
        cα2 * sα1   (sα3 * sα2 * sα1 + cα3 * cα1)   (cα3 * sα2 * sα1 - sα3 * cα1)
        -sα2        sα3 * cα2                       cα3 * cα2
    ]

    return R
end

"""
    rotation_matrix_zyx(q0, q1, q2, q3)

Calculate the rotation matrix ``R_{ab}`` that transforms a vector in b frame to a frame
(``v_{a} = R_{ab} v_{b}``) given the quaternions ``q_0, q_1, q_2, q_3``.
"""
function rotation_matrix_zyx(q0, q1, q2, q3)

    q02, q12, q22, q32 = q0*q0, q1*q1, q2*q2, q3*q3

    R = @SMatrix [
        (q02+q12-q22-q32)    2*(q1*q2 - q0*q3)    2*(q1*q3 + q0*q2)
        2*(q1*q2 + q0*q3)    (q02-q12+q22-q32)    2*(q2*q3 - q0*q1)
        2*(q1*q3 - q0*q2)    2*(q2*q3 + q0*q1)    (q02-q12-q22+q32)
    ]

    return R
end

"""
    quaternions(ψ, θ, ϕ)

Calculate quaternion representation given the Euler angles (ψ, θ, ϕ) (rad).
"""
function quaternions(ψ, θ, ϕ)

    s_ψ2, c_ψ2 = sin(ψ/2), cos(ψ/2)
    s_θ2, c_θ2 = sin(θ/2), cos(θ/2)
    s_ϕ2, c_ϕ2 = sin(ϕ/2), cos(ϕ/2)

    q0 = c_ψ2*c_θ2*c_ϕ2 + s_ψ2*s_θ2*s_ϕ2
    q1 = c_ψ2*c_θ2*s_ϕ2 - s_ψ2*s_θ2*c_ϕ2
    q2 = c_ψ2*s_θ2*c_ϕ2 + s_ψ2*c_θ2*s_ϕ2
    q3 = s_ψ2*c_θ2*c_ϕ2 - c_ψ2*s_θ2*s_ϕ2

    return @SVector [q0, q1, q2, q3]
end

"""
    euler_angles(q0, q1, q2, q3)

Calculate Euler angles (ψ, θ, ϕ) (rad) given the quaternions ``q_0, q_1, q_2, q_3``.
"""
function euler_angles(q0, q1, q2, q3)

    ψ = atan(2 * (q1*q2 + q0*q3), q0*q0 + q1*q1 - q2*q2 - q3*q3)
    θ = asin(-2 * (q1*q3 - q0*q2))
    ϕ = atan(2 * (q2*q3 + q0*q1), q0*q0 - q1*q1 - q2*q2 + q3*q3)

    return @SVector [mod2pi(ψ), θ, ϕ]
end

"""
    body2horizon(x, y, z, ψ, θ, ϕ)

Transform the vector coordintes (x, y, z) given in body axis to local horizon given the
Euler angles (ψ, θ, ϕ) (rad).
"""
function body2horizon(x, y, z, ψ, θ, ϕ)
    v = @SVector [x, y, z]
    rv = rotation_matrix_zyx(ψ, θ, ϕ) * v
    return rv
end

"""
    horizon2body(x, y, z, ψ, θ, ϕ)

Transform the vector coordintes (x, y, z) given in local horizon axis to body given the
Euler angles (ψ, θ, ϕ) (rad).
"""
function horizon2body(x, y, z, ψ, θ, ϕ)
    v = @SVector [x, y, z]
    rv = transpose(rotation_matrix_zyx(ψ, θ, ϕ)) * v
    return rv
end

"""
    body2horizon(x, y, z, q0, q1, q2, q3)

Transform the vector coordintes (x, y, z) given in body axis to local horizon given the
quaternions ``q_0, q_1, q_2, q_3``.
"""
function body2horizon(x, y, z, q0, q1, q2, q3)
    v = @SVector [x, y, z]
    rv = rotation_matrix_zyx(q0, q1, q2, q3) * v
    return rv
end

"""
    horizon2body(x, y, z, q0, q1, q2, q3)

Transform the vector coordintes (x, y, z) given in local horizon axis to body given the
quaternions ``q_0, q_1, q_2, q_3``.
"""
function horizon2body(x, y, z, q0, q1, q2, q3)
    v = @SVector [x, y, z]
    rv = transpose(rotation_matrix_zyx(q0, q1, q2, q3)) * v
    return rv
end

"""
    wind2body(x, y, z, α, β)

Transform the vector coordintes (x, y, z) given in wind axis to body given the angle of
attack (α) and the angle of sideslip (β) (rad).
"""
function wind2body(x, y, z, α, β)
    v = @SVector [x, y, z]
    rv = transpose(rotation_matrix_zyx(-β, α, 0)) * v
    return rv
end

"""
    body2wind(x, y, z, α, β)

Transform the vector coordintes (x, y, z) given in body axis to wind given the angle of
attack (α) and the angle of sideslip (β) (rad).
"""
function body2wind(x, y, z, α, β)
    v = @SVector [x, y, z]
    rv = rotation_matrix_zyx(-β, α, 0) * v
    return rv
end

"""
    wind2horizon(x, y, z, χ, γ, μ)

Transform the vector coordintes (x, y, z) given in wind axis to local horizon given the
velocity angles (χ, γ, μ) (rad).
"""
function wind2horizon(x, y, z, χ, γ, μ)
    v = @SVector [x, y, z]
    rv = rotation_matrix_zyx(χ, γ, μ) * v
    return rv
end

"""
    horizon2wind(x, y, z, χ, γ, μ)

Transform the vector coordintes (x, y, z) given in local horizon axis to wind given the
velocity angles (χ, γ, μ) (rad).
"""
function horizon2wind(x, y, z, χ, γ, μ)
    v = @SVector [x, y, z]
    rv = transpose(rotation_matrix_zyx(χ, γ, μ)) * v
    return rv
end

"""
    ecef2horizon(x, y, z, lat, lon)

Transform the vector coordintes (x, y, z) given in ECEF (Earth Fixed Earth Centered)
coordinates to local horizon coordinates using geodetic latitude and longitude (rad).
"""
function ecef2horizon(x, y, z, lat, lon)
    v = @SVector [x, y, z]
    rv = transpose(rotation_matrix_zyx(lon, -lat - π/2, 0)) * v
end

"""
    horizon2ecef(x, y, z, lat, lon)

Transform the vector coordintes (x, y, z) given in local horizon axis to ECEF (Earth Fixed
Earth Centered) using geodetic latitude and longitude (rad).
"""
function horizon2ecef(x, y, z, lat, lon)
    v = @SVector [x, y, z]
    rv = rotation_matrix_zyx(lon, -lat - π/2, 0) * v
end


text:
############################################################################
### Types
############################################################################

struct SNLLS{N1, N2, N3, I1, I2, I3, I6, M1, M2, M3, M4, N4, I4, M5, I5, V1} <: SemImply

    q_directed::N1
    q_undirected::N2
    size_σ::N3

    A_indices_linear::I1
    A_indices_cartesian::I2
    S_indices::I3
    σ_indices::I6

    A_pre::M1
    I_A::M2
    G::M3
    ∇G::M4

    q_mean::N4
    M_indices::I4
    G_μ::M5
    G_μ_indices::I5

    start_val::V1

end

############################################################################
### Constructors
############################################################################

function SNLLS(;
    ram_matrices,
    start_val = start_fabin3,
    gradient = true,
    kwargs...)

    A, S, F, M, parameters =
        ram_matrices.A, ram_matrices.S, ram_matrices.F, ram_matrices.M, ram_matrices.parameters

    n_var, n_nod = size(F)

    A = Matrix(A)
    S = Matrix(S)
    F = Matrix(F); F = convert(Matrix{Float64}, F)

    if check_constants(S)
        @error "constant variance/covariance parameters different from zero are not allowed in SNLLS"
    end

    σ_indices = findall(isone, LowerTriangular(ones(n_var, n_var)))

    # store the indices of the parameters
    A_indices = get_parameter_indices(parameters, A)
    S_indices = get_parameter_indices(parameters, S; index_function = eachindex_lower, linear_indices = true)

    A_indices = [convert(Vector{Int}, indices) for indices in A_indices]
    S_indices = [convert(Vector{Int}, indices) for indices in S_indices]

    A_pars, S_pars = get_partition(A_indices, S_indices)

    # dimension of undirected and directed parameters
    q = size(start_val)
    q_undirected = length(S_pars)
    q_directed = length(A_pars)

    A_indices_linear = A_indices[A_pars]
    A_indices_cartesian = linear2cartesian.(A_indices_linear, [A])

    S_indices = linear2cartesian.(S_indices, [S])
    S_indices = S_indices[S_pars]

    # A matrix
    A_pre = zeros(size(A)...)
    set_constants!(A, A_pre)

    A_rand = copy(A_pre)
    randpar = rand(length(start_val))

    fill_matrix(
        A_rand,
        A_indices_linear,
        randpar)

    acyclic = isone(det(I-A_rand))

    # check if A is lower or upper triangular
    if iszero(A_rand[.!tril(ones(Bool, size(A)...))])
        A_pre = LowerTriangular(A_pre)
    elseif iszero(A_rand[.!tril(ones(Bool, size(A)...))'])
        A_pre = UpperTriangular(A_pre)
    elseif acyclic
        @info "Your model is acyclic, specifying the A Matrix as either Upper or Lower Triangular can have great performance benefits."
    end

    I_A = zeros(n_nod, n_nod)

    size_σ = Int(0.5*(n_var^2+n_var))

    if !isnothing(M)

        if check_constants(M)
            @error "constant mean parameters different from zero are not allowed in SNLLS"
        end

        M_indices = get_parameter_indices(parameters, M)
        M_indices = [convert(Vector{Int}, indices) for indices in M_indices]
        M_pars = get_partition(M_indices)
        M_indices = M_indices[M_pars]
        q_mean = length(M_pars)

        G = zeros(size_σ+n_var, q_undirected+q_mean)

        G_μ = zeros(n_var, q_mean)

        G_μ_indices = CartesianIndices((size_σ .+ (1:n_var), q_undirected .+ (1:q_mean)))
        # TODO: analyze sparsity pattern of G

        if gradient
            ∇G = zeros((size_σ+n_var)*(q_undirected+q_mean), q_directed)
            # TODO: analyze sparsity pattern of ∇G
        else
            ∇G = nothing
        end

    else

        q_mean = nothing
        M_indices = nothing
        G_μ = nothing
        G_μ_indices = nothing

        G = zeros(size_σ, q_undirected)
        # TODO: analyze sparsity pattern of G

        if gradient
            ∇G = zeros(size_σ*q_undirected, q_directed)
            # TODO: analyze sparsity pattern of ∇G
        else
            ∇G = nothing
        end

    end

    return SNLLS(
        q_directed,
        q_undirected,
        size_σ,

        A_indices_linear,
        A_indices_cartesian,
        S_indices,
        σ_indices,

        A_pre,
        I_A,
        G,
        ∇G,

        q_mean,
        M_indices,
        G_μ,
        G_μ_indices,

        start_val
    )
end

############################################################################
### functors
############################################################################

function (imply::SNLLS)(par, F, G, H, model)

    fill_matrix(
        imply.A_pre,
        imply.A_indices_linear,
        par)

    imply.I_A .= inv(I - imply.A_pre)

    fill_G!(
        imply.G,
        imply.q_undirected,
        imply.size_σ,
        imply.S_indices,
        imply.σ_indices,
        imply.I_A)

    if !isnothing(imply.M_indices)

        FI_A = imply.I_A[1:Int(model.observed.n_man), :]

        fill_G_μ!(
            imply.G_μ,
            imply.q_mean,
            imply.M_indices,
            FI_A)

        copyto!(imply.G, imply.G_μ_indices, imply.G_μ, CartesianIndices(imply.G_μ))

    end

    if G
        fill_∇G!(
            imply.∇G,
            imply.q_undirected,
            imply.size_σ,
            imply.q_directed,
            imply.S_indices,
            imply.σ_indices,
            imply.A_indices_cartesian,
            imply.I_A,
            imply.q_mean,
            imply.M_indices,
            model.observed.n_man)
    end

end

############################################################################
### additional functions
############################################################################

function fill_G!(G, q_undirected, size_σ, S_indices, σ_indices, I_A)

    fill!(G, zero(eltype(G)))

    for s in 1:q_undirected
        for ind in S_indices[s]
            l, k = ind[1], ind[2]
            # rows
            for r in 1:size_σ
                i, j = σ_indices[r][1], σ_indices[r][2]
                G[r, s] += I_A[i, l]*I_A[j, k]
                if l != k
                    G[r, s] += I_A[i, k]*I_A[j, l]
                end
            end
        end
    end

end

function fill_G_μ!(G_μ, q_mean, M_indices, FI_A)

    fill!(G_μ, zero(eltype(G_μ)))

    for i in 1:q_mean
        for j in M_indices[i]
            G_μ[:, i] .+= FI_A[:, j]
        end
    end

end

function fill_∇G!(∇G, q_undirected, size_σ, q_directed, S_indices, σ_indices, A_indices, I_A, q_mean::Nothing, M_indices::Nothing, n_var)

    fill!(∇G, zero(eltype(∇G)))

    for s in 1:q_undirected

        for c ∈ S_indices[s]
            l, k = c[1], c[2]

            for r in 1:size_σ
                i, j = σ_indices[r][1], σ_indices[r][2]
                t = (s-1)*size_σ + r

                for m in 1:q_directed
                    u, v = A_indices[m][1][1], A_indices[m][1][2]
                    ∇G[t, m] += I_A[i, u]*I_A[v, l]*I_A[j, k] + I_A[i, l]*I_A[j, u]*I_A[v, k]
                    if l != k
                        ∇G[t, m] += I_A[i, u]*I_A[v, k]*I_A[j, l] + I_A[i, k]*I_A[j, u]*I_A[v, l]
                    end
                end

            end

        end

    end

end

function fill_∇G!(∇G, q_undirected, size_σ, q_directed, S_indices, σ_indices, A_indices, I_A, q_mean, M_indices, n_var)

    fill!(∇G, zero(eltype(∇G)))

    size_G_rows = size_σ + Int(n_var)

    for s in 1:q_undirected

        for c ∈ S_indices[s]
            l, k = c[1], c[2]

            for r in 1:size_σ
                i, j = σ_indices[r][1], σ_indices[r][2]
                t = (s-1)*size_G_rows + r

                for m in 1:q_directed
                    u, v = A_indices[m][1][1], A_indices[m][1][2]
                    ∇G[t, m] += I_A[i, u]*I_A[v, l]*I_A[j, k] + I_A[i, l]*I_A[j, u]*I_A[v, k]
                    if l != k
                        ∇G[t, m] += I_A[i, u]*I_A[v, k]*I_A[j, l] + I_A[i, k]*I_A[j, u]*I_A[v, l]
                    end
                end
            end
        end
    end

    for s in 1:q_mean

        for k ∈ M_indices[s]

            for r in 1:Int(n_var)
                t =  (q_undirected + s - 1)*size_G_rows + size_σ + r

                for m in 1:q_directed
                    u, v = A_indices[m][1][1], A_indices[m][1][2]
                    ∇G[t, m] += I_A[r, u]*I_A[v, k]
                end
            end
        end
    end

end

############################################################################
### Pretty Printing
############################################################################

function Base.show(io::IO, struct_inst::SNLLS)
    print_type_name(io, struct_inst)
    print_field_types(io, struct_inst)
end


text:
@testset "leadtime.jl" begin
    lt = LeadTime(2, 3, 10)
    @test state(lt) == [10,10,10]
    @test state_size(lt) == 3
    dest = []
    InventoryModels.push!(lt, 5, dest)
    @test state(lt) == [10,10,10,5]
    @test only(lt.pull_orders) == (dest => 5)
    InventoryModels.dispatch!(lt)
    @test dest == [10, lt]
    @test state(lt) == [10,10,5]
    @test InventoryModels.reward!(lt) == -25*2
    @test isempty(lt.pull_orders)
    @test InventoryModels.inventory_position(lt) == 25
    reset!(lt)
    @test state(lt) == [10,10,10]
    @test isempty(lt.pull_orders)

    lt = LeadTime(3, Uniform(5,15))
    @test all(5 .< state(lt) .< 15)
    @test state_size(lt) == 3
    next = first(state(lt))
    dest = []
    InventoryModels.push!(lt, 5, dest)
    @test state(lt)[4] == 5
    @test only(lt.pull_orders) == (dest => 5)
    InventoryModels.dispatch!(lt)
    @test dest == [next, lt]
    @test state(lt)[3] == 5
    InventoryModels.reward!(lt) == 0
    @test isempty(lt.pull_orders)
    reset!(lt)
    @test all(5 .< state(lt) .< 15)
    @test isempty(lt.pull_orders)
end


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

DOCUMENTER_DEBUG=true

using Documenter, MultiStochGrad

makedocs(
    format = Documenter.HTML(prettyurls = false),
    sitename = "MultiStochGrad",
    pages = Any[
        "Introduction" => "INTRO.md",
        "MultiStochGrad.jl " => "index.md",
        "scsg.jl" => "scsg.md",
    ]
)


text:
#
#  CellListFunctions.jl
#  CollagenPolymer
#
#  Created by Christopher Revell on 03/08/2020.
#
#

# Functions from https://discourse.julialang.org/t/cell-list-algorithm-is-slower-than-double-for-loop/36621/2

module CellListFunctions

# Import Julia packages
using LinearAlgebra
using Dictionaries
using StaticArrays

# Import local program modules

function find_pairs(nParticles, pos, s, nGrid, neighbourCells)

    pairsList = Tuple{Int64, Int64}[]           # Array of tuples storing neighbour pairs
    boundaryList = Tuple{Int64, Int64, Int64}[]

    # Allocate all particles in matrix pos to grid points
    cellLists = gridAllocate(pos,nParticles,s)

    for key in keys(cellLists)
        # Find neighbour pair list from particles in the same grid cell and neighbouring grid cells
        sameCellPairs!(cellLists[key], pairsList, pos, s)
        loopNeighbour!(key, neighbourCells, cellLists, pairsList, pos, s)

        # Check edge grid cells to find list of particles interacting with boundary
        for jj=1:3
    		if key[jj]==1
    			for kk in cellLists[key]
    				push!(boundaryList,(kk,jj,1))
    			end
    		end
    		if key[jj]==nGrid
    			for kk in cellLists[key]
    				push!(boundaryList,(kk,jj,nGrid))
    			end
    		end
    	end

    end

    return pairsList,boundaryList

end

# Allocate all particles in matrix pos to grid points. Return Dictionary cellLists mapping (x,y,z) indices to list of particles.
function gridAllocate(pos, N, interactionThresh)

    cellLists = Dictionary{SVector{3,Int64},Vector{Int64}}()

    for i = 1:N
        indexSVector = ceil.(Int64,pos[i]/interactionThresh)
        if indexSVector ∉ keys(cellLists)
            insert!(cellLists,indexSVector,[i])
        else
            push!(cellLists[indexSVector],i)
        end
    end

    return cellLists
end

# Find list of neighbouring cells around given cell index
function getForwardCell!(neighbourCells, index)
    neighbourCells[1]  = (index[1]+1, index[2]  , index[3]  )
    neighbourCells[2]  = (index[1]+1, index[2]  , index[3]-1)
    neighbourCells[3]  = (index[1]+1, index[2]  , index[3]+1)
    neighbourCells[4]  = (index[1]+1, index[2]+1, index[3]  )
    neighbourCells[5]  = (index[1]+1, index[2]+1, index[3]-1)
    neighbourCells[6]  = (index[1]+1, index[2]+1, index[3]+1)
    neighbourCells[7]  = (index[1]+1, index[2]-1, index[3]  )
    neighbourCells[8]  = (index[1]+1, index[2]-1, index[3]-1)
    neighbourCells[9]  = (index[1]+1, index[2]-1, index[3]+1)
    neighbourCells[10] = (index[1]  , index[2]  , index[3]+1)
    neighbourCells[11] = (index[1]  , index[2]+1, index[3]  )
    neighbourCells[12] = (index[1]  , index[2]+1, index[3]+1)
    neighbourCells[13] = (index[1]  , index[2]+1, index[3]-1)

    return neighbourCells
end

function sameCellPairs!(cellList, pairsList, pos, interactionThresh)
    len = length(cellList)
    if len >= 2
        for i=1:len-1
            for j=i+1:len
                addOrNot!(cellList[i], cellList[j], pos, interactionThresh, pairsList)
            end
        end
    end
    return pairsList
end

#
function loopNeighbour!(key, neighbourCells, cellLists, pairsList, pos,interactionThresh)
    getForwardCell!(neighbourCells, key)
    for neighbourIndex in neighbourCells
        if neighbourIndex ∈ keys(cellLists)
            findPair!(cellLists[key], cellLists[neighbourIndex], pairsList, pos, interactionThresh)
        end
    end
    return pairsList
end

# Function to loop over all particles in two given grid points and then call function  to test whether each particle pair is in interaction range
function findPair!(cellList1, cellList2, pairsList, pos, interactionThresh)
    for P1 in cellList1
        for P2 in cellList2
            addOrNot!(P1, P2, pos, interactionThresh, pairsList)
        end
    end
    return pairsList
end

# Test whether two points P1 and P2 are within separation limit interactionThresh; if so, add tuple (P1,P2) to vector of pairs pairsList
function addOrNot!(P1, P2, pos, interactionThresh, pairsList)
    dx = pos[P1] - pos[P2]
    dist = norm(dx)
    if dist <= interactionThresh
        push!(pairsList, (P1,P2)) # Add this pair to array of neighbour pairs if within range
    end
    return pairsList
end

export find_pairs

end

max_stars_repo_name: chris-revell/collagen-model-jl

max_issues_repo_name: chris-revell/collagen-model-jl

max_forks_repo_name: chris-revell/collagen-model-jl


text:
push!(LOAD_PATH, "../LF_common/")
using m_LF3d
include("../LF_common/sparse_LF3d.jl")
include("build_nabla2_x.jl")

const ⊗ = kron

function test_main( NN::Array{Int64} )

  Nx = NN[1]
  Ny = NN[2]
  Nz = NN[3]

  AA = [0.0, 0.0, 0.0]
  BB = [16.0, 16.0, 16.0]

  # Initialize LF
  LF = init_LF3d_c( NN, AA, BB, verbose=true )

  @printf("Using loop")
  @time Lx_v1 = build_nabla2_x(LF)
  #println(Lx_v1)

  #@printf("Using Kronecker product directly")
  #@time Lx_v2 = LF.LFx.D2jl ⊗ speye(Ny) ⊗ speye(Nz)

  # Check whether there are any differences between the two matrices
  #println(Lx_v1 - Lx_v2)
  #println(Lx_v2)

end

test_main([101,101,101])


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

"""
    AbstractScalarShape{T} <: AbstractValueShape
"""
abstract type AbstractScalarShape{T} <: AbstractValueShape end

export AbstractScalarShape

@inline Base.size(::AbstractScalarShape) = ()
@inline Base.length(::AbstractScalarShape) = 1

@inline default_unshaped_eltype(shape::AbstractScalarShape{T}) where {T<:Real} =
    default_datatype(T)

@inline default_unshaped_eltype(shape::AbstractScalarShape{Complex}) = default_datatype(Real)

@inline default_unshaped_eltype(shape::AbstractScalarShape{<:Complex{T}}) where {T} =
    default_unshaped_eltype(_valshapeoftype(T))

@inline shaped_type(shape::AbstractScalarShape{<:Real}, ::Type{T}) where {T<:Real} = T

@inline shaped_type(shape::AbstractScalarShape{<:Complex}, ::Type{T}) where {T<:Real} = Complex{T}

"""
    ScalarShape{T} <: AbstractScalarShape{T}

An `ScalarShape` describes the shape of scalar values of a given type.

Constructor:

    ScalarShape{T::Type}()

T may be an abstract type of Union, or a specific type, e.g.

    ScalarShape{Real}()
    ScalarShape{Integer}()
    ScalarShape{Float32}()
    ScalarShape{Complex}()

Scalar shapes may have a total number of degrees of freedom
(see [`totalndof`](@ref)) greater than one, e.g. shapes of complex-valued
scalars:

    totalndof(ScalarShape{Real}()) == 1
    totalndof(ScalarShape{Complex}()) == 2

See also the documentation of [`AbstractValueShape`](@ref).
"""
struct ScalarShape{T} <: AbstractScalarShape{T} end

export ScalarShape

import Base.<=
@inline <=(a::ScalarShape{T}, b::ScalarShape{U}) where {T,U} = T<:U

@inline _valshapeoftype(T::Type{<:Number}) = ScalarShape{T}()

@inline totalndof(::ScalarShape{T}) where {T <: Real} = 1

@inline function totalndof(::ScalarShape{T}) where {T <: Any}
    if @generated
        fieldtypes = ntuple(i -> fieldtype(T, i), Val(fieldcount(T)))
        field_flatlenghts = sum(U -> totalndof(_valshapeoftype(U)), fieldtypes)
        l = prod(field_flatlenghts)
        quote $l end
    else
        fieldtypes = ntuple(i -> fieldtype(T, i), Val(fieldcount(T)))
        field_flatlenghts = sum(U -> totalndof(_valshapeoftype(U)), fieldtypes)
        l = prod(field_flatlenghts)
        l
    end
end

(shape::ScalarShape{T})(::UndefInitializer) where {T<:Number} = zero(default_datatype(T))

function unshaped(x::Union{T,AbstractArray{T,0}}, shape::ScalarShape{U}) where {T<:Real,U<:Real}
    T <: U || throw(ArgumentError("Element type $T of scalar value not compatible with type $U of given scalar shape"))
    unshaped(x)
end

replace_const_shapes(f::Function, shape::ScalarShape) = shape

const ScalarAccessor{T} = ValueAccessor{ScalarShape{T}} where {T}

@inline view_idxs(idxs::AbstractUnitRange{<:Integer}, va::ScalarAccessor{<:Real}) = first(idxs) + va.offset

# ToDo: view_idxs for scalars with dof greater than 1 (complex, etc.)

Base.@propagate_inbounds function _bcasted_view(data::AbstractArrayOfSimilarVectors{<:Real,N}, va::ScalarAccessor) where N
    flat_data = flatview(data)
    idxs = view_idxs(axes(flat_data, 1), va)
    colons = map(_ -> :, Base.tail(axes(flat_data)))
    view(flat_data, idxs, colons...)
end

Base.copy(instance::VSBroadcasted2{N,typeof(getindex), AbstractArrayOfSimilarVectors{<:Real,N},Ref{<:ScalarAccessor}}) where N =
    copy(_bcasted_view(instance.args[1], instance.args[2][]))

Base.copy(instance::VSBroadcasted2{N,typeof(view), AbstractArrayOfSimilarVectors{<:Real,N},Ref{<:ScalarAccessor}}) where N =
    _bcasted_view(instance.args[1], instance.args[2][])


text: include("newton.jl")