text:
include(joinpath(dirname(@__FILE__), "PlainDocs.jl"))
include(joinpath(dirname(@__FILE__), "ExampleAside.jl"))
include(joinpath(dirname(@__FILE__), "NoNewline.jl"))
include(joinpath(dirname(@__FILE__), "Lexicon134.jl"))

Docile.Cache.register_module(joinpath(dirname(@__FILE__), "PlainDocs.jl"))
Docile.Cache.register_module(joinpath(dirname(@__FILE__), "ExampleAside.jl"))
Docile.Cache.register_module(joinpath(dirname(@__FILE__), "NoNewline.jl"))
Docile.Cache.register_module(joinpath(dirname(@__FILE__), "Lexicon134.jl"))

import PlainDocs
import ExampleAside
import NoNewline

facts("PlainDocs.") do

    metadata = Docile.Interface.metadata(PlainDocs)
    entries  = Docile.Interface.entries(metadata)

    context("General.") do

        @fact length(entries) --> 87

        @fact Docile.Interface.metadata(metadata) --> @compat(
            Dict{Symbol, Any}(
                :format  => :md,
                :exports => Set([:PlainDocs]),
                :manual  => String[],
                )
            )

        package_data = Docile.Cache.getpackage(PlainDocs)
        module_data  = Docile.Cache.getmodule(PlainDocs)

        @fact Docile.Cache.getmeta(package_data) --> @compat(
            Dict{Symbol, Any}(
                :format => Docile.Formats.MarkdownFormatter,
                )
            )
        @fact Docile.Cache.getmeta(module_data) --> @compat(
            Dict{Symbol, Any}(
                :exports => Set([:PlainDocs]),
                )
            )

    end

    context("Functions.") do

        @fact rawdocs(entries, fmeth(PlainDocs.f_0)) --> "f_0/0"
        @fact rawdocs(entries, fmeth(PlainDocs.f_1)) --> "f_1/1"
        @fact rawdocs(entries, fmeth(PlainDocs.f_2)) --> "f_2/1"
        @fact rawdocs(entries, fmeth(PlainDocs.f_3)) --> "f_3/1"
        @fact rawdocs(entries, fmeth(PlainDocs.f_4)) --> "f_4/1"

        @fact rawdocs(entries, fmeth(PlainDocs.f_5)) --> "f_5"
        @fact rawdocs(entries, fmeth(PlainDocs.f_6)) --> "f_6"
        @fact rawdocs(entries, fmeth(PlainDocs.f_7)) --> "f_7"
        @fact rawdocs(entries, fmeth(PlainDocs.f_8)) --> "f_8"
        @fact rawdocs(entries, fmeth(PlainDocs.f_9)) --> "f_9"
        @fact rawdocs(entries, fmeth(PlainDocs.f_10)) --> "f_10"
        @fact rawdocs(entries, fmeth(PlainDocs.f_11)) --> "f_11"
        @fact rawdocs(entries, fmeth(PlainDocs.f_12)) --> "f_12"
        @fact rawdocs(entries, fmeth(PlainDocs.f_13)) --> "f_13"
        @fact rawdocs(entries, fmeth(PlainDocs.f_14)) --> "f_14"
        @fact rawdocs(entries, fmeth(PlainDocs.f_15)) --> "f_15"
        @fact rawdocs(entries, fmeth(PlainDocs.f_16)) --> "f_16"
        @fact rawdocs(entries, fmeth(PlainDocs.f_17)) --> "f_17"
        @fact rawdocs(entries, fmeth(PlainDocs.f_18)) --> "f_18"
        @fact rawdocs(entries, fmeth(PlainDocs.f_19)) --> "f_19"

        @fact rawdocs(entries, fmeth(PlainDocs.f_20)) --> "f_20"
        @fact rawdocs(entries, fmeth(PlainDocs.f_21)) --> "f_21"
        @fact rawdocs(entries, fmeth(PlainDocs.f_22)) --> "f_22"
        @fact rawdocs(entries, fmeth(PlainDocs.f_23)) --> "f_23"
        @fact rawdocs(entries, fmeth(PlainDocs.f_24)) --> "f_24"
        @fact rawdocs(entries, fmeth(PlainDocs.f_25)) --> "f_25"
        @fact rawdocs(entries, fmeth(PlainDocs.f_26)) --> "f_26"
        @fact rawdocs(entries, fmeth(PlainDocs.f_27)) --> "f_27"
        @fact rawdocs(entries, fmeth(PlainDocs.f_28)) --> "f_28"
        @fact rawdocs(entries, fmeth(PlainDocs.f_29)) --> "f_29"

        @fact rawdocs(entries, fmeth(PlainDocs.f_30)) --> "f_30"
        @fact rawdocs(entries, fmeth(PlainDocs.f_31)) --> "f_31"
        @fact rawdocs(entries, fmeth(PlainDocs.f_32)) --> "f_32"
        @fact rawdocs(entries, fmeth(PlainDocs.f_33)) --> "f_33"

        @fact rawdocs(entries, PlainDocs.f_1) --> "g_f_1"
        @fact rawdocs(entries, PlainDocs.f_2) --> "g_f_2"

        @fact rawdocs(entries, PlainDocs.f_34) --> "f_34"
        @fact rawdocs(entries, PlainDocs.f_35) --> "f_35"

        for m in methods(PlainDocs.f_36)
            @fact rawdocs(entries, m) --> "f_36"
        end

        @fact rawdocs(entries, PlainDocs.f_37)  --> "f_37"
        @fact rawdocs(entries, PlainDocs.f_37!) --> "f_37"

        @fact rawdocs(entries, fmeth(PlainDocs.f_38)) --> "f_38/f_39"
        @fact rawdocs(entries, fmeth(PlainDocs.f_39)) --> "f_38/f_39"

        @fact rawdocs(entries, PlainDocs.f_40) --> "f_40"
        @fact rawdocs(entries, PlainDocs.f_40!) --> "f_40"

        @fact rawdocs(entries, meth(PlainDocs.f_41, PlainDocs.tup())) --> "f_41/f_42"
        @fact rawdocs(entries, meth(PlainDocs.f_41, PlainDocs.tup(Any))) --> "f_41/f_42"
        @fact rawdocs(entries, meth(PlainDocs.f_42, PlainDocs.tup())) --> "f_41/f_42"
        @fact rawdocs(entries, meth(PlainDocs.f_42, PlainDocs.tup(Any))) --> "f_41/f_42"

        @fact rawdocs(entries, meth(PlainDocs.f_43, PlainDocs.tup(Any, Any, Any))) --> "f_43"

        @fact rawdocs(entries, PlainDocs.f_44) --> "f_44, f_45"
        @fact rawdocs(entries, PlainDocs.f_45) --> "f_44, f_45"

    end

    context("Globals.") do

        @fact rawdocs(entries, :G_M_1) --> "G_M_1"
        @fact rawdocs(entries, :G_M_2) --> "G_M_2"

        @fact rawdocs(entries, :G_C_1) --> "G_C_1"
        @fact rawdocs(entries, :G_C_2) --> "G_C_2"

    end

    context("Macros.") do

        @fact rawdocs(entries, macrofunc(PlainDocs, "m_1")) --> "m_1"
        @fact rawdocs(entries, macrofunc(PlainDocs, "m_2")) --> "m_2"

    end

    context("Types.") do

        @fact rawdocs(entries, PlainDocs.T_A_1) --> "T_A_1"
        @fact rawdocs(entries, PlainDocs.T_A_2) --> "T_A_2"

        @fact rawdocs(entries, :T_TA_1) --> "T_TA_1"
        @fact rawdocs(entries, :T_TA_2) --> "T_TA_2"

        @fact rawdocs(entries, PlainDocs.T_M_1) --> "T_M_1"
        @fact rawdocs(entries, PlainDocs.T_M_2) --> "T_M_2"

        @fact rawdocs(entries, PlainDocs.T_I_1) --> "T_I_1"
        @fact rawdocs(entries, PlainDocs.T_I_2) --> "T_I_2"

        @fact rawdocs(entries, PlainDocs.T_IC_1) --> "T_IC_1"

        @fact rawdocs(entries, meth(PlainDocs.T_IC_1, PlainDocs.tup())) --> "T_IC_1/0"
        @fact rawdocs(entries, meth(PlainDocs.T_IC_1, PlainDocs.tup(Any))) --> "T_IC_1/1"

        @fact rawdocs(entries, PlainDocs.T_IC_2) --> "T_IC_2"

        @fact rawdocs(entries,
                      meth(PlainDocs.T_IC_2{Any},
                           PlainDocs.tup())
                      ) --> "T_IC_2/0"
        @fact rawdocs(entries,
                      meth(PlainDocs.T_IC_2{Any},
                           PlainDocs.tup(Any))
                      ) --> "T_IC_2/1-2"
        @fact rawdocs(entries,
                      meth(PlainDocs.T_IC_2{Any},
                           PlainDocs.tup(Any, Any))
                      ) --> "T_IC_2/1-2"

        @fact rawdocs(entries, PlainDocs.T_IC_3) --> "T_IC_3"

        @fact rawdocs(entries,
                      meth(PlainDocs.T_IC_3{Real},
                           PlainDocs.tup(Integer)
                           )
                      ) --> "T_IC_3/1"
        @fact rawdocs(entries,
                      meth(PlainDocs.T_IC_3{Real},
                           PlainDocs.tup(Real,Real)
                           )
                      ) --> "T_IC_3/2"
        @fact rawdocs(entries,
                      meth(PlainDocs.T_IC_3{Real},
                           PlainDocs.tup(Type{Matrix{Real}},
                                         Vector{Real},
                                         Vararg{Int})
                           )
                      ) --> "T_IC_3/3"

        @fact rawdocs(entries, PlainDocs.BT_1)             --> "BT_1"
        @fact docsmeta(entries, :category, PlainDocs.BT_1) --> :bitstype

        @fact rawdocs(entries, PlainDocs.BT_2)             --> "BT_2"
        @fact docsmeta(entries, :category, PlainDocs.BT_2) --> :bitstype

    end

end

facts("Base") do

    context("Cache Base and it's submodules.") do

        @fact isa(Docile.Cache.getparsed(Base), ObjectIdDict) --> true

    end

end

## 'function <name> end' syntax tests. ##

if VERSION >= v"0.4-dev+4989"
    include(joinpath(dirname(@__FILE__), "FunctionSyntax", "FunctionSyntax.jl"))
    Docile.Cache.register_module(joinpath(dirname(@__FILE__), "FunctionSyntax", "FunctionSyntax.jl"))
    import FunctionSyntax

    facts("Function Syntax.") do
        metadata = Docile.Interface.metadata(FunctionSyntax)
        entries  = Docile.Interface.entries(metadata)

        @fact length(entries) --> 2

        @fact Docile.Interface.metadata(metadata) --> @compat(
            Dict{Symbol, Any}(
                :format  => :md,
                :exports => Set([:FunctionSyntax]),
                :manual  => String[]
                )
            )

        @fact rawdocs(entries, FunctionSyntax.f_1) --> "f_1"
        @fact rawdocs(entries, FunctionSyntax.f_2) --> "f_2\n"

        @fact docsmeta(entries, :category, FunctionSyntax.f_1) --> :function
        @fact docsmeta(entries, :category, FunctionSyntax.f_2) --> :function
    end
end

facts("Misc.") do

    context("Module field") do

        obj = Docile.Cache.objects(ExampleAside)[1]
        @fact obj.mod --> ExampleAside

    end

    context("No Newline.") do

        obj = Docile.Cache.objects(NoNewline)[1]
        @fact Docile.Cache.getraw(NoNewline, obj) --> "f"

    end

    context("Lexicon Issue #134.") do

        objs = Docile.Cache.objects(Lexicon134)
        @fact length(objs) --> 2

    end

end


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

using Markdown
using InteractiveUtils

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

# ╔═╡ c07d1eba-ea4c-4ef3-82cd-7db338e40c34
using LinearAlgebra, PlutoUI, Plots

# ╔═╡ a15c0e18-c71c-4e22-9818-50f26b1b6907
module Solution

	using LinearAlgebra

	export solve_quadratic, quadratic_ls, gradient_descent

	abstract type Tracker end

	struct NoTrack <: Tracker
	end

	notrack = NoTrack()

	struct PathTrack{T<:Any} <: Tracker
		xsteps::Array{T,1}
		PathTrack(x::T) where {T} = new{T}(copy(T[x]))
	end

	trace(::NoTrack, x) = nothing
	trace(tracker::PathTrack, x) = push!(tracker.xsteps, copy(x))

	nsteps(tracker::PathTrack) = length(tracker.xsteps) - 1

	function getx(tracker::PathTrack)
		n = length(tracker.xsteps)
		return [x[1] for x in tracker.xsteps], [x[2] for x in tracker.xsteps]
	end

	"""
		solve_quadratic(p::Real, q::Real[, r::Real=0.0])

	Finds the minimizer of an 1-D quadratic system. If the quadratic term `p` is
	negative, the function raises an error.

	Inputs:
		- `p`, `q`, `r`: the terms of the 1-D quadratic system

	Output:
		- xstar: the minimizer, a number
	"""
	function solve_quadratic(p::Real, q::Real, r::Real=0.0)
		@assert p > 0.0
		return - q / p
	end

	"""
		solve_quadratic(P::AbstractMatrix, q::AbstractVector, r::Real=0)

	Finds the minimizer of an N-D quadratic system.
	P is assumed to be a symmetric positive-definite matrix.

	Inputs:
		- P, q, r: the terms of the nD quadratic system

	Output:
		- xstar: the minimizer, an (n x 1) vector
	"""
	function solve_quadratic(P::AbstractMatrix, q::AbstractVector, r::Real=0.0)
		return - P \ q
	end

	"""
		quadratic_ls(P::AbstractMatrix, q::AbstractVector, Δx::AbstractVector,
										x::AbstractVector)

	Find the exact step size that minimized a quadratic system in
	a given point x for a given search direction Dx

	Inputs:
		- P, q: the terms of the nD quadratic system
		- x: starting point
		- Δx: search direction

	Output:
		- t: optimal step size
	"""
	function quadratic_ls(P::AbstractMatrix, q::AbstractVector, Δx::AbstractVector,
										x::AbstractVector)
		t = solve_quadratic(Δx' * P * Δx, Δx' * P * x + Δx ⋅ q)
		return t
	end

	"""
		gradient_descent(P::AbstractArray, q::AbstractVector,
				x₀::AbstractVector; β::Real=0.0, ϵ::Real=1e-6,
				tracker::Tracker=notrack)

	Computes the minimizes of a quadratic system using gradient descent. Optionally
	provide momentum.

	Inputs:
		- P, q: the terms of the nD quadratic system
		- x₀: starting point
		- ϵ: convergence parameter
		- β: momentum parameter
		- track: bool to indicate whether to store the path

	Outputs:
		- xstar: the found minimum
	"""
	function gradient_descent(P::AbstractArray, q::AbstractVector,
				x₀::AbstractVector; β::Real=0.0, ϵ::Real=1e-6,
				track=false)
		@assert 0 ≤ β < 1
		x = x₀  # initial value
		Δx = zero(x)  # pre-allocate a vector for the gradient
		track && (path = [copy(x)])
		while true
			Δx .= (1.0 - β) * (- P * x .- q) .+ β *  Δx
			if norm(Δx) < ϵ
				break
			end
			# determine stepsize using exact line search
			t = quadratic_ls(P, q, Δx, x)
			# perform step
			x .+= t * Δx
			track && push!(path, copy(x))
		end
		if track
			return x, path
		else
			return x
	end
	end

end

# ╔═╡ 752ced8e-1b80-11ec-298c-050f319f7228
md"""
# Quadratic optimization

*STMO*

**Michiel Stock**

![](https://github.com/MichielStock/STMO/blob/master/chapters/02.Quadratic/Figures/logo.png?raw=true)
"""

# ╔═╡ c7c39058-1035-4b34-8209-cc71837b6b38
md"""
## Motivation

Quadratic systems are essential! Firstly, we can often closely approximate systems close to their maximum by a quadratic system. Studying the minimization of quadratic systems can teach us about the minimization of general convex functions. Quadratic systems are also crucial in their own right! Many statistical models, graph problems, molecular models, etc. can be formulated as quadratic systems:

- least-square minimization problems
- inference using multivariate normal distributions
- molecular modeling using spring-mass systems
- signal recovery

"""

# ╔═╡ 0a65fdd9-5579-4a60-a826-d438b6a49caf
md"""
## Warming up: one-dimensional quadratic systems

In the scalar case, a quadratic function is given by

$$f(x) = \frac{1}{2}px^2+qx +r\,,$$

with $p>0$ (we will shortly see why).

Our optimization problem is given by:

$$\min_x\,\frac{1}{2}px^2+qx +r\,.$$

This can easily be solved by setting the first order derivative equal to zero:

$$\frac{\mathrm{d}f(x)}{\mathrm{d}x} = px + q$$

$$px^\star+q = 0 \Leftrightarrow x^\star=\frac{-q}{p}\,.$$

To show that this is the sole minimizer of $f(x)$, we have to prove that the second-order derivative is positive at this point. This means that at that point, the derivative of the function is increasing: a little to the left the function is increasing, a little to the right and the function is decreasing. We have

$$\left.\frac{\mathrm{d}^2f(x)}{\mathrm{d}x^2}\right|_{x^\star} = p\,,$$

so if $p>0$ then $x^\star$ is the minimizer of $f(x)$.

**Assignment 1**

Complete the code for solving the 1-D quadratic system. Use it to find the minimum of
$$4x^2+8x + 2\,.$$
Check the solution graphically.
"""

# ╔═╡ e2309eb6-449d-497e-84f9-2e496f29f5a5
"""
    solve_quadratic(p::Real, q::Real, r::Real=0.0)

Finds the minimizer of a 1-D quadratic system,
raises an error if there is no minimizer (p<0)

Inputs:
    - p, q, r: the coefficients of the 1D quadratic system

Output:
    - xstar: the minimizer
"""
function solve_quadratic(p::Real, q::Real, r::Real=0.0)
    @assert missing
    return missing
end

# ╔═╡ 0a9c4cc2-375b-4e0e-85ea-bc5992a92d07
let
	p = missing
	q = missing
	r = missing

	xstar = missing
end

# ╔═╡ a7f4f2fe-8a4d-4dcf-af55-d87298b5cace
# make a plot of the quadratic curve and its solution

# ╔═╡ 60e49a4d-a497-4edc-bde1-88161e9edd5b
md"""
## Towards $n$-dimensional quadratic systems

Let us directly move from the one-dimensional case to the $n$-dimensional case. We will use vector notation:

$$\mathbf{x} = \begin{bmatrix}
       x_1 \\ \vdots \\ x_n
     \end{bmatrix} \in \mathbb{R}^n\,.$$

A general $n$-dimensional linear system is given by:

$$f(\mathbf{x}) =\frac{1}{2} \mathbf{x}^\top P \mathbf{x} + \mathbf{q}^\top\mathbf{x} + r\,,$$

with $P$ an $n\times n$ symmetric matrix, $\mathbf{q}$ an $n$-dimensional vector and $r$ a scalar. Below, the contours of a two-dimensional example are shown.
"""

# ╔═╡ 3850436e-d1d3-463c-8dde-790b40582c4e
let
	P = [4 1; 1 2]
	q = [-1; 1]
	r = -10

	# define the function
	f_quadr(x) = 0.5x' * P * x + q' * x + r

	contourf(-10:0.1:10, -10:0.1:10, (x, y)->f_quadr([x, y]), color=:speed)

	xlabel!("x1"); ylabel!("x2")

end

# ╔═╡ f125cad0-ee03-4474-b535-c1f24c5045a2
md"**Question 1**

Why is $P$ symmetric?"

# ╔═╡ ef216f91-4ee1-4b2c-8ebf-d34fb18d3a1a
md"""
So we want to solve the problem:

$$\min_\mathbf{x}\,\frac{1}{2}\mathbf{x}^\top P \mathbf{x} + \mathbf{q}^\top\mathbf{x} + r\,.$$

The concept of a derivative is extended towards higher dimensions using the *gradient* operator:

$$\nabla_\mathbf{x} = \begin{bmatrix}
       \frac{\partial \, }{\partial x_1} \\ \vdots \\ \frac{\partial \, }{\partial x_n}
     \end{bmatrix}\,,$$

so that the gradient of $f(\mathbf{x})$ is given by:

$$\nabla_\mathbf{x} f(\mathbf{x}) = \begin{bmatrix}
       \frac{\partial f(\mathbf{x}) }{\partial x_1} \\ \vdots \\ \frac{\partial f(\mathbf{x}) }{\partial x_n}
     \end{bmatrix}\,.$$

From now on, we will drop the subscript in the gradient when clear from context. For those not familiar to vector calculus, the most useful rules are given below. Here, $a, b$ and $c$ are scalars, $f$ and $g$ are arbitrary differentiable functions and $\mathbf{b}$ is a vector.

| rule | example                  |
| :------------- | :------------- |
| linearity      | $\nabla_\mathbf{x}(a f(\mathbf{x}) +b g(\mathbf{x})) = a\nabla_\mathbf{x} f(\mathbf{x}) +b\nabla_\mathbf{x} g(\mathbf{x})$       |
| product rule | $\nabla_\mathbf{x}(f(\mathbf{x}) g(\mathbf{x})) = g(\mathbf{x})\nabla_\mathbf{x} f(\mathbf{x}) + f(\mathbf{x})\nabla_\mathbf{x} g(\mathbf{x})$|
|chain rule | $\nabla_{\mathbf{x}} f(g(\mathbf{x})) = \frac{\partial f}{\partial g}\mid_{g(\mathbf{x})} \nabla_\mathbf{x} f(\mathbf{x})$|
| quadratic term | $\nabla_\mathbf{x} \left(\frac{1}{2}\mathbf{x}^\top A\mathbf{x}\right)= A\mathbf{x}$|
|linear term| $\nabla_\mathbf{x} (\mathbf{b}^\top\mathbf{x})=\mathbf{b}$|
|constant term |$\nabla_\mathbf{x} c = 0$ |
"""

# ╔═╡ de7ef8bd-04c1-4c86-80ef-4e1f339f8113
md"The gradient of the quadratic function is

$$\nabla f(\mathbf{x})=P\mathbf{x} +\mathbf{q}\,.$$

Setting this to zero gives

$$\mathbf{x}^\star=-P^{-1}\mathbf{q}\,.$$

> *Even though the solution contains the inverse of a matrix, it is seldom a good idea to compute a matrix inverse. Instead, use a solver for the linear system $A\mathbf{x}=\mathbf{b}$ (numerically stable).*

How do we know that $\mathbf{x}^\star$ is the minimizer of the quadratic system? For this we have to extend the concept of a second order derivative to $n$ dimensions. We define the *Hessian* as:

$$\nabla^2 f(\mathbf{x}) = \begin{bmatrix}
\frac{\partial^2 f(\mathbf{x})}{\partial {x_{1}^2}} & \frac{\partial^2 f(\mathbf{x})}{\partial x_1 x_2} & \ldots &  \frac{\partial^2 f(\mathbf{x})}{\partial x_1 x_n}\\
\frac{\partial^2 f(\mathbf{x})}{\partial x_1 x_2} & \frac{\partial^2 f(\mathbf{x})}{\partial {x_2}^2} & \ldots & \vdots \\
\vdots & \vdots & \ddots & \vdots \\
\frac{\partial^2 f(\mathbf{x})}{\partial x_1 x_n} & \frac{\partial^2 f(\mathbf{x})}{\partial x_2 x_n} & \ldots & \frac{\partial^2 f(\mathbf{x})}{\partial x_n^2}
\end{bmatrix}\,.$$

For the quadratic system, this boils down to

$$\nabla^2 f(\mathbf{x}) = P\,.$$

The condition for $\mathbf{x}^\star$ to be the minimizer of $f(\mathbf{x})$ is that the Hessian should be *positive-definite* in that point.

> A symmetric $n\times n$ matrix $A$ is positive-definite (in symbols: $A\succ0$), if for any vector $\mathbf{z}\in\mathbb{R}^n$:
> $$\mathbf{z}^\top A \mathbf{z} > 0\,.$$

A matrix is positive-definite if (and only if) all its eigenvalues as positive. In Julia, we can easily check whether a matrix is positive-define using `isposdef`:"

# ╔═╡ 8494232c-26ca-4a90-bd35-128a5bde892d
A = [5 -1; -1 2]

# ╔═╡ 9c47ce39-f015-4d51-afe7-3c6cb946d1f0
isposdef(A)

# ╔═╡ 5bd99c99-d60b-4998-8381-c516426607f5
B = [0.5 3; 3 1]

# ╔═╡ 14987a6a-a984-4b04-b951-2db4e79bcac5
isposdef(B)

# ╔═╡ 00220013-46c0-4269-a85b-125255c52de1
md"""
However, this is an expensive operation (it is based on attempting a Cholesky decomposition), so one should avoid performing unnecessary checks for positive definiteness.

A point $\mathbf{x}^\star$ at which the gradient vanishes is a minimizer if and only if

$$\nabla^2 f(\mathbf{x})|_{\mathbf{x}^\star} \succ 0\,.$$

So, for the quadratic problem, $x^\star$ is the unique minimizer iff $P\succ 0$. This means that along every direction $\mathbf{v}\in \mathbb{R}^n$ to project $\mathbf{x}$, the problem reduces to a one-dimensional quadratic function with a positive second-order constant:

$$x_v = \mathbf{v}^\top \mathbf{x}\\
f'(x_v) = \frac{1}{2}x_v (\mathbf{v}^\top P \mathbf{v}) x_v + (\mathbf{v}^\top \mathbf{q})x_v + r\,,$$

where $\mathbf{v}^\top P \mathbf{v}>0$ if $P\succ 0$, which in turn implies that $f'(x_v)$ has a minimizer.

If $P\succ 0$, the quadratic system is a *convex* function with a single minimizer. In many problems, $P$ is positive-definite, so there is a well-defined solution. We will develop this further in Chapter 2.
"""

# ╔═╡ 16f6f2e8-90d2-499c-a5c8-e69297fb9c42
md"""
**Assignment 2**

Complete the code for solving the $n$-D quadratic system. Use it to find the minimum of

$$f(\mathbf{x}) = \mathbf{x}^\top\begin{bmatrix}4 & 1 \\ 1 & 2\end{bmatrix}\mathbf{x} + \begin{bmatrix}3 \\ 1\end{bmatrix}^\top\mathbf{x} + 12\,.$$
"""

# ╔═╡ 6ff5822a-1c4e-491b-b215-b3a1698e2048
let
	P = missing
	q = missing
	r = missing

	xstar = missing
end

# ╔═╡ da307ee7-55ce-4978-a2d6-7cf3e42a9eca
md"""

**Question 2**

Consider $L_2$ regularized ridge regression:

$$\min_\mathbf{x}\, \frac{1}{2}(\mathbf{y} - B\mathbf{x})^\top(\mathbf{y} - B\mathbf{x}) + \frac{c}{2}\cdot \mathbf{x}^\top\mathbf{x}\,,$$

with $c>0$. Write this in the standard form of a quadratic system and show that it is convex. Give the expression for the minimizer.

## Time and memory complexity of the exact solution

The exact solution for a convex quadratic system hinges on solving a $n\times n$ linear system. Conventional solvers for linear systems have a time complexity of $\mathcal{O}(n^3)$. Solving such systems is doable for problems of moderate size ($n<10000$), but becomes infeasible for large-scale problems (on a standard computer).
"""

# ╔═╡ 35dce8ae-f61b-4036-baab-ff14796b30be
md"""
Storing an $n\times n$ matrix also has a memory requirement of $\mathcal{O}(n^2)$. When $n$ is too large, this cannot fit in main memory. In the remainder of this chapter, we will consider the case when $P$ is too large to work with, while matrix-vector products $P\mathbf{x}$ *can* be computed. Some examples of when such settings occur:

- **low-rank**: $P=B^\top B$, with $B\in \mathbb{R}^{n\times p}$, with $p\ll n$.
- **sparse**: $P$ is a very sparse matrix.
- **structured**: $P$ has a special structure so that $P\mathbf{x}$ can be computed on the fly, e.g. $P_{ij}=i^2j^3$.
- **block**: $P$ is a sparse block matrix (blocks can be loaded and processed independently).
"""

# ╔═╡ 2e6ab057-669b-4d6b-bff0-4e34056217a4
md"""## Descent methods

Instead of computing the solution of a convex quadratic system in one step, we will use *descent methods*. Here, a minimizing sequence $\mathbf{x}^{(k)},\, k=1,\dots$, where

$$\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} +t^{(k)}\Delta \mathbf{x}^{(k)}\,,$$

with $t^{(k)}\geq 0$ called the *step size* (in machine learning often called *learning rate*) and $\Delta \mathbf{x}^{(k)}$ called the *search direction*. Proper descent methods have that

$$f(\mathbf{x}^{(k+1)}) < f(\mathbf{x}^{(k)})\,,$$

except when $\mathbf{x}^{(k)}$ is optimal. For this property to hold, the search direction should satisfy

$$(\Delta \mathbf{x}^{(k)})^\top \nabla f(\mathbf{x}) < 0\,.$$

![Descent and ascent step.](https://github.com/MichielStock/STMO/blob/master/chapters/02.Quadratic/Figures/descent_step.png?raw=true)

"""

# ╔═╡ 42541537-2fc5-41cf-9479-e2325c55d67d
md"""
### General descent algorithm

Below is the general pseudocode of a general descent method:

> **given** a starting point $\mathbf{x}$
>
> **repeat**
>> 1. Determine descent direction $\Delta \mathbf{x}$
>> 2. *Line search*. Choose $t>0$.
>> 3. *Update*. $\mathbf{x}:=\mathbf{x} + t \Delta \mathbf{x}$.
>
> **until** stopping criterion is reached.
>
> **Output**: $\mathbf{x}$

Usually, the convergence criterion is of the form

$$||\nabla f(\mathbf{x})|| < \nu\,.$$

The step size can be chosen in several ways:

- **exact**: $t=\arg\min_{s\geq 0}\, f(\mathbf{x}+s\Delta \mathbf{x})$.
- **approximate**: choose a $t$ that only approximately minimizes $f(\mathbf{x}+s\Delta \mathbf{x})$.
- **decaying**: choose some decaying series, e.g. $t = \frac{1}{\alpha+k}$.
- **constant**: a constant step size (often works fine in practice).

For quadratic systems we can compute the exact step size, as this amounts to a simple one-dimensional quadratic problem:

$$t=\arg\min_{s\geq 0}\, \frac{1}{2}(\mathbf{x}+s\Delta \mathbf{x})^\top P (\mathbf{x}+s\Delta \mathbf{x}) + (\mathbf{x}+s\Delta \mathbf{x})^\top \mathbf{q} + r$$

$$=\arg\min_{s\geq 0}\, \frac{1}{2}s^2(\Delta\mathbf{x})^\top P \Delta\mathbf{x} + s((\Delta \mathbf{x})^\top P\mathbf{x}+(\Delta \mathbf{x})^\top\mathbf{q}) +\text{constant}$$

This can be solved as:

$$t = \frac{-(\Delta\mathbf{x})^\top P \mathbf{x}-(\Delta\mathbf{x})^\top\mathbf{q}}{(\Delta\mathbf{x})^\top P \Delta\mathbf{x}}$$

**Assignment 3**

Complete the code for the exact line search for quadratic systems.
"""

# ╔═╡ 7931547d-7e74-45a7-ab83-3212d705dca3
"""
    quadratic_ls(P::AbstractMatrix, q::AbstractVector, Δx::AbstractVector,
                                    x::AbstractVector)

Find the exact step size that minimized a quadratic system in
a given point x for a given search direction Dx

Inputs:
    - P, q: the terms of the nD quadratic system
    - x: starting point
    - Δx: search direction

Output:
    - t: optimal step size
"""
function quadratic_ls(P::AbstractMatrix, q::AbstractVector, Δx::AbstractVector,
                                    x::AbstractVector)
    t = missing  # compute t
    return t
end

# ╔═╡ 05236a85-dd10-4c10-ab9a-96cee2390d6d
md"""
## Gradient descent

A natural choice for the search direction is the negative gradient:

$$\Delta \mathbf{x} = - \nabla f(\mathbf{x})\,.$$

Remember, for the quadratic system, the gradient was

$$\nabla f(\mathbf{x})=P\mathbf{x} + \mathbf{q}\,,$$

so

$$\Delta \mathbf{x} = - P\mathbf{x} - \mathbf{q}\,.$$

### Gradient descent algorithm for quadratic systems

> **given** a starting point $\mathbf{x}$
> **repeat**
>> 1. *Descent direction*. $\Delta \mathbf{x} := - P\mathbf{x} - \mathbf{q}$
>> 2. *Line search*. Choose optimal $t>0$.
>> 3. *Update*. $\mathbf{x}:=\mathbf{x} + t \Delta \mathbf{x}$.
> **until** stopping criterion is reached.

**Assignment 4**

Complete the code for the gradient descent algorithm. Solve the previous quadratic system and compare it with the analytic solution. Start at $\mathbf{x}^{(0)}=[0,0]^\top$. How many steps do you need for the algorithm to converge? You can access this information from a `Tracker` structure using the function `length`.
"""

# ╔═╡ 91e61db9-57b8-4576-9960-bd74984f57b0
"""
    gradient_descent(P::AbstractArray, q::AbstractVector,
            x₀::AbstractVector; ϵ::Real=1e-6,
            tracker::Tracker=notrack)

Computes the minimizes of a quadratic system using gradient descent. Optionally
provide momentum.

Inputs:
    - P, q: the terms of the nD quadratic system
    - x₀: starting point
    - ϵ: convergence parameter
    - tracker: object of the type `Tracker` to save the steps

Outputs:
    - xstar: the found minimum
"""
function gradient_descent(P::AbstractArray, q::AbstractVector,
            x₀::AbstractVector; ϵ::Real=1e-4)

    x = x0  # initial value
    n_steps = 0
    while true
        Dx = missing  # compute GD direction
        if missing
            break
        end
        t = missing  # step size
        missing  # perform step
    end
    return x
  end

# ╔═╡ 3244fd32-30a3-4e5a-9884-dfb89ce867c7
# solve the system with gradient descent

# ╔═╡ 4b7fbd19-8b0f-4097-95a4-778c48833dcc
@bind γ Slider(1:0.1:20, default=10)

# ╔═╡ 8c80e876-55de-414d-81d4-292231b281ad
md"""
### Illustration

We illustrate the gradient descent algorithm for the following system:

$$\min_{x_1,x_2}\, \frac{1}{2} (x_1^2 + \gamma x_2^2)\,$$

with $\gamma$ set to $γ.
"""

# ╔═╡ 9359bca7-f702-4b9c-a55d-a2a85bda62b7
md"""
## Convergence analysis

We can study the convergence of the gradient descent algorithm by using eigenvalue decomposition. The matrix $P$ can be written as:

$$P = U\Lambda U^\top\,,$$

with

-  $\Lambda=\text{diag}(\lambda_1,\ldots,\lambda_n)$ a matrix with the eigenvalues on the diagonal (sorted from small to large);
-  $U = [\mathbf{u}_1, \ldots, \mathbf{u}_n]$, a matrix with the corresponding eigenvectors.

Note that because $P\succ 0$, all eigenvalues are real and positive, and all eigenvectors form a real orthonormal basis.

Consider the following linear transformation:

$$\mathbf{z}^{(k)}= U^\top ( \mathbf{x}^{\star}-\mathbf{x}^{(k)})\,,$$

which allows us to rewrite the error in closed-form:

$$f(\mathbf{x}^{(k)}) - f(\mathbf{x}^\star) = \frac{1}{2}\sum_{i=1}^n (1-t\lambda_i)^{2k}\lambda_i[(\mathbf{u}_i)^\top(\mathbf{x}^{(0)}-\mathbf{x}^\star)]^2\,.$$

Here, we see that:

1. The error decomposes in independent terms in the eigenspace.
2. The convergence of each term is determined by the *rate*: $|1-t\lambda_i|$. Convergence occurs as a geometric series.
3. The total number of steps until convergence is determined by either the smallest and largest eigenvalue.
4. Optimal value for fixed step size is $t=\frac{1}{\lambda_1+\lambda_n}$.
"""

# ╔═╡ ec6ef0e2-3589-4c88-90fa-4083739142ce
md"""
Furthermore, it can be shown that if we use an exact line search for the step size, the error $f(\mathbf{x}^{(k)}) - f(\mathbf{x}^\star)\leq \epsilon$ we need fewer than

$$\frac{\log((f(\mathbf{x}^{(0)}) - f(\mathbf{x}^\star))/\epsilon)}{\log(1/c)}\,,$$

with $c=1-\frac{\lambda_1}{\lambda_n}<1$. The quantity $\kappa=\frac{\lambda_n}{\lambda_1}$ is called the *condition number* and largely determines the convergence. We observe:

- The quality of the initial guess $f(\mathbf{x}^{(k)}) - f(\mathbf{x}^\star$) has only a logarithmic impact on the number of steps required.
- Only a few extra steps are needed to decrease $\epsilon$ with one order of magnitude.
- If the condition number is large, then $\log(1/c)\approx 1/\kappa$. Large condition numbers require more steps.

Below, we illustrate this bound for different condition numbers.
"""

# ╔═╡ 31937c09-736a-40b6-babe-2d386ddb430e
md"""
## Gradient descent with momentum

> *While finding the gradient of an objective function is a splendid idea, [descending] the gradient directly may not be.* ~ David J.C. MacKay

Even on simple quadratic problems as discussed here, gradient descent often takes a surprisingly large number of steps to converge. This is because the gradient does not necessarily points in the general direction of the minimum. For convex problems, we are only guaranteed that the gradient points in the half-space of the minimum - a rather weak guarantee! Many improvements in gradient descent have been devised. We will briefly discuss a small modification that can lead to a substantial increase in performance.

"""

# ╔═╡ e206f6f7-c198-438b-99df-7c2942c4cc15
md"""
### Steps with memory

$$\Delta \mathbf{x}^{(k+1)} = \beta \Delta \mathbf{x}^{(k)} - (1-\beta)\nabla f(\mathbf{x}^{(k)})$$

$$\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + t^{(k)}\Delta \mathbf{x}^{(k+1)}\,,$$

with $\beta\in[0,1)$ called the *momentum parameter*.

### Gradient descent algorithm with momentum

> **given** a starting point $\mathbf{x}$,  $\beta\in[0,1)$
> **initialize** $\Delta \mathbf{x}= \mathbf{0}$
> **repeat**
>> 1. *Descent direction*. $\Delta \mathbf{x} := \beta \Delta \mathbf{x}- (1-\beta)(P\mathbf{x}+\mathbf{q})$
>> 2. *Line search*. Choose optimal $t>0$.
>> 3. *Update*. $\mathbf{x}:=\mathbf{x} + t \Delta \mathbf{x}$.
> **until** stopping criterion is reached.
> **Output** $\mathbf{x}$
"""

# ╔═╡ 44c75163-e7b7-4071-8d61-3a609f01ffbd
@bind β Slider(0:0.05:1, show_value=true, default=0.4)

# ╔═╡ 1619f8f3-e638-4ff4-ade0-764840c11215
β

# ╔═╡ 58f001c8-460f-42fd-9f96-c2cc9126873b
md"""
## Illustration of momentum

**Assignment 5**

1. Complete the code for gradient descent with momentum. Use it find the solution for the above system, also starting at $\mathbf{x}=[0,0]^\top$. Set $\beta=0.1$. Do you see an improvement?
2. Compare both algorithms for minimizing the following system:

$$f(\mathbf{x}) = \frac{1}{2}\mathbf{x}^\top\begin{bmatrix}500 & 2 \\ 2 & 1\end{bmatrix}\mathbf{x} + \begin{bmatrix}-40 \\100 \end{bmatrix}^\top\mathbf{x} -5\,,$$

at $\mathbf{x}_0= [0, 0]^\top$. Does momentum increase the speed now?
"""

# ╔═╡ ae161416-fae9-4d75-81a4-dfacb82db055
"""
    gradient_descent_momentum(P::AbstractArray, q::AbstractVector,
            x₀::AbstractVector; β::Real=0.0, ϵ::Real=1e-6,
            tracker::Tracker=notrack)

Computes the minimizes of a quadratic system using gradient descent. Optionally
provide momentum.

Inputs:
    - P, q: the terms of the nD quadratic system
    - x₀: starting point
    - ϵ: convergence parameter
    - β: momentum parameter
    - tracker: object of the type `Tracker` to save the steps

Outputs:
    - xstar: the found minimum
"""
function gradient_descent_momentum(P::AbstractArray, q::AbstractVector,
            x₀::AbstractVector; β::Real=0.0, ϵ::Real=1e-6)
    @assert 0 ≤ β < 0
    x = x₀  # initial value
    Δx = zero(x)  # pre-allocate a vector for the gradient
    while true
        Δx .= missing
        if missing
            break
        end
        # determine stepsize using exact line search
        t = quadratic_ls(P, q, Δx, x)
        # perform step
        x .+= missing
    end
    return x
end

# ╔═╡ 62524ce5-75f1-4ddb-a69c-a219c1c7c50e
md"""
## Conjugated gradient descent

Conjugated gradient descent is an important method for approximately minimizing quadratic systems. It converges much faster than simple gradient descent and does not have a hyperparameter as with momentum updates.

The main idea is that in every step the current search direction and all previous search directions are conjugate with respect to $P$, i.e.

$$(\Delta\mathbf{x}^{(k)})^\top  P\Delta\mathbf{x}^{(l)}=0 \quad\forall  k> l\,.$$

- The search directions can be computed \alert{without} storing the previous directions.
- Solves an $n$-dimensional quadratic system exactly in $n$ steps (provided there are no numerical errors).
- Convergence is

$$|\mathbf{e}^{(k)}|_P \le \left(\frac{\kappa -1}{\kappa + 1} \right)^k|\mathbf{e}^{(0)}|_P\,$$

with $|\mathbf{e}^{(k)}|_P = (\mathbf{x}^\star - \mathbf{x}^{(k)})^\top P (\mathbf{x}^\star - \mathbf{x}^{(k)})$.

The interested reader is referred to '[An Introduction to the Conjugate Gradient Method Without the Agonizing Pain](http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf)' by Jonathan Richard Shewchuk for an in-depth overview.
"""

# ╔═╡ 2198ca9d-7522-47e1-925a-14c7fdbc7bc5
md"""
## Exercise: signal recovery

As a practical example of minimizing quadratic systems, let us consider a simple signal recovery problem. Consider an $n$-dimensional real vector $\mathbf{x}=[x_1,\ldots,x_n]^\top$. Rather than observing this vector directly, we have $m$ noisy measurements at random indices (indices drawn randomly with replacement from $\{1,\ldots,n\}$): $\mathcal{O} = \{(i_j, y_j)\mid j=1,\ldots,m\}$. These measurements are stored in an $m$-dimensional vector $\mathbf{y}$. Can we recover $\mathbf{x}$ from $\mathbf{y}$?
"""

# ╔═╡ 28577dc0-0521-4c07-bf9d-338d4b582bf8
md"""
If $m<n$, then we do not have a single measurement for every element of $\mathbf{x}$. Even if $m>n$, it is likely that some elements of $\mathbf{x}$ are not observed due to chance (for large $n$, if $n=m$, then about 37% of the elements will not be sampled). Recovering $\mathbf{x}$ from $\mathbf{y}$ is an impossible assumption if we do not make some assumptions, this seems an impossible problem in general.

If we assume that the different values of $\mathbf{x}$ are on a line, then we can make a *smoothness* assumption: elements of $\mathbf{x}$ for which the indices are close, likely will have similar values. This idea is expressed in the following minimization problem:

$$\min_\mathbf{x}\, \frac{1}{2}\sum_{(i_j, y_j)\in \mathcal{O}}(y_j-{x}_{i_j})^2 + \frac{C}{2} \mathbf{x}^\top K^{-1}\mathbf{x}\,,$$

with $K^{-1}$ an inverse kernel (or covariance matrix) and $C$ a tuning hyperparameter. The matrix $K^{-1}$ encodes how the different elements of $\mathbf{x}$ are related, constructing such a matrix is a topic in machine learning (see course Predictive Modelling). For our purposes, we have chosen this matrix as such that elements should have values closes to each other. Hence, the minimization problem has two terms:

- a data-fitting term to make sure that the recovered vector $\mathbf{x}$ matches the observations,
- a regularization term to ensure the smoothness of the solution.

The parameter $C$ determines the trade-off between the two terms. See below for the heatmap of $K$.
"""

# ╔═╡ 98e95644-913f-46d9-9deb-d78f7102db3b
md"""
The problem can written purely in matrix notation by using the $(m\times n)$ *bookkeeping matrix* $R$ for which $R_{ij}=1$ if the the $j$-th element of $\mathbf{y}$ corresponds to the $i$-th element of $\mathbf{x}$ and $R_{ij}=0$ otherwise. Hence, the compact matrix form is:

$$\min_\mathbf{x}\, \frac{1}{2}(\mathbf{y}-R\mathbf{x})^\top(\mathbf{y}-R\mathbf{x}) + \frac{C}{2} \mathbf{x}^\top K^{-1}\mathbf{x}\,.$$
"""

# ╔═╡ 7cb64961-d94b-4fe7-913d-50bd094553f2
md"""
**Assignments**

1. Write the minimization problem in the standard form.
2. Use the function `generate_noisy_measurements` to generate $m=100$ noisy measurements (standard deviation is 1, default) of a vector with dimensionality $n=1000$. Use the functions `make_connection_matrix` and `make_bookkeeping` to generate the associated matrices $K^{-1}$ and $R$.
3. Use $C=1$, generate $\mathbf{x}^\star$ using the closed-form solution, using gradient descent and gradient descent with momentum. How many steps do the two descent methods need to converge? Use a vector of zeros as the initial point.
4. Minimize the system for values of $C=1\times 10^{-2}, 1\times 10^{-1}, 1, 10, 100$.  Use for momentum $\beta=0, 0.1, 0.2,\ldots, 0.9$. Make a table of the number of steps needed to reach convergence for the different values of $C$ and $\beta$. Make a plot with the different $\mathbf{x^\star}$ for different values of $C$. Note that the result should be (nearly) the same whether you use momentum or not, only the number of steps will differ.
"""

# ╔═╡ 8e74c0e1-a2de-4a3b-9e9e-41c006ea4c7f

# ╔═╡ 3e1f98e8-fea3-40bf-a9dd-9af36afc7877

# ╔═╡ 742b079b-092c-4907-9029-2aa6916c8d2e

# ╔═╡ 313f361f-a93d-403a-856f-5ffa9c2557c7
md"""
## References

- Boyd, S. and Vandenberghe, L., '*[Convex Optimization](https://web.stanford.edu/~boyd/cvxbook/bv_cvxbook.pdf)*'. Cambridge University Press (2004)
- Goh, G. '[*Why Momentum really works*](https://distill.pub/2017/momentum/)' (2017)
"""

# ╔═╡ e5864656-9465-49aa-b997-2c237d16f15b
md"Here are some utility functions and the solutions."

# ╔═╡ 5beb181b-454c-4314-84e2-ee1d0f220c71
begin
	myblue = "#304da5"
	mygreen = "#2a9d8f"
	myyellow = "#e9c46a"
	myorange = "#f4a261"
	myred = "#e76f51"
	myblack = "#50514F"

	mycolors = [myblue, myred, mygreen, myorange, myyellow]
end;

# ╔═╡ dbe09c4b-7ec2-415d-b1ed-c689dca7e8d1
let
	solvetime = Float64[]
	invtime = Float64[]
	sizes = [2, 5, 10, 50, 100, 500, 1000, 5000, 10000]

	P = randn(2, 2)
	q = randn(2)

	# warmup
	@elapsed(P \ q)
	@elapsed(inv(P) * q)

	for n in sizes
		P = randn(n, n) |> x -> x' * x + I
		q = randn(n)
		push!(solvetime, @elapsed(P \ q))
		push!(invtime, @elapsed(inv(P) * q))
	end

	plot(sizes, solvetime, color=myblue, label="time for solving", lw=2, xscale=:log10,
								yscale=:log10)
	plot!(sizes, invtime, color=myred, label="time for inverting P", lw=2)
	xlabel!("n")
	ylabel!("seconds")
end

# ╔═╡ 07f1caf6-db49-4d9f-8ff1-d6e8e6f19db7
let
	P = [1 0; 0 γ]
	q = [0, 0]
	r = 0
	x₀ = [10.0, 1.0]

	f(x1, x2) = [x1, x2] |> x -> x' * P * x / 2 + q' * x + r

	xstar, path = Solution.gradient_descent(P, q, copy(x₀), track=true)

	contourf(-11:0.1:11, -5:0.1:5, f, color=:speed)
	plot!(first.(path), last.(path) , color=myorange, label="GD", lw=2)
	xlabel!("x1"); ylabel!("x2")
end

# ╔═╡ b3f5a694-dd90-47dd-a2e8-e238f3e8bdc7
let
	λs = [0.1, 1.5, 1.6, 1.8, 2]

	x₀ = ones(5)
	t = 2 / (maximum(λs) + minimum(λs))
	# error xstar = 0, x₀ = 1
	errors = [(1.0 - t * λ)^(2k) * λ * x₀ᵢ^2 for k in 0:1000, (λ, x₀ᵢ) in zip(λs, x₀)]
	cumerrors = cumsum(errors, dims=2)
	plot(cumerrors[:,1], fillrange=0, xscale=:log10, color=myblue, fillcolor=myblue, label="lam 1 = $(λs[1])")
	plot!(cumerrors[:,2], fillrange=cumerrors[:,1], color=myred, fillcolor=myred, label="lam 2 = $(λs[2])")
	plot!(cumerrors[:,3], fillrange=cumerrors[:,2], color=mygreen, fillcolor=mygreen, label="lam 3 = $(λs[3])")
	plot!(cumerrors[:,4], fillrange=cumerrors[:,3], color=myorange, fillcolor=myorange, label="lam 4 = $(λs[4])")
	plot!(cumerrors[:,5], fillrange=cumerrors[:,4], color=myyellow, fillcolor=myyellow, label="lam 5 = $(λs[5])")
	xlabel!("k+1")
	ylabel!("f(x) - f(x*)")
	title!("Error of gradient descent\n(cumulative eigencomponents)")
end

# ╔═╡ 0170ef68-cf32-46a2-ac39-ff971517dfb4
let
	λs = [0.1, 1.5, 1.6, 1.8, 2]

	x₀ = ones(5)
	t = 2 / (maximum(λs) + minimum(λs))
	# error xstar = 0, x₀ = 1
	errors = [(1.0 - t * λ)^(2k) * λ * x₀ᵢ^2 for k in 0:1000, (λ, x₀ᵢ) in zip(λs, x₀)]
	cumerrors = cumsum(errors, dims=2)

	plot(errors, ylims=(1e-15,7), lw=2, xscale=:log10, yscale=:log10,
				 color = [myblue myred myorange mygreen myyellow],
				 label=[" lambda  $i = $(λs[i])" for j in 1:1, i in 1:5])
	xlabel!("k+1")
	ylabel!("f(x) - f(x*)")
	title!("Error of gradient descent\n(individual contributions)")
end

# ╔═╡ 41c04f53-091e-4e7f-a920-5df2a3366db6
let
	κs = [1.1, 2.5, 5, 50, 100]

	bestt(κ) = 2 / (1 + κ)
	errors = [(1 - bestt(κ))^(2k) + (1-bestt(κ) * κ)^(2k) * κ for k in 1:10000, κ in κs]

	plot(errors, color = reshape(mycolors, 1, :),
			labels = [" kappa = $κ" for i in 1:1, κ in κs],
			xscale=:log10, yscale=:log10, ylims=(1e-10, 100), lw=2)

	κ2c(κ) = 1 - 1 / κ
	bounds = [exp(log(1+κ) - k * log(1/κ2c(κ))) for k in 1:10000, κ in κs]

	plot!(bounds, lw=2, ls=:dash, color = reshape(mycolors, 1, :), label="")
	ylabel!("f(x) - f(x*)")
	xlabel!("k + 1")
	title!("Convergence (-) and bound (--) of GD\n for different condition numbers")
end

# ╔═╡ 575f94d7-dadb-4217-8d54-68aa5ee07236
let
	P = [1 0; 0 γ]
	q = [0, 0]
	r = 0
	x₀ = [10.0, 1.0]

	f(x1, x2) = [x1, x2] |> x -> x' * P * x / 2 + q' * x + r

	# without momentum
	xstar, path_gd = Solution.gradient_descent(P, q, copy(x₀), track=true)
	# with momentum
	xstar, path_gdm = Solution.gradient_descent(P, q, copy(x₀); β, track=true)

	contourf(-11:0.1:11, -5:0.1:5, f, color=:speed)
	plot!(first.(path_gd), last.(path_gd) , color=myorange, label="GD", lw=2)
	plot!(first.(path_gdm), last.(path_gdm) , color=myblue, label="GDM", lw=2)
	xlabel!("x1"); ylabel!("x2")
end

# ╔═╡ f03d9959-321d-4bd6-b062-1bae6fb28d64
let
	P = [1 0; 0 γ]
	q = [0, 0]
	r = 0
	x₀ = [10.0, 1.0]

	f(x) = x' * P * x / 2 + q' * x + r

	# without momentum
	xstar, path_gd = Solution.gradient_descent(P, q, copy(x₀), track=true)
	# with momentum
	xstar, path_gdm = Solution.gradient_descent(P, q, copy(x₀); β, track=true)

	plot(f.(path_gd), label="GD", color=myorange, lw=2, yaxis=:log)
	plot!(f.(path_gdm), label="GDM", color=myblue, lw=2)
	xlabel!("iteration + 1")
	ylabel!("error")
end

# ╔═╡ ed98b0d2-2493-4bb7-aceb-eabab7074c9a

"""Compute a quadratic function."""
fquad(x::Vector, P::Matrix, q::Vector, r::Real=0.0) = 0.5x' * P * x + q' * x + r

# ╔═╡ 5f0cc86d-0566-4269-b7a9-9a76c7ee5da0

"""Plot a 1-D quadratic function."""
function plot_quadratic(p::Real, q::Real, r::Real, (xl, xu), stepsize=0.1; kwargs...)
    return plot(x -> p*x^2 + q*x+r, xl:stepsize:xu, xlabel="\$x\$"; kwargs...)
end

# ╔═╡ 5fba2a78-1d80-4796-9dac-04c299670231
"""Plot a 2-D quadratic function."""
function plot_quadratic(P::AbstractMatrix, q::AbstractVector, r::Real,
                            (x1l, x1u), (x2l, x2u), stepsize=0.1; kwargs...)
    fun = (x1, x2) -> [x1,x2] |> x -> 0.5x' * P * x + q' * x + r
    return contour(x1l:stepsize:x1u, x2l:stepsize:x2u, fun, xlabel="\$x_1\$",
                                            ylabel="\$x_2\$"; kwargs...)
end

# ╔═╡ 6590053d-fe91-4401-88e2-03bf29b068ce

# ╔═╡ b165fb28-2671-4b0e-8a46-9a008bb6af46
signalfun(x, n) = 3sin(x * 2 * pi / n) +
            2cos(x * 4 * pi / n) +
            sin(x * 4 * pi / n) + 0.8 * cos(x * 12 * pi / n)

# ╔═╡ c139044f-83f3-4582-b301-e41ded50da4c
"""
    generate_noisy_measurements(m, n; σ=1.0)

Generate noisy measurements according to some function f

Inputs:
    - m : number of observations
    - n : dim of x
    - σ : normally distributed noise (default = 1)

Output:
    - y : vector of noisy measurements
    - ind : vector of indices
"""
function generate_noisy_measurements(m, n; σ=1.0)
    ind = rand(1:n, m)
    y = signalfun.(ind, n) .+ randn(m) * σ
    return y, ind
end

# ╔═╡ 7894fa71-96a6-4f9e-9bd5-e715975b2fb5
begin
	n, m = 1000, 100
	y, ind = generate_noisy_measurements(m, n);
	plot(x -> signalfun(x, n), 1:n, label=" x_i ", color=myblue,
            xlabel=" j ", ylabel="value", lw=2)
	scatter!(ind, y, color=myorange, label=" y_j")
end

# ╔═╡ 5eafacee-4c7b-4cbd-adbb-631223341785
"""
    make_connection_matrix(n, γ=100)

Generates the kernel matrix and the inverse kernel matrix

Uses γ as a characteristic length scale of the radial basis kernel,
assumes periodic boundaries.

Inputs:
    - n : number of points
    - γ : length scale

Output:
    - K
    - Kinv
"""
function make_connection_matrix(n, γ=100)
    # cyclic distance
    d(i, j) = min((i-j)^2, ((i-n)-j)^2, ((j-n)-i)^2)
    # distance matrix
    D = [d(i, j) for i in 1:n, j in 1:n]
    # kernel matrix
    K = exp.(-D / γ^2) + 1e-2I
    return K, inv(K)
end

# ╔═╡ bb203652-f242-4d21-97d6-82419b18daad
K, Kinv = make_connection_matrix(n)

# ╔═╡ d9630c7b-5840-4a87-bde7-c16629665cb2
heatmap(log10.(K), title="Heatmap of K")

# ╔═╡ 1fe7b368-ad79-42b8-a682-94b59cb2f0f6
"""
    make_bookkeeping(I, n)

Constructs the bookkeeping matrix R

Inputs:
    - ind : the indices of the measurements
    - n : dimensionality of signal vector

Output:
    - R : m x n bookkeeping matrix
"""
function make_bookkeeping(ind, n)
    m = length(ind)
    R = zeros(Int, m, n)
    for (i, j) in enumerate(ind)
        R[i,j] = 1
    end
    return R
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
Plots = "~1.22.1"
PlutoUI = "~0.7.10"
"""

# ╔═╡ 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"

[[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 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[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 = "c2178cfbc0a5a552e16d097fae508f2024de61a3"
uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
version = "0.59.0"

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

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

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

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

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

[[HypertextLiteral]]
git-tree-sha1 = "72053798e1be56026b81d4e2682dbe58922e5ec9"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.0"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[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 = "2537ed3c0ed5e03896927187f5f2ee6a4ab342db"
uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
version = "1.0.14"

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

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

[[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 = "7ad0dfa8d03b7bcf8c597f59f5292801730c55b8"
uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c"
version = "0.4.1"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[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 = "1162ce4a6c4b7e31e0e6b14486a6986951c73be9"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.5.2"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[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:
# ╟─752ced8e-1b80-11ec-298c-050f319f7228
# ╠═c07d1eba-ea4c-4ef3-82cd-7db338e40c34
# ╟─c7c39058-1035-4b34-8209-cc71837b6b38
# ╟─0a65fdd9-5579-4a60-a826-d438b6a49caf
# ╠═e2309eb6-449d-497e-84f9-2e496f29f5a5
# ╠═0a9c4cc2-375b-4e0e-85ea-bc5992a92d07
# ╠═a7f4f2fe-8a4d-4dcf-af55-d87298b5cace
# ╟─60e49a4d-a497-4edc-bde1-88161e9edd5b
# ╠═3850436e-d1d3-463c-8dde-790b40582c4e
# ╟─f125cad0-ee03-4474-b535-c1f24c5045a2
# ╟─ef216f91-4ee1-4b2c-8ebf-d34fb18d3a1a
# ╟─de7ef8bd-04c1-4c86-80ef-4e1f339f8113
# ╠═8494232c-26ca-4a90-bd35-128a5bde892d
# ╠═9c47ce39-f015-4d51-afe7-3c6cb946d1f0
# ╠═5bd99c99-d60b-4998-8381-c516426607f5
# ╠═14987a6a-a984-4b04-b951-2db4e79bcac5
# ╟─00220013-46c0-4269-a85b-125255c52de1
# ╟─16f6f2e8-90d2-499c-a5c8-e69297fb9c42
# ╠═6ff5822a-1c4e-491b-b215-b3a1698e2048
# ╟─da307ee7-55ce-4978-a2d6-7cf3e42a9eca
# ╟─dbe09c4b-7ec2-415d-b1ed-c689dca7e8d1
# ╟─35dce8ae-f61b-4036-baab-ff14796b30be
# ╟─2e6ab057-669b-4d6b-bff0-4e34056217a4
# ╟─42541537-2fc5-41cf-9479-e2325c55d67d
# ╠═7931547d-7e74-45a7-ab83-3212d705dca3
# ╟─05236a85-dd10-4c10-ab9a-96cee2390d6d
# ╠═91e61db9-57b8-4576-9960-bd74984f57b0
# ╠═3244fd32-30a3-4e5a-9884-dfb89ce867c7
# ╟─8c80e876-55de-414d-81d4-292231b281ad
# ╟─4b7fbd19-8b0f-4097-95a4-778c48833dcc
# ╟─07f1caf6-db49-4d9f-8ff1-d6e8e6f19db7
# ╟─9359bca7-f702-4b9c-a55d-a2a85bda62b7
# ╟─b3f5a694-dd90-47dd-a2e8-e238f3e8bdc7
# ╟─0170ef68-cf32-46a2-ac39-ff971517dfb4
# ╟─ec6ef0e2-3589-4c88-90fa-4083739142ce
# ╟─41c04f53-091e-4e7f-a920-5df2a3366db6
# ╟─31937c09-736a-40b6-babe-2d386ddb430e
# ╟─e206f6f7-c198-438b-99df-7c2942c4cc15
# ╠═1619f8f3-e638-4ff4-ade0-764840c11215
# ╟─44c75163-e7b7-4071-8d61-3a609f01ffbd
# ╟─575f94d7-dadb-4217-8d54-68aa5ee07236
# ╟─f03d9959-321d-4bd6-b062-1bae6fb28d64
# ╟─58f001c8-460f-42fd-9f96-c2cc9126873b
# ╠═ae161416-fae9-4d75-81a4-dfacb82db055
# ╟─62524ce5-75f1-4ddb-a69c-a219c1c7c50e
# ╟─2198ca9d-7522-47e1-925a-14c7fdbc7bc5
# ╠═7894fa71-96a6-4f9e-9bd5-e715975b2fb5
# ╟─28577dc0-0521-4c07-bf9d-338d4b582bf8
# ╠═bb203652-f242-4d21-97d6-82419b18daad
# ╠═d9630c7b-5840-4a87-bde7-c16629665cb2
# ╠═98e95644-913f-46d9-9deb-d78f7102db3b
# ╟─7cb64961-d94b-4fe7-913d-50bd094553f2
# ╠═8e74c0e1-a2de-4a3b-9e9e-41c006ea4c7f
# ╠═3e1f98e8-fea3-40bf-a9dd-9af36afc7877
# ╠═742b079b-092c-4907-9029-2aa6916c8d2e
# ╟─313f361f-a93d-403a-856f-5ffa9c2557c7
# ╟─e5864656-9465-49aa-b997-2c237d16f15b
# ╠═5beb181b-454c-4314-84e2-ee1d0f220c71
# ╠═ed98b0d2-2493-4bb7-aceb-eabab7074c9a
# ╠═5f0cc86d-0566-4269-b7a9-9a76c7ee5da0
# ╠═5fba2a78-1d80-4796-9dac-04c299670231
# ╠═6590053d-fe91-4401-88e2-03bf29b068ce
# ╠═a15c0e18-c71c-4e22-9818-50f26b1b6907
# ╠═b165fb28-2671-4b0e-8a46-9a008bb6af46
# ╠═c139044f-83f3-4582-b301-e41ded50da4c
# ╠═5eafacee-4c7b-4cbd-adbb-631223341785
# ╠═1fe7b368-ad79-42b8-a682-94b59cb2f0f6
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

ext: jl

max_stars_repo_name: mathematiguy/STMO

max_issues_repo_name: mathematiguy/STMO

max_forks_repo_name: mathematiguy/STMO


text:
include("libnnpack_types.jl")
include("error.jl")
include("libnnpack.jl")
include("performance.jl")
include("interface.jl")

const depsjl_path = joinpath(dirname(@__FILE__), "..", "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("NNPACK not installed properly, run Pkg.build(\"NNlib\"), restart Julia and try again")
end

const shared_threadpool_dict = Dict{UInt64, Base.RefValue}()

"""
    is_nnpack_available()

Checks if the current hardware is supported by NNPACK.
"""
function is_nnpack_available()
    check_deps() isa Nothing || return false
    status = nnp_initialize()
    if status == nnp_status_unsupported_hardware
        return false
    else
        return true
    end
end

"""
    allocate_threadpool()

Allocates several threadpool based on the upper limit on the number of threads for the machine.
Allows NNPACK to intelligently choose which threadpool to use for getting the best
performance.
"""
function allocate_threadpool()
    global NNPACK_CPU_THREADS = NNPACK_CPU_THREADS > 8 ? UInt64(8) : floor(log2(NNPACK_CPU_THREADS))
    for i in 1:Int(NNPACK_CPU_THREADS)
        threads = UInt64(2^i)
        push!(shared_threadpool_dict, threads => Ref(pthreadpool_create(threads)))
    end
end

@init begin
    check_deps()
    status = nnp_initialize()
    if status == nnp_status_unsupported_hardware
        @warn "Hardware is unsupported by NNPACK so falling back to default NNlib"
    end
    try
        global NNPACK_CPU_THREADS = parse(UInt64, ENV["NNPACK_CPU_THREADS"])
    catch
        # Sys.CPU_THREADS should be a better default if we are tuning the benchmark suite on
        # a particular machine. However, we fix the runtime threadpool here to have a max of
        # 4 threads so anything above will be ignored anyways
        global NNPACK_CPU_THREADS = UInt64(4)
    end
    allocate_threadpool()
end


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script to compile data for some descriptive statistics
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

using DataFrames, CSV, XLSX, LinearAlgebra, Statistics, StatsBase, StatsPlots, Dates

# other scripts
dir_home = "x:/VIVES/1-Personal/Florian/git/2022_01_comext/src/"
include(dir_home * "functions.jl")

# location of data input/output (io)
dir_io = "C:/Users/u0148308/data/comext/"
dir_dropbox = "C:/Users/u0148308/Dropbox/BREXIT/"

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# total number of products BE exports to WORLD, UK

function n_products(df::DataFrame, declarants::Vector{String}, partners::Vector{String}, flow::String)

    # subset dataframe
    df = subset(df, :DECLARANT_ISO => ByRow(x -> x in declarants), :FLOW => ByRow(x -> x == flow))
    subset!(df, :PRODUCT_NC => ByRow(x -> x != "TOTAL")) # take out total
    subset!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> !ismissing(x))) # take out missing values
    subset!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> x != 0)) # take out 0 (there are some introduced by adding EU)
    df = df[:, [:DECLARANT_ISO, :PARTNER_ISO, :FLOW, :PRODUCT_NC, :PERIOD, :VALUE_IN_EUROS]]

    # WORLD
    cols_grouping = ["DECLARANT_ISO", "FLOW", "PRODUCT_NC", "PERIOD"]
    gdf = groupby(df, cols_grouping)
    df_WORLD = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)
    row_WORLD = [df_WORLD.DECLARANT_ISO[1] "WORLD" flow df_WORLD.PERIOD[1] length(unique(df_WORLD.PRODUCT_NC)) sum(df_WORLD.VALUE_IN_EUROS)]

    # GB
    df_GB = subset(df, :PARTNER_ISO => ByRow(x -> x in partners))
    cols_grouping = ["DECLARANT_ISO", "FLOW", "PRODUCT_NC", "PERIOD"]
    gdf = groupby(df_GB, cols_grouping)
    df_GB = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)
    row_GB = [df_GB.DECLARANT_ISO[1] "GB" flow df_GB.PERIOD[1] length(unique(df_GB.PRODUCT_NC)) sum(df_GB.VALUE_IN_EUROS)]

    cols_names = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD", "PRODUCT_NC_COUNT", "VALUE_IN_EUROS"]
    output = DataFrame(vcat(row_WORLD, row_GB), cols_names)

    return output
end

#a = n_products(df, ["BE"], ["GB"], "exports")

# timespan
years = string.(2001:2021)
months = lpad.(1:12, 2, '0')

df_fig5 = DataFrame(DECLARANT_ISO=String[], PARTNER_ISO=String[], FLOW=String[], PERIOD=Int64[],
                    PRODUCT_NC_COUNT=Int64[], VALUE_IN_EUROS=Float64[])

for i in years
    for j in months

        #import and clean data
        df = initial_cleaning(i, j)

        # append data for figures/tables
        append!(df_fig5, n_products(df, ["BE"], ["GB"], "exports"))

        println(" ✓ Data for $j/$i has been successfully added. \n")
    end
end

CSV.write(dir_io * "clean/" * "df_fig5" * ".csv", df_fig5)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

path = dir_io * "clean/" * "df_fig5" * ".csv"
df = CSV.read(path, DataFrame)

# formatting
transform!(df, [:DECLARANT_ISO, :PARTNER_ISO, :FLOW] .=> ByRow(string), renamecols=false)
transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)

# average product value
transform!(df, [:VALUE_IN_EUROS, :PRODUCT_NC_COUNT] => ByRow( (x, n) -> x/n) => :AVG_PROD_VALUE)

# HP filter
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
transform!(df, [:AVG_PROD_VALUE] => (x -> HP(x, 20)) => :AVG_PROD_VALUE_HP)

# product count
p = @df df plot(:DATE, :PRODUCT_NC_COUNT,
        group=:PARTNER_ISO, lw=2, legend=:topleft, ylabel="products", title="Belgian export product count: \n World vis-a-vis Great Britain")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
savefig(p, dir_io * "clean/images/fig5/" * "fig5_" * "WORLD_GB" * "_product_count" * ".png") # export image
savefig(p, dir_dropbox * "results/images/fig5/" * "fig5_" * "WORLD_GB" * "_product_count" * ".png") # export image

# average product value
p = @df df plot(:DATE, :AVG_PROD_VALUE,
        group=:PARTNER_ISO, lw=2, legend=:topleft, ylabel="euros", title="Average Belgian export value per product: \n Total vis-a-vis Great Britain")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
savefig(p, dir_io * "clean/images/fig5/" * "fig5_" * "WORLD_GB" * "_avg_product_value" * ".png") # export image
savefig(p, dir_dropbox * "results/images/fig5/" * "fig5_" * "WORLD_GB" * "_avg_product_value" * ".png") # export image

# HP
p = @df df plot(:DATE, :AVG_PROD_VALUE_HP,
        group=:PARTNER_ISO, lw=2, legend=:topleft, ylabel="euros", title="Average Belgian export value per product (HP): \n Total vis-a-vis Great Britain")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
savefig(p, dir_io * "clean/images/fig5/" * "fig5_" * "WORLD_GB" * "_avg_product_value_HP" * ".png") # export image
savefig(p, dir_dropbox * "results/images/fig5/" * "fig5_" * "WORLD_GB" * "_avg_product_value_HP" * ".png") # export image


text:
using HDF5

"""
Create a lock file
"""
lock_create(filename) = open(filename,"w")

const LOCK_SH = 1
const LOCK_EX = 2
const LOCK_NB = 4
const LOCK_UN = 8
has_job     = false
job_info    = nothing
job_results = nothing

flock(fd::Int, op::Int) = ccall(:flock, Int, (Int, Int), fd, op)

lock_read(lock::Int)   = flock(lock, LOCK_SH)
lock_write(lock::Int)  = flock(lock, LOCK_EX)
lock_unlock(lock::Int) = flock(lock, LOCK_UN)

"""
Aquire a job from the task scheduler
"""
function get_job(task_file)
    h5open(task_file, "r+") do file
        id = 1
        while(true)
            #println('*')
            addr   = string("/sched/",id,"/taken")
            script = string("/sched/",id,"/script")
            func   = string("/sched/",id,"/func")
            args   = string("/sched/",id,"/args")
            if(!exists(file, addr))
                #println(":(")
                global has_job
                has_job = false
                break
            elseif(read(file[addr]) == 0)
                #println(":)")
                println("job id = $id")
                global has_job
                global job_info
                has_job  = true
                arg_name = read(file[args])
                arg_val  = map(x->read(file[string("/input/",id,"/",x)]),arg_name)
                job_info = (id, read(file[script]), read(file[func]), arg_val)
                o_delete(file, addr)
                file[addr] = 1
                break
            end
            id += 1
        end
    end

    #println("get_job...")
    #has_job = false
end

"""
Run the job
"""
function perform_job()
    #println("Performing job: ", job_info)
    reload(job_info[2])
    global job_results
    job_results = eval(parse(string(job_info[3], "(job_info[4]...)")))
    #println("job results = ", job_results)
end

"""
Put the results into the task file
"""
function dump_results(task_file)
    h5open(task_file, "r+") do file
        if(!isa(job_results, Tuple))
            file[string("/results/",job_info[1], "/result")] = job_results
        else
            for i=1:length(job_results)
                file[string("/results/",job_info[1], "/result/",i)] = job_results[i]
            end
        end
    end
end

"""
Consume jobs from task_file, using lock_file to avoid concurency issues until no
jobs remain
"""
function child_main(task_file, lock_file)
    #println("Main...");
    global has_job
    lock    = lock_create(lock_file)
    lock_id = fd(lock)
    while(true)
        lock_write(lock_id)
        get_job(task_file)
        lock_unlock(lock_id)
        if(!has_job)
            break
        end

        perform_job()

        lock_write(lock_id)
        dump_results(task_file)
        lock_unlock(lock_id)
    end
end


text:
using HELICS
using Test
const h = HELICS

@testset "Versions" begin

    @test h.helicsGetVersion() isa String
    @test VersionNumber(split(h.helicsGetVersion())[1]) == h.HELICS_LIBRARY_VERSION

end

@testset "ValueFederate" begin
    include("valuefederate.jl")
end

@testset "MessageFederate" begin
    include("messagefederate.jl")
end

@testset "CombinationFederate" begin
    include("combinationfederate.jl")
end

@testset "MessageFilter" begin
    include("messagefilter.jl")
end

@testset "API" begin
    include("api.jl")
end

@testset "Query" begin
    include("query.jl")
end


text:
# Setup plotting environment
using Pkg;
Pkg.add("Plots");
using Plots;
Pkg.add("PyPlot");
Pkg.add("PyCall");
Pkg.add("LaTeXStrings");
pyplot();

# Define a Gaussian shaped function
# Note - we put loop inside function so it can operate on Array object and
# so it can execute efficiently. Calling function for each element of Array
# is more expensive.
ndist(x,μ,σ)=[1/√(2*π*σ^2)*exp(-((xx-μ)^2)/(2*σ^2)) for xx in x]

# Define a list of values and apply Gaussian function to them, showing timing.
x=collect(-10:0.1:10);
μ=0.;σ=π;@time yy=ndist(x,μ,σ)
plot(x,yy)

Pkg.add("Statistics")
using Statistics;
print(Statistics.std(yy))

Pkg.add("InteractiveUtils")
using InteractiveUtils;
@code_native ndist(x,μ,σ);

max_stars_repo_name: christophernhill/julia-intro

max_issues_repo_name: christophernhill/julia-intro

max_forks_repo_name: christophernhill/julia-intro


text:
module test_bukdu_controllers

using Test # @test_throws @test
using Bukdu # routes get ApplicationController Routing redirect_to Router

struct FrontController <: ApplicationController
    conn::Conn
end
index(c::FrontController) = redirect_to(c.conn, "/path")

@test_throws Routing.AbstractControllerError routes(() -> get("/", ApplicationController, index))

routes() do
    get("/", FrontController, index)
end
result = Router.call(get, "/")
@test result.resp.status == 302
@test result.route.action === index

Routing.empty!()

end # module test_bukdu_controllers