text:
######################################################################################################
# Copy from https://bitbucket.org/wmkoolen/tidnabbil/src/master/purex_games_paper/ by Wouter Koolen  #
#  * Baselines are all implemented by Wouter Koolen.                                                 #
######################################################################################################
import Distributions: Multinomial
include("../utilities/regret.jl");
include("../utilities/tracking.jl");
include("../utilities/envelope.jl");
###################################################################### Ours ######################################################################
"""
Ours: Frank-Wolfe based Sampling
"""
struct FWSampling end
long(sr::FWSampling) = "FW-Sampling";
abbrev(sr::FWSampling) = "FWS";
mutable struct FWSamplingState
    x;
    FWSamplingState(K) = new([1.0/K for i=1:K]);
end
function start(sr::FWSampling, N)
    FWSamplingState(length(N));
end
function nextsample(sr::FWSamplingState, pep, istar, ξ, N, S)
    K = length(N);
    t = sum(N);
    r = t^(-9.0/10)/K;
    hw = N./t; # emp. sampled rate
    hμ = S./N; # emp. estimates
    hi = argmax(hμ);
    z = [0.0 for i=1:K];
    if !hμ_in_lambda(hμ, hi, K) || is_complete_square(floor(Int, t/K))
        z = [1.0/K for i=1:K];
    else
        f, ∇f, fidx = compute_f_∇f_standard_bai(sr.x, hμ, ξ, hi, r, K);
        if length(fidx) == 1 # best challenger
            challenger_idx = argmax(∇f[fidx[1]]);
            z = [(challenger_idx==j) ? 1 : 0 for j=1:K];
        else # solve LP of the zero-sum matrix game
            Σ = [[(i==j) ? 1 : 0 for j=1:K]-sr.x for i=1:K];
            A = [[Σ[i]'∇f[j] for i=1:K] for j in fidx]; # construct payoff matrix
            z = solveZeroSumGame(A, K, length(fidx));
        end
    end
    setfield!(sr, :x, sr.x*((t-1.0)/t) + z*1.0/t);
    return argmax(sr.x ./ hw);
end

###################################################################### Baselines ######################################################################
"""
Uniform sampling implemented by Wouter Koolen
"""
struct RoundRobin end
long(sr::RoundRobin) = "Uniform";
abbrev(sr::RoundRobin) = "RR";
function start(sr::RoundRobin, N)
    return sr;
end
function nextsample(sr::RoundRobin, pep, istar, ξ, N, S)
    return 1+(sum(N) % length(N));
end

"""
Track and Stop (Garivier and Kaufmann, 2016) implemented by Wouter Koolen
"""
struct TrackAndStop
    TrackingRule;
end
long(sr::TrackAndStop) = "Track-and-Stop " * abbrev(sr.TrackingRule);
abbrev(sr::TrackAndStop) = "T-" * abbrev(sr.TrackingRule);

struct TrackAndStopState
    t;
    TrackAndStopState(TrackingRule, N) = new(
        ForcedExploration(TrackingRule(N))
    );
end
function start(sr::TrackAndStop, N)
    TrackAndStopState(sr.TrackingRule, N);
end
function nextsample(sr::TrackAndStopState, pep, istar, ξ, N, S)
    K = length(N);
    t = sum(N);
    # oracle at ξ (the closest feasible bandit model)
    _, w = oracle(pep, ξ);
    # tracking
    return track(sr.t, N, w);
end

"""
Optimistic Track and Stop (Degenne, Koolen and Ménard, 2019) implemented by Wouter Koolen
Wouter's implementation has a bug: `optimistic_oracle` fails infrequently (<5 failures in our 3000 repeated experiments).
Our workaround to this is to switch to `oracle` if `optimistic_oracle` fails.
"""
struct OptimisticTrackAndStop
    TrackingRule;
end
long(sr::OptimisticTrackAndStop) = "Optimistic TaS " * abbrev(sr.TrackingRule);
abbrev(sr::OptimisticTrackAndStop) = "O-" * abbrev(sr.TrackingRule);
struct OptimisticTrackAndStopState
    t;
    OptimisticTrackAndStopState(TrackingRule, N) = new(
        TrackingRule(N)
    );
end
function start(sr::OptimisticTrackAndStop, N)
    OptimisticTrackAndStopState(sr.TrackingRule, N);
end
function nextsample(sr::OptimisticTrackAndStopState, pep, istar, ξ, N, S)
    K = length(N);
    t = sum(N);
    # Optimistic oracle at ξ (the closest feasible bandit model)
    w = minimum(ξ);
    _, w = optimistic_oracle(pep, ξ, N);
    try
        _, w = optimistic_oracle(pep, ξ, N);
    catch e
        _, w = oracle(pep, ξ);
    end
    # tracking
    return track(sr.t, N, w);
end

"""
Gradient Ascent algorithm (Ménard, 2019) implemented by Wouter Koolen.
"""
struct Menard
    TrackingRule;
    scale;
end
long(sr::Menard) = "Menard " * abbrev(sr.TrackingRule);
abbrev(sr::Menard) = "M-" * abbrev(sr.TrackingRule);

function start(sr::Menard, N)
    MenardState(sr.TrackingRule, sr.scale, N);
end
struct MenardState
    h;
    t;
    MenardState(TrackingRule, scale, N) = new(
        FixedShare(length(N), S=scale),
        TrackingRule(N)
    );
end
function nextsample(sr::MenardState, pep, istar, ξ, N, S)
    K = length(N);
    t = sum(N);
    hμ = S./N; # emp. estimates
    # query the learner
    w = act(sr.h);
    # best response λ-player
    _, (_, λs), (_, _) = glrt(pep, w, hμ);
    # gradient
    ∇ = [d(getexpfam(pep, k), hμ[k], λs[k]) for k in eachindex(hμ)];
    # update learner
    incur!(sr.h, -∇);
    # tracking
    return track(sr.t, N, w);
end

"""
AdaHedge vs Best-Response (Section 3.1 in Degenne, Koolen and Ménard, 2019) implemented by Wouter Koolen
"""
struct DaBomb
    TrackingRule;
    M;
end
long(sr::DaBomb) = "DaBomb " * abbrev(sr.TrackingRule);
abbrev(sr::DaBomb) = "D-" * abbrev(sr.TrackingRule);
struct DaBombState
    hs; # one online learner per answer
    t;
    DaBombState(TrackingRule, N, M) = new(
        map(x -> AdaHedge(length(N)), 1:M),
        TrackingRule(N)
    );
end
function start(sr::DaBomb, N)
    DaBombState(sr.TrackingRule, N, sr.M);
end
# optimistic gradients
function optimistic_gradient(pep, hμ, t, N, λs)
    [let dist = getexpfam(pep, k),
     ↑ = dup(dist, hμ[k], log(t)/N[k]),
     ↓ = ddn(dist, hμ[k], log(t)/N[k])
     max(d(dist, ↑, λs[k]), d(dist, ↓, λs[k]), log(t)/N[k])
     end
     for k in eachindex(hμ)];
end
function nextsample(sr::DaBombState, pep, istar, ξ, N, S)
    K = length(N);
    t = sum(N);
    hμ = S./N; # emp. estimates
    # query the learner
    w = act(sr.hs[istar]);
    # best response λ-player
    _, (_, λs), (_, ξs) = glrt(pep, w, hμ);
    ∇ = optimistic_gradient(pep, hμ, t, N, λs);
    incur!(sr.hs[istar], -∇);
    # tracking
    return track(sr.t, N, w);
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 9c340f8b-877d-41b9-baca-2698cc2b5988
using Plots

# ╔═╡ 99adb0e2-1b94-11ec-1351-79707f6016af
md"""
# Bracket search

*STMO*

**Michiel Stock**
"""

# ╔═╡ 6b30fe20-f321-4cb8-a2fc-0a2820d71bc0
md"As a first introduction to mathematical optimization, we will study two *bracketing methods*. Bracketing methods can be used to minimize scalar function with a single input variable. These algorithms identify an interval $[a, b]$ containing the desired minimum. We will assume that the functions are *unimodal*, i.e. they contain a single minimum."

# ╔═╡ bdf6c059-db72-4d90-b98f-2e77b7d89f37
md"""
## Finding an initial bracket

The following function will generate from an initial $x$ an interval containing a local minimum. It works by moving one of the limits of the interval until the function value at that limit starts to increase.
"""

# ╔═╡ 662b4d13-505c-4202-b104-0db80d8b1822
function bracket_minimum(f, x=0.0; s=1e-2, k=2.0)
  a, ya = x, f(x)  # lower limit of the interval
  b, yb = x, f(x + s)  #  upper limit of the interval
  if yb > ya  # default assumes f is decreasing in a, flip if not the case
    a, b = b, a
    ya, yb = yb, ya
    s = -s  # go to the left instead of right
  end
  while true
    c, yc = b + s, f(b + s)
    if yc > yb
      return a < c ? (a, c) : (c, a)
    end
    a, ya = b, yb
    b, yb = c, yc
    s *= k
  end
end

# ╔═╡ 42131c7f-2980-4bc9-be00-64e25ecd238e
md"Example of finding an initial bracket, starting from point $x$."

# ╔═╡ 143801cb-25ee-4a49-882a-a8ab377a25f8
md"""
## Bisection method

The bisection method is technically not a minimization algorithm, but a *root finding method*, i.e., it can be used to find an $x$ for a given function $g: \mathbb{R}\rightarrow \mathbb{R}$ such that

$$g(x) = 0\,.$$

We can use the bisection method to find a minimum of a function $f(x)$ by finding points where the first derivative $f'(x)$ is equal to 0.

> **Question**: in addition to the first derivative being equal to zero, which other criteria should hold for $x^\star$ to be a minimizer?

The bisection method departs from an initial bracket $[a, b]$, where $g(a)$ and $g(b)$ have opposing signs. If $g(x)$ is a continuous function, the *intermediate value theorem* states that there is at least one $x\in[a,b]$ such that $g(x)$.

In every step of the bisection method, the interval is cut into half. The midpoint $x = (a+b)/2$ is computed, and a new bracket is formed from the midpoint and the side that contains zero.

"""

# ╔═╡ 0eac24ee-344d-4ac6-ad3d-d3c29e98da09
md"""
The bisection procedure is repeated until the length of the interval is smaller than a small $\epsilon > 0$. This $\epsilon$ is called a *convergence parameter*, because it determines when the algorithm will halt. The pseudocode of the bisection method is given below.

> **given** $g(x)$, the derivative of a function, initial interval $[a, b]$, tolerance $\epsilon$.
>
> **while** $b - a > \epsilon$
>> 1. *Determine midpoint*. $x:=(a+b)/2$
>> 2. *Update*.
>>> **if $g(x)=0$**: $a:=x$, $b:=x$
>>>
>>> **else if $\text{sign}(g(x)) = \text{sign}(g(a))$**: $a:=x$
>>>
>>> **else**: $b:=x$
>
> **Output**: $[a, b]$

Since every step shrinks the interval with a factor two, it is easy to show that this procedure stops within

$$\log_2\left(\frac{b-a}{\epsilon}\right)$$

iterations.
"""

# ╔═╡ 524e2c11-e527-44ec-9510-92bb5905ad7a
md"""
**Assignment 1** Complete the code for the bisection method. Use it to find the minimum of

$$f(x) = \log(e^{8x-3}+2e^{0.2x^2-x+1})$$

Check your results graphically.
"""

# ╔═╡ ff4c2d0b-3b10-41c6-8579-c86710a54bbe
function bisection(g, a, b; ϵ=1e-3)
  @assert a < b "a should be smaller than b"

  # a is minizer?
  g(a) == 0 && return a, a
  # b is minizer?
  g(b) == 0 && return b, b

  while missing
    missing
  end
  return a, b
end

# ╔═╡ 727f53c3-2036-4474-8225-1e9f32c0343c
f(x) = log(exp(8x-3)+2exp(0.2x^2-x+1))

# ╔═╡ 5eecd08d-0567-440b-97a3-121882800d79
bracket_minimum(f, k=1.2)

# ╔═╡ 1a2ff096-fcd5-4b5f-901b-4282627202ac
md"> To type `f′` just type `f\prime<tab>`. Whenever you are unsure how to type a unicode character you can call it using the help, e.g. `?ϵ`."

# ╔═╡ 90e525e6-ac62-4299-ac42-cda344074051
f′(x) = missing  # use Wolfram alpha for now

# ╔═╡ 62d84643-d26f-475b-8431-e8f4b82f77a5

# ╔═╡ fc9563d5-2de3-4e11-9836-d41c1c33073a

# ╔═╡ 17bd8b62-a5df-42ab-9e5c-0a6d1fd6ad6f
md"""
## Quadratic fit search

As we will see in the next lecture, finding the minimum of quadratic functions of the form

$$q(x) = p_1 + p_2 x + p_3 x^2$$

is trivial. The *quadratic fit search method* will approximate a function by a quadratic function and computes the minimizer of that function.

Given three bracketing points $a<b<c$ and their respective evaluations $y_a=f(a), y_b=f(b)$ and $y_c=f(c)$, we can fit a quadratic curve by solving the following system of equations for the coefficients:

$$\begin{bmatrix}y_a\\ y_b \\ y_c\end{bmatrix} = \begin{bmatrix}1 & a & a^2 \\ 1 & b & b^2 \\ 1 & c & c^2 \end{bmatrix} \begin{bmatrix}p_1\\ p_2 \\ p_3\end{bmatrix}\,.$$

For example, in Julia, we have:
"""

# ╔═╡ 21a3085c-609d-47de-8a96-7303187e5aa4
a, b, c = -2, -1, 3

# ╔═╡ 7beb8b85-4e35-4e13-896c-3fffc992ea13
ya, yb, yc = f.((a, b, c))

# ╔═╡ bff05d34-060f-4860-9c70-70c6cfea6721
p1, p2, p3 = [1 a a^2; 1 b b^2; 1 c c^2] \ [ya, yb, yc]

# ╔═╡ 7f10c92e-e778-43c7-aaa4-38e6150773d7
md"""We can even dispense solving the system and use a closed-form fit based on three points:

$$q(x) = y_a\frac{(x-b)(x-c)}{(a-b)(a-c)} +  y_b\frac{(x-a)(x-c)}{(b-a)(b-c)}+ y_c\frac{(x-a)(x-b)}{(c-a)(c-b)}\,.$$

The unique minimum of this quadratic curve is computed by:

$$x^\star = \frac{1}{2} \frac{y_a(b^2-c^2) + y_b(c^2-a^2)+ y_c(a^2-b^2)}{y_a(b-c) + y_b(c-a)+ y_c(a-b)}\,.$$

We have implemented these two formulas in `quadratic_fit` and `quadratic_fit_min`. See below for an illustration.
"""

# ╔═╡ 55934af9-5dd1-43c2-bf24-9430df6cdbd1
function quadratic_fit(a, b, c, ya, yb, yc)
  @assert a != b != c
  return x -> ya * (x-b) * (x-c) / ((a-b) * (a-c)) +
              yb * (x-a) * (x-c) / ((b-a) * (b-c)) +
              yc * (x-a) * (x-b) / ((c-a) * (c-b))
end

# ╔═╡ d513a33a-c878-4a18-a5e6-309cc7f2937d
function quadratic_fit_min(a, b, c, ya, yb, yc)
  return 0.5 * (ya * (b^2 - c^2) + yb * (c^2 - a^2) + yc * (a^2 - b^2)) /
        (ya * (b-c) + yb * (c-a) + yc * (a-b))
end

# ╔═╡ 693872b2-9223-4782-b9fe-724bcf8bf9ba
x = quadratic_fit_min(a, b, c, ya, yb, yc)

# ╔═╡ 5e5b17b2-b1a4-4c72-b588-7c03ae6788c7
md"""
Given the minimizer $x^\star$ of the quadratic approximation, we can update the interval closer to this value. The quadratic fit search repeats this approach $n$ times. The pseudocode is given below.

> **given** $f(x)$, the function to be minimized, three increasing values $a, b, c$ with $f(a) > f(b)$ and $f(c) > f(b)$, the number of steps $n$.
>
> **repeat** $n$ times
>> 1. *Fit quadratic*
>> 2. *Determine minimizer of $q(x)$*
>> 3. *Update bracket*:
>>> **if $x\in [a,b]$**: $a, b, c := a, x, b$
>>> **else**: $a, b, c := b, x, c$
>
> **Output**: $a, b, c$
"""

# ╔═╡ dfb42ca1-6dbb-4668-90d7-4e5cd1365c99
md"**Assignment 2**

Complete the code for the quadratic fit search method. Use it again to find the minimizer of the provided $f(x)$.
"

# ╔═╡ ed1b64ee-f2a8-4584-9f39-9906306453f4
function quadratic_fit_search(f, a, b, c, n)
  @assert a < b < c "a, b, c not consecutive "
  @assert f(a) > f(b) && f(c) > f(b) "f(b) should be less than f(a) and f(c)"

  for i in 1:n
    missing  # complete me
  end
  return a, b, c
end

# ╔═╡ b0d790f7-3593-4965-aa3b-f2c728d8c155

# ╔═╡ 46f2c4a8-dcde-4fa8-8b08-cc66b1d3b4c7

# ╔═╡ 26036622-b1e7-480c-9a60-fde35c9d6a8b
md"""The coming chapters will discuss the minimization of quadratic functions in much greater detail.

# Exercises

1. Find the maximum of $\sqrt[\leftroot{-2}\uproot{2}x]{x}$.
2. Given a chemical system for which a chemical with concentration is degraded by enzymic conversion and chemical decay. The rate of removal is given by

$$r(x) = 0.02 + 1.2\frac{x}{10+x}\,.$$

The system receives a constant input with a concentration of 0.3. What is the steady-state concentration?"""

# ╔═╡ 412c4bc9-87c3-4e1b-9968-b2f8ebb05d3a

# ╔═╡ b106794f-c997-499a-8bc6-cc2117105f06

# ╔═╡ a093d875-235f-4606-b6ff-85cf05a0f779

# ╔═╡ 57847e93-2316-4d1d-912d-646b0d8b8d24
md"""
## References

- Kochenderfer, M. J. and Wheeler, T., '*Algorithms for Optimization*'. MIT Press (2019)
"""

# ╔═╡ c486a4da-d50c-49f5-98a0-46b417469b06
begin
	myblue = "#304da5"
	mygreen = "#2a9d8f"
	myyellow = "#e9c46a"
	myorange = "#f4a261"
	myred = "#e76f51"
	myblack = "#50514F"

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

# ╔═╡ a3259f42-3ab2-4cca-9179-ecf9c019aa4e
let
	f(x) = 0.003x^4 + 8x^3 - 3x - 8
	(a, b) = bracket_minimum(f, k=1.2)
	plot(f, -0.5:0.01:0.5, xlabel="x", label="f(x)", color=myblue, lw=2)
	vline!([a, b], label="", color=myorange, lw=2)
	vline!([0], label="", color=mygreen, lw=2)
	xticks!([a, b, 0], ["a", "b", "x"])
	yticks!(:none)
end

# ╔═╡ 33b799e2-70a1-4ec2-a0c7-16fb02bce10c
let
	g(x) = tanh(sin(x))

	a, b = -0.75, 1.2

	plot(g, -pi/2:0.1:pi/2, label="g(x)", lw=2, color=myblue, xlabel="x", legend=:left)
	scatter!([a], [g(a)], color=myorange, label="a")
	scatter!([b], [g(b)], color=mygreen, label="b")
	vline!([(b+a)/2], color=myred, label="x = (b-a)/2", lw=2)
	hline!([0], label="", color=myblack, ls=:dash)
	title!("Illustration of the bisection method")
end

# ╔═╡ 98c07443-db0d-4c86-aef4-ceeafd5a95ef
let
	plot(f, -3:0.1:5, label="f(x)", lw=2, xlabel="x", color=myblue)
	scatter!([a, b, c], [ya, yb, yc], label="a, b, c", color=myorange)
	plot!(quadratic_fit(a, b, c, ya, yb, yc), -3:0.1:5, label=" q(x) ", lw=2, color=mygreen)
	vline!([x], label="x", lw=2, color=myred)
end

# ╔═╡ 7a5b2671-44b1-4382-a9a6-3ed86ecf4fae
module Solution
	export bisection, quadratic_fit_search

	function bisection(g, a, b; ϵ=1e-3)
	  @assert a < b "a should be smaller than b"

	  # a is minizer?
	  g(a) == 0 && return a, a
	  # b is minizer?
	  g(b) == 0 && return b, b

	  ga,gb = g(a),g(b)
	  while (b-a)/2 > ϵ
			x = (a+b)/2
			gx = g(x)
			gx == 0 && return x,x
			if sign(gx) == sign(ga)
				a,ga = x,gx
			else
				b,gb = x,gx
			end
	  end
	  return a, b
	end

	function quadratic_fit_search(f, a, b, c, n)
	  @assert a < b < c "a, b, c not consecutive "
	  @assert f(a) > f(b) && f(c) > f(b) "f(b) should be less than f(a) and f(c)"

	  for i in 1:n
		x = quadratic_fit_min(a,b,c,f(a),f(b),f(c))
			if a < x < b
				a,b,c = a,x,b
				else b < x < c
				a,b,c = b,x,c
			end
	  end
	  return a, b, c
	end
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

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

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

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

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

[[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:
# ╟─99adb0e2-1b94-11ec-1351-79707f6016af
# ╠═9c340f8b-877d-41b9-baca-2698cc2b5988
# ╟─6b30fe20-f321-4cb8-a2fc-0a2820d71bc0
# ╟─bdf6c059-db72-4d90-b98f-2e77b7d89f37
# ╠═662b4d13-505c-4202-b104-0db80d8b1822
# ╟─42131c7f-2980-4bc9-be00-64e25ecd238e
# ╟─5eecd08d-0567-440b-97a3-121882800d79
# ╠═a3259f42-3ab2-4cca-9179-ecf9c019aa4e
# ╟─143801cb-25ee-4a49-882a-a8ab377a25f8
# ╟─33b799e2-70a1-4ec2-a0c7-16fb02bce10c
# ╟─0eac24ee-344d-4ac6-ad3d-d3c29e98da09
# ╟─524e2c11-e527-44ec-9510-92bb5905ad7a
# ╠═ff4c2d0b-3b10-41c6-8579-c86710a54bbe
# ╠═727f53c3-2036-4474-8225-1e9f32c0343c
# ╟─1a2ff096-fcd5-4b5f-901b-4282627202ac
# ╠═90e525e6-ac62-4299-ac42-cda344074051
# ╠═62d84643-d26f-475b-8431-e8f4b82f77a5
# ╠═fc9563d5-2de3-4e11-9836-d41c1c33073a
# ╟─17bd8b62-a5df-42ab-9e5c-0a6d1fd6ad6f
# ╠═21a3085c-609d-47de-8a96-7303187e5aa4
# ╠═7beb8b85-4e35-4e13-896c-3fffc992ea13
# ╠═bff05d34-060f-4860-9c70-70c6cfea6721
# ╟─7f10c92e-e778-43c7-aaa4-38e6150773d7
# ╠═55934af9-5dd1-43c2-bf24-9430df6cdbd1
# ╠═d513a33a-c878-4a18-a5e6-309cc7f2937d
# ╠═693872b2-9223-4782-b9fe-724bcf8bf9ba
# ╟─98c07443-db0d-4c86-aef4-ceeafd5a95ef
# ╟─5e5b17b2-b1a4-4c72-b588-7c03ae6788c7
# ╟─dfb42ca1-6dbb-4668-90d7-4e5cd1365c99
# ╠═ed1b64ee-f2a8-4584-9f39-9906306453f4
# ╠═b0d790f7-3593-4965-aa3b-f2c728d8c155
# ╠═46f2c4a8-dcde-4fa8-8b08-cc66b1d3b4c7
# ╟─26036622-b1e7-480c-9a60-fde35c9d6a8b
# ╠═412c4bc9-87c3-4e1b-9968-b2f8ebb05d3a
# ╠═b106794f-c997-499a-8bc6-cc2117105f06
# ╠═a093d875-235f-4606-b6ff-85cf05a0f779
# ╟─57847e93-2316-4d1d-912d-646b0d8b8d24
# ╟─c486a4da-d50c-49f5-98a0-46b417469b06
# ╟─7a5b2671-44b1-4382-a9a6-3ed86ecf4fae
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

ext: jl

max_stars_repo_name: MichielStock/STMO

max_issues_repo_name: MichielStock/STMO

max_forks_repo_name: MichielStock/STMO


text:
using CompScienceMeshes
using BEAST
using LinearAlgebra
using FastBEAST
using IterativeSolvers

CM = CompScienceMeshes

function test_beast_laplace_singlelayer(h; threading=:single,
    farquaddata=nothing, svdrecompress=false)

    Γ = CM.meshsphere(1, h)

    X = lagrangecxd0(Γ)

    𝒱 = Helmholtz3D.singlelayer(wavenumber=0.0)

    hmat = hassemble(𝒱,X,X, treeoptions=BoxTreeOptions(nmin=50), threading=threading, svdrecompress=svdrecompress)

    mat = assemble(𝒱,X,X)
    return mat, hmat
end

mat, hmat_single = test_beast_laplace_singlelayer(0.1)

@test nnz(hmat_single) == 3916760

@test compressionrate(hmat_single) > 0.3
@test estimate_reldifference(hmat_single, mat) ≈ 0 atol=1e-4

mat, hmat_multi = test_beast_laplace_singlelayer(0.1, threading=:multi)

@test nnz(hmat_multi) == 3916760

@test compressionrate(hmat_multi) > 0.3
@test estimate_reldifference(hmat_multi, mat) ≈ 0 atol=1e-4
@test compressionrate(hmat) > 0.3

mat, hmat_single = test_beast_laplace_singlelayer(0.1, farquaddata=quaddata)

@test nnz(hmat_single) == 3916760

@test compressionrate(hmat_single) > 0.3
@test estimate_reldifference(hmat_single, mat) ≈ 0 atol=1e-3

mat, hmat_multi = test_beast_laplace_singlelayer(0.1, threading=:multi, farquaddata=quaddata)

@test nnz(hmat_multi) == 3916760

@test compressionrate(hmat_multi) > 0.3
@test estimate_reldifference(hmat_multi, mat) ≈ 0 atol=1e-3
@test compressionrate(hmat) > 0.3

mat, hmat_svdmulti = test_beast_laplace_singlelayer(0.1, threading=:multi, farquaddata=quaddata, svdrecompress=true)

@test nnz(hmat_svdmulti) == 3356811

@test compressionrate(hmat_svdmulti) > 0.3
@test estimate_reldifference(hmat_svdmulti, mat) ≈ 0 atol=1e-3
@test compressionrate(hmat) > 0.3

#function test_beast_laplace_singlelayer_manufactured(h)

# ##

# mutable struct Monopole{T,P} <: BEAST.Functional
#     location::P
#     wavenumber::T
#     amplitude::T
# end

# function Monopole(p,k,a = 1)
#     T = promote_type(eltype(p), typeof(k), typeof(a))
#     P = similar_type(typeof(p), T)
#     Monopole{T,P}(p,k,a)
# end

# monopolehh3d(;
#     location    = error("missing arguement `location`"),
#     wavenumber   = error("missing arguement `wavenumber`"),
#     ) = Monopole(location, wavenumber)

# function (f::Monopole)(r)
#     p = f.location
#     k = f.wavenumber
#     a = f.amplitude
#     return a * 1/norm(r-p)
# end

#     h=0.2
#     Γ = CM.meshsphere(1, h)

#     X = lagrangecxd0(Γ)

#     κ = 0.0
#     𝒱 = Helmholtz3D.singlelayer(wavenumber=κ)
#     q = Monopole(SVector(0.35,0.35,0.25), 0.0) #monopolehh3d(location = SVector(0.35,0.35,0.25), wavenumber = 0.0)
#     f = assemble(strace(q,Γ), X)
#     hmat = hassemble(𝒱,X,X, nmin=50)
#     mat = assemble(𝒱,X,X)

# ##
#     hx, hch = cg(hmat,f, log=true,reltol=1e-4)

#     x, ch = cg(mat,f, log=true,reltol=1e-4)
# @test norm(hx - x)/norm(x) ≈ 0 atol=1e-3
# ##
#     Φ, Θ = [0.0], range(0,stop=π,length=100)
#     pts = [point(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)) for ϕ in Φ for θ in Θ]
#     ffd = potential(𝒱, pts, x, X)

# #end

# ##
# h=0.02
# Γ = CM.meshsphere(1, h)

# X = lagrangecxd0(Γ)

# κ = 0.0
# 𝒱 = Helmholtz3D.singlelayer(wavenumber=κ)
# hmat = hassemble(𝒱,X,X, nmin=100)
# println(compressionrate(hmat))

max_stars_repo_name: sbadrian/FastBEAST

max_issues_repo_name: sbadrian/FastBEAST

max_forks_repo_name: sbadrian/FastBEAST


text:
abstract type AbstractReshaper end

"""
    StaticReshaper()

Reshaper that constructs [`SizedAbstractArray`](@ref).
"""
struct StaticReshaper <: AbstractReshaper end

"""
    make_reshape(reshaper::AbstractReshaper, ::Type{Size}, data) where Size

Reshape array `data` to size `Size` using method provided by `reshaper`.
"""
function make_reshape(reshaper::AbstractReshaper, ::Type{Size}, data) where {Size}
    return error("make_reshape is not defined for reshaper of type $(typeof(reshaper)), size $(Size) and data of type $(typeof(data)).")
end
function make_reshape(::StaticReshaper, ::Type{Size}, data) where {Size}
    return SizedAbstractArray{Size}(data)
end

"""
    ArrayReshaper()

Reshaper that constructs `Base.ReshapedArray`.
"""
struct ArrayReshaper <: AbstractReshaper end

function make_reshape(::ArrayReshaper, ::Type{Size}, data) where {Size}
    return reshape(data, size_to_tuple(Size))
end

"""
    ShapeSpecification(reshapers, manifolds::Manifold...)

A structure for specifying array size and offset information for linear
storage of points and tangent vectors on the product manifold of `manifolds`.

The first argument, `reshapers`, indicates how a view representing a point
in the [`ProductArray`](@ref) will be reshaped. It can either be an object
of type `AbstractReshaper` that will be applied to all views or a tuple
of such objects that will be applied to subsequent manifolds.

Two main reshaping methods are provided by types [`StaticReshaper`](@ref)
that is faster for manifolds represented by small arrays (up to about 100
elements) and [`ArrayReshaper`](@ref) that is faster for larger arrays.

For example, consider the shape specification for the product of
a sphere and group of rotations:

```julia-repl
julia> M1 = Sphere(2)
Sphere{2}()

julia> M2 = Manifolds.Rotations(2)
Manifolds.Rotations{2}()

julia> reshaper = Manifolds.StaticReshaper()
Manifolds.StaticReshaper()

julia> shape = Manifolds.ShapeSpecification(reshaper, M1, M2)
Manifolds.ShapeSpecification{(1:3, 4:7),Tuple{Tuple{3},Tuple{2,2}},
  Tuple{Manifolds.StaticReshaper,Manifolds.StaticReshaper}}(
  (Manifolds.StaticReshaper(), Manifolds.StaticReshaper()))
```

`TRanges` contains ranges in the linear storage that correspond to a specific
manifold. `Sphere(2)` needs three numbers and is first, so it is allocated the
first three elements of the linear storage (`1:3`). `Rotations(2)` needs four
numbers and is second, so the next four numbers are allocated to it (`4:7`).
`TSizes` describe how the linear storage must be reshaped to correctly
represent points. In this case, `Sphere(2)` expects a three-element vector, so
the corresponding size is `Tuple{3}`. On the other hand, `Rotations(2)`
expects two-by-two matrices, so its size specification is `Tuple{2,2}`.
"""
struct ShapeSpecification{TRanges,TSizes,TReshapers}
    reshapers::TReshapers
end

function ShapeSpecification(reshapers, manifolds::Manifold...)
    sizes = map(m -> representation_size(m), manifolds)
    lengths = map(prod, sizes)
    ranges = UnitRange{Int64}[]
    k = 1
    for len in lengths
        push!(ranges, k:(k + len - 1))
        k += len
    end
    TRanges = tuple(ranges...)
    TSizes = Tuple{map(s -> Tuple{s...}, sizes)...}
    if isa(reshapers, AbstractReshaper)
        rtuple = map(m -> reshapers, manifolds)
        return ShapeSpecification{TRanges,TSizes,typeof(rtuple)}(rtuple)
    end
    return ShapeSpecification{TRanges,TSizes,typeof(reshapers)}(reshapers)
end

"""
    ProductArray(shape::ShapeSpecification, data)

An array-based representation for points and tangent vectors on the
product manifold. `data` contains underlying representation of points
arranged according to `TRanges` and `TSizes` from `shape`.
Internal views for each specific sub-point are created and stored in `parts`.
"""
struct ProductArray{
    TM<:ShapeSpecification,
    T,
    N,
    TData<:AbstractArray{T,N},
    TV<:Tuple,
    TReshaper,
} <: AbstractArray{T,N}
    data::TData
    parts::TV
    reshapers::TReshaper
end

# The two-argument version of this constructor is substantially faster than
# the generic one.
function ProductArray(
    M::Type{ShapeSpecification{TRanges,Tuple{Size1,Size2},TReshapers}},
    data::TData,
    reshapers,
) where {TRanges,Size1,Size2,TReshapers,T,N,TData<:AbstractArray{T,N}}
    views = (
        make_reshape(reshapers[1], Size1, view(data, TRanges[1])),
        make_reshape(reshapers[2], Size2, view(data, TRanges[2])),
    )
    return ProductArray{M,T,N,TData,typeof(views),typeof(reshapers)}(data, views, reshapers)
end
function ProductArray(
    M::Type{ShapeSpecification{TRanges,Tuple{Size1,Size2,Size3},TReshapers}},
    data::TData,
    reshapers,
) where {TRanges,Size1,Size2,Size3,TReshapers,T,N,TData<:AbstractArray{T,N}}
    views = (
        make_reshape(reshapers[1], Size1, view(data, TRanges[1])),
        make_reshape(reshapers[2], Size2, view(data, TRanges[2])),
        make_reshape(reshapers[3], Size3, view(data, TRanges[3])),
    )
    return ProductArray{M,T,N,TData,typeof(views),typeof(reshapers)}(data, views, reshapers)
end
function ProductArray(
    M::Type{ShapeSpecification{TRanges,TSizes,TReshapers}},
    data::TData,
    reshapers,
) where {TRanges,TSizes,TReshapers,T,N,TData<:AbstractArray{T,N}}
    views = map(
        (size, range, reshaper) -> make_reshape(reshaper, size, view(data, range)),
        size_to_tuple(TSizes),
        TRanges,
        reshapers,
    )
    return ProductArray{M,T,N,TData,typeof(views),typeof(reshapers)}(data, views, reshapers)
end
ProductArray(M::ShapeSpecification, data) = ProductArray(typeof(M), data, M.reshapers)

@doc raw"""
    prod_point(M::ShapeSpecification, pts...)

Construct a product point from product manifold `M` based on point `pts`
represented by [`ProductArray`](@ref).

# Example
To construct a point on the product manifold $S^2 × ℝ^2$
from points on the sphere and in the euclidean space represented by,
respectively, `[1.0, 0.0, 0.0]` and `[-3.0, 2.0]` you need to construct shape
specification first. It describes how linear storage of `ProductArray`
corresponds to array representations expected by `Sphere(2)` and `Euclidean(2)`.

    M1 = Sphere(2)
    M2 = Euclidean(2)
    reshaper = Manifolds.StaticReshaper()
    Mshape = Manifolds.ShapeSpecification(reshaper, M1, M2)

Next, the desired point on the product manifold can be obtained by calling
`Manifolds.prod_point(Mshape, [1.0, 0.0, 0.0], [-3.0, 2.0])`.
"""
function prod_point(M::ShapeSpecification, pts...)
    data = mapreduce(vcat, pts) do pt
        return reshape(pt, :)
    end
    # Array(data) is used to ensure that the data is mutable
    # `mapreduce` can return `SArray` for some arguments
    return ProductArray(M, Array(data))
end

@doc raw"""
    submanifold_component(M::Manifold, p, i::Integer)
    submanifold_component(M::Manifold, p, ::Val(i)) where {i}
    submanifold_component(p, i::Integer)
    submanifold_component(p, ::Val(i)) where {i}

Project the product array `p` on `M` to its `i`th component. A new array is returned.
"""
submanifold_component(::Any...)
@inline function submanifold_component(M::Manifold, p, i::Integer)
    return submanifold_component(M, p, Val(i))
end
@inline submanifold_component(M::Manifold, p, i::Val) = submanifold_component(p, i)
@inline submanifold_component(p, ::Val{I}) where {I} = p.parts[I]
@inline submanifold_component(p, i::Integer) = submanifold_component(p, Val(i))

@doc raw"""
    submanifold_components(M::Manifold, p)
    submanifold_components(p)

Get the projected components of `p` on the submanifolds of `M`. The components are returned in a Tuple.
"""
submanifold_components(::Any...)
@inline submanifold_components(M::Manifold, p) = submanifold_components(p)
@inline submanifold_components(p) = p.parts

function Base.BroadcastStyle(
    ::Type{<:ProductArray{ShapeSpec}},
) where {ShapeSpec<:ShapeSpecification}
    return Broadcast.ArrayStyle{ProductArray{ShapeSpec}}()
end

function Base.similar(
    bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{ProductArray{ShapeSpec}}},
    ::Type{ElType},
) where {ShapeSpec,ElType}
    A = find_pv(bc)
    return ProductArray(ShapeSpec, similar(A.data, ElType), A.reshapers)
end

Base.dataids(x::ProductArray) = Base.dataids(x.data)

"""
    find_pv(x...)

`A = find_pv(x...)` returns the first `ProductArray` among the arguments.
"""
@inline find_pv(bc::Base.Broadcast.Broadcasted) = find_pv(bc.args)
@inline find_pv(args::Tuple) = find_pv(find_pv(args[1]), Base.tail(args))
@inline find_pv(x) = x
@inline find_pv(a::ProductArray, rest) = a
@inline find_pv(::Any, rest) = find_pv(rest)

Base.size(x::ProductArray) = size(x.data)

Base.@propagate_inbounds Base.getindex(x::ProductArray, i) = getindex(x.data, i)

Base.@propagate_inbounds Base.setindex!(x::ProductArray, val, i) = setindex!(x.data, val, i)

function Base.:+(
    v1::ProductArray{ShapeSpec},
    v2::ProductArray{ShapeSpec},
) where {ShapeSpec<:ShapeSpecification}
    return ProductArray(ShapeSpec, v1.data + v2.data, v1.reshapers)
end

function Base.:-(
    v1::ProductArray{ShapeSpec},
    v2::ProductArray{ShapeSpec},
) where {ShapeSpec<:ShapeSpecification}
    return ProductArray(ShapeSpec, v1.data - v2.data, v1.reshapers)
end
function Base.:-(v::ProductArray{ShapeSpec}) where {ShapeSpec<:ShapeSpecification}
    return ProductArray(ShapeSpec, -v.data, v.reshapers)
end

function Base.:*(
    a::Number,
    v::ProductArray{ShapeSpec},
) where {ShapeSpec<:ShapeSpecification}
    return ProductArray(ShapeSpec, a * v.data, v.reshapers)
end

number_eltype(a::ProductArray) = number_eltype(a.data)
number_eltype(::Type{TPA}) where {TM,TData,TPA<:ProductArray{TM,TData}} = eltype(TData)

function _show_component(io::IO, v; pre = "", head = "")
    sx = sprint(show, "text/plain", v, context = io, sizehint = 0)
    sx = replace(sx, '\n' => "\n$(pre)")
    return print(io, head, pre, sx)
end

function _show_component_range(io::IO, vs, range; pre = "", sym = "Component ")
    for i in range
        _show_component(io, vs[i]; pre = pre, head = "\n$(sym)$(i) =\n")
    end
    return nothing
end

function _show_product_repr(io::IO, x; name = "Product representation", nmax = 4)
    n = length(x.parts)
    print(io, "$(name) with $(n) submanifold component$(n == 1 ? "" : "s"):")
    half_nmax = div(nmax, 2)
    pre = "  "
    sym = " Component "
    if n ≤ nmax
        _show_component_range(io, x.parts, 1:n; pre = pre, sym = sym)
    else
        _show_component_range(io, x.parts, 1:half_nmax; pre = pre, sym = sym)
        print(io, "\n ⋮")
        _show_component_range(io, x.parts, (n - half_nmax + 1):n; pre = pre, sym = sym)
    end
    return nothing
end

function Base.show(io::IO, ::MIME"text/plain", x::ProductArray)
    return _show_product_repr(io, x; name = "ProductArray")
end

function allocate(x::ProductArray{ShapeSpec}) where {ShapeSpec<:ShapeSpecification}
    return ProductArray(ShapeSpec, allocate(x.data), x.reshapers)
end
function allocate(
    x::ProductArray{ShapeSpec},
    ::Type{T},
) where {ShapeSpec<:ShapeSpecification,T}
    return ProductArray(ShapeSpec, allocate(x.data, T), x.reshapers)
end
allocate(p::ProductArray, ::Type{T}, s::Size{S}) where {S,T} = Vector{T}(undef, S)
allocate(p::ProductArray, ::Type{T}, s::Integer) where {T} = Vector{T}(undef, s)

"""
    ProductRepr(parts)

A more general but slower representation of points and tangent vectors on
a product manifold.

# Example:

A product point on a product manifold `Sphere(2) × Euclidean(2)` might be
created as

    ProductRepr([1.0, 0.0, 0.0], [2.0, 3.0])

where `[1.0, 0.0, 0.0]` is the part corresponding to the sphere factor
and `[2.0, 3.0]` is the part corresponding to the euclidean manifold.
"""
struct ProductRepr{TM<:Tuple}
    parts::TM
end

ProductRepr(points...) = ProductRepr{typeof(points)}(points)

function number_eltype(x::ProductRepr)
    return typeof(reduce(+, one(number_eltype(eti)) for eti in x.parts))
end

allocate(x::ProductRepr) = ProductRepr(map(allocate, submanifold_components(x))...)
function allocate(x::ProductRepr, ::Type{T}) where {T}
    return ProductRepr(map(t -> allocate(t, T), submanifold_components(x))...)
end
allocate(p::ProductRepr, ::Type{T}, s::Size{S}) where {S,T} = Vector{T}(undef, S)
allocate(p::ProductRepr, ::Type{T}, s::Integer) where {S,T} = Vector{T}(undef, s)

function Base.copyto!(x::ProductRepr, y::ProductRepr)
    map(copyto!, submanifold_components(x), submanifold_components(y))
    return x
end

function Base.:+(v1::ProductRepr, v2::ProductRepr)
    return ProductRepr(map(+, submanifold_components(v1), submanifold_components(v2))...)
end

function Base.:-(v1::ProductRepr, v2::ProductRepr)
    return ProductRepr(map(-, submanifold_components(v1), submanifold_components(v2))...)
end
Base.:-(v::ProductRepr) = ProductRepr(map(-, submanifold_components(v)))

Base.:*(a::Number, v::ProductRepr) = ProductRepr(map(t -> a * t, submanifold_components(v)))
Base.:*(v::ProductRepr, a::Number) = ProductRepr(map(t -> t * a, submanifold_components(v)))

Base.:/(v::ProductRepr, a::Number) = ProductRepr(map(t -> t / a, submanifold_components(v)))

function Base.convert(::Type{ProductRepr{TPR}}, x::ProductRepr) where {TPR<:Tuple}
    if @isdefined TPR
        return ProductRepr(convert(TPR, submanifold_components(x)))
    else
        return x
    end
end

function Base.show(io::IO, ::MIME"text/plain", x::ProductRepr)
    return _show_product_repr(io, x; name = "ProductRepr")
end

ManifoldsBase._get_vector_cache_broadcast(::ProductRepr) = Val(false)


text:
using Zygote

export PrioritizedDQNLearner

"""
    PrioritizedDQNLearner(;kwargs...)

See paper: [Prioritized Experience Replay](https://arxiv.org/abs/1511.05952)

# Keywords

- `approximator`::[`AbstractQApproximator`](@ref): used to get Q-values of a state.
- `target_approximator`::[`AbstractQApproximator`](@ref): similar to `approximator`, but used to estimate the target (the next state).
- `loss_fun`: the loss function.
- `γ::Float32=0.99f0`: discount rate.
- `batch_size::Int=32`
- `update_horizon::Int=1`: length of update ('n' in n-step update).
- `min_replay_history::Int=32`: number of transitions that should be experienced before updating the `approximator`.
- `update_freq::Int=4`: the frequency of updating the `approximator`.
- `target_update_freq::Int=100`: the frequency of syncing `target_approximator`.
- `stack_size::Union{Int, Nothing}=4`: use the recent `stack_size` frames to form a stacked state.
- `default_priority::Float64=100.`: the default priority for newly added transitions.
"""
Base.@kwdef mutable struct PrioritizedDQNLearner{Tq<:AbstractQApproximator, Tt<:AbstractQApproximator,Tf, Ts<:Union{Int, Nothing}} <: AbstractLearner
    approximator::Tq
    target_approximator::Tt
    loss_fun::Tf
    stack_size::Ts = 4
    γ::Float32 = 0.99f0
    batch_size::Int = 32
    update_horizon::Int = 1
    min_replay_history::Int = 32
    update_freq::Int = 1
    target_update_freq::Int = 100
    update_step::Int = 0
    default_priority::Float64 = 100.0

    function PrioritizedDQNLearner(approximator::Tq, target_approximator::Tt, loss_fun::Tf, stack_size::Ts, args...) where {Tq, Tt, Tf, Ts}
        copyto!(approximator, target_approximator)  # force sync
        new{Tq,Tt, Tf, Ts}(approximator, target_approximator, loss_fun, stack_size, args...)
    end
end

function update!(learner::PrioritizedDQNLearner{<:NeuralNetworkQ}, batch)
    learner.update_step += 1
    learner.update_step % learner.update_freq == 0 || return nothing

    Q, Qₜ, γ, loss_fun, update_horizon, batch_size = learner.approximator,
        learner.target_approximator,
        learner.γ,
        learner.loss_fun,
        learner.update_horizon,
        learner.batch_size
    states, rewards, terminals, next_states = map(x->to_device(Q, x), (batch.states, batch.rewards, batch.terminals, batch.next_states))
    actions = CartesianIndex.(batch.actions, 1:batch_size)
    priorities = Vector{Float32}()

    gs = gradient(Q) do
        q = batch_estimate(Q, states)[actions]
        q′ = dropdims(maximum(batch_estimate(Qₜ, next_states); dims = 1), dims = 1)
        G = rewards .+ γ^update_horizon .* (1 .- terminals) .* q′

        batch_losses = loss_fun(G, q)
        priorities = (Zygote.dropgrad(batch_losses) .+ 1f-10)
        mean(batch_losses)
    end

    update!(Q, gs)

    if learner.update_step % learner.target_update_freq == 0
        copyto!(Qₜ, Q)
    end

    priorities |> to_host
end

function extract_transitions(
    buffer::CircularTurnBuffer{PRTSA},
    learner::PrioritizedDQNLearner,
)
    if length(buffer) > learner.min_replay_history
        inds = sample_indices(buffer, learner.batch_size, learner.update_horizon, learner.stack_size)
        inds, extract_SARTS(buffer, inds, learner.γ, learner.update_horizon, learner.stack_size)
    else
        nothing
    end
end


text:
using NLS_Solver
using Documenter

DocMeta.setdocmeta!(NLS_Solver, :DocTestSetup, :(using NLS_Solver); recursive=true)

makedocs(;
         modules=[NLS_Solver],
         authors="vincent-picaud <picaud.vincent@gmail.com> and contributors",
         repo="https://github.com/vincent-picaud/NLS_Solver.jl/blob/{commit}{path}#L{line}",
         sitename="NLS_Solver.jl",
         format=Documenter.HTML(;
                                prettyurls=get(ENV, "CI", "false") == "true",
                                canonical="https://vincent-picaud.github.io/NLS_Solver.jl",
                                assets=String[],
                                ),
    pages=[
        "index.md",
        "Getting Started" => [
            "unconstrained_nls.md",
            "bound_constrained_nls.md",
            "nonlinear_regressions.md",
        ],
        "api.md",
    ],
)

deploydocs(;
           repo="github.com/vincent-picaud/NLS_Solver.jl.git",
           devbranch="main"
)


text:
@inbounds for i = ix:ix-1+n
    x[i] = $op( a, y[i] )
end


text:
export sinkhorn, ot_dist, empirical_sinkhorn, dtw, emd, ot_plot1D, empirical_emd
# roadmap: implement all related functions from https://github.com/rflamary/POT/blob/master/ot/bregman.py

@doc raw"""
    sinkhorn(a::Union{PyObject, Array{Float64}}, b::Union{PyObject, Array{Float64}}, M::Union{PyObject, Array{Float64}};
    reg::Float64 = 1.0, iter::Int64 = 1000, tol::Float64 = 1e-9, method::String="sinkhorn")

Computes the optimal transport with Sinkhorn algorithm. The mathematical formulation is

```math
\begin{aligned}
\arg\min_P &\ \left(P, M\right) + \lambda \Omega(\Gamma)\\
\text{s.t.} &\ \Gamma 1 = a\\
&\ \Gamma^T 1 = b\\
& \Gamma \geq 0
\end{aligned}
```
Here $\Omega$ is the entropic regularization. Note if $\lambda$ is very small, the algorithm may encounter numerical instabilities.

The implementation are adapted from https://github.com/rflamary/POT.
"""
function sinkhorn(a::Union{PyObject, Array{Float64}}, b::Union{PyObject, Array{Float64}}, M::Union{PyObject, Array{Float64, 2}};
    reg::Union{PyObject,Float64} = 1.0, iter::Int64 = 1000, tol::Float64 = 1e-9, method::String="sinkhorn", returnall::Bool=false)
    if isa(a, Array)
        @assert sum(a)≈1.0
        @assert all(a .>= 0)
    end
    if isa(b, Array)
        @assert sum(b)≈1.0
        @assert all(b .>= 0)
    end
    a = convert_to_tensor(a)
    b = convert_to_tensor(b)
    M = convert_to_tensor(M)
    reg = convert_to_tensor(reg)
    iter = convert_to_tensor(iter)
    tol = convert_to_tensor(tol)
    if method=="lp"
        pth = install("OTNetwork")
        lp = load_op_and_grad(pth, "ot_network"; multiple=true)
        if returnall
            P, loss = lp(a, b, M, iter)
            P = set_shape(P, (length(a), length(b)))
            P, loss
        end
        return lp(a, b, M, iter)[2]
    else
        METHODIDS = Dict(
            "sinkhorn"=>0,
            "greenkhorn"=>1
        )

        if !haskey(METHODIDS, method)
            error("$method not implemented")
        end
        METHODID = METHODIDS[method]

        sk = load_system_op("sinkhorn_knopp"; multiple=true)
        if returnall
            return sk(a,b,M,reg,iter,tol, constant(METHODID))
        end
        return sk(a,b,M,reg,iter,tol, constant(METHODID))[2]
    end
end

@doc raw"""
    emd(a::Union{PyObject, Array{Float64}}, b::Union{PyObject, Array{Float64}}, M::Union{PyObject, Array{Float64}};
    iter::Int64 = 1000, tol::Float64 = 1e-9, returnall::Bool=false)

Computes the Earth Mover's Distance, which is defined as

$$D(M) = \sum_{i=1}^m \sum_{j=1}^n M_{ij} d_{ij}$$

Here $M \in \mathbb{R}^{m\times n}$ is the ground distance matrix. The algorithm solves the following optimization problem

$$\begin{aligned}\min_{M} &\ D(M)\\\text{s.t.} & \ \sum_{i=1}^m M_{ij} = b_j\\ &\ \sum_{j=1}^n M_{ij} = a_i \end{aligned}$$

The internal solver for the optimization problem is a netflow solver. The algorithm requires $\sum_i a_i = \sum_j b_j = 1$.
"""
function emd(a::Union{PyObject, Array{Float64}}, b::Union{PyObject, Array{Float64}}, M::Union{PyObject, Array{Float64}};
    iter::Int64 = 1000, tol::Float64 = 1e-9, returnall::Bool=false)
    return sinkhorn(a, b, M, reg = 0.0, iter = iter, tol = tol, method = "lp",
            returnall = returnall)
end

"""
    ot_plot1D(a::Array{Float64, 1}, b::Array{Float64, 1}, M::Array{Float64, 2})

Plots the optimal transport matrix for 1D distributions.
"""
function ot_plot1D(a::Array{Float64, 1}, b::Array{Float64, 1}, M::Array{Float64, 2})
    @assert length(a) == size(M, 1)
    @assert length(b) == size(M, 2)
    pl = require_import(:PyPlot)
    pl.figure(4, figsize = (5, 5))

    xa = 0:size(M, 1)-1
    xb = 0:size(M, 2)-1

    ax1 = pl.subplot(222)
    pl.plot(xb, b, "r")
    pl.yticks(())
    pl.title("Target")

    ax2 = pl.subplot(223)
    pl.plot(a, xa, "b")
    pl.gca().invert_xaxis()
    pl.gca().invert_yaxis()
    pl.xticks(())
    pl.title("Source")

    pl.subplot(224, sharex=ax1, sharey=ax2)
    pl.imshow(M, interpolation="nearest")
    pl.axis("off")

    pl.xlim((0, size(M, 2)-1))
    pl.tight_layout()
    pl.subplots_adjust(wspace=0., hspace=0.2)

end

@doc raw"""
    empirical_sinkhorn(x::Union{PyObject, Array{Float64}}, y::Union{PyObject, Array{Float64}};
        reg::Union{PyObject,Float64} = 1.0, iter::Int64 = 1000, tol::Float64 = 1e-9, method::String="sinkhorn", dist::Function=dist, returnall::Bool=false)

Computes the empirical Sinkhorn distance with sinkhorn algorithm. Here $x$ and $y$ are samples from two distributions.

- `reg` (default = 1.0): entropy regularization parameter
- `tol` (default = 1e-9), `iter` (default = 1000): tolerance and max iterations for the Sinkhorn algorithm
- `dist` (default = 2): Integer or Function, the distance function between two samples; if `dist` is integer, $L-dist$ norm is used.
- `returnall`: returns (`TransportMatrix`, `Loss`) if true; otherwise, only `Loss` is returned.

The implementation are adapted from https://github.com/rflamary/POT.
"""
function empirical_sinkhorn(x::Union{PyObject, Array{Float64}}, y::Union{PyObject, Array{Float64}};
    reg::Union{PyObject,Float64} = 1.0, iter::Int64 = 1000, tol::Float64 = 1e-9,
    method::String="sinkhorn", dist::Union{Integer,Function}=2, returnall::Bool=false, normalized::Bool = false)
    x, y = convert_to_tensor([x,y], [Float64, Float64])
    length(size(x))==1 && (x = reshape(x, :, 1))
    length(size(y))==1 && (y = reshape(y, :, 1))
    if isa(dist, Integer)
        d = copy(dist)
        dist = (x,y) -> ot_dist(x, y, d)
    end
    M = dist(x, y)
    if normalized
        M = M/maximum(M)
    end
    a = tf.ones(tf.shape(x)[1], dtype=tf.float64)/cast(Float64, tf.shape(x)[1])
    b = tf.ones(tf.shape(y)[1], dtype=tf.float64)/cast(Float64, tf.shape(y)[1])
    sinkhorn(a, b, M; reg=reg, iter=iter, tol=tol, method=method, returnall=returnall)
end

"""
    empirical_emd(x::Union{PyObject, Array{Float64}}, y::Union{PyObject, Array{Float64}};
        iter::Int64 = 1000, tol::Float64 = 1e-9, dist::Union{Integer,Function}=2, returnall::Bool=false)

Same as [`empirical_sinkhorn`](@ref), except that the Earth Mover Distance is computed.
"""
function empirical_emd(x::Union{PyObject, Array{Float64}}, y::Union{PyObject, Array{Float64}};
     iter::Int64 = 1000, tol::Float64 = 1e-9, dist::Union{Integer,Function}=2, returnall::Bool=false, normalized::Bool = false)
     empirical_sinkhorn(x, y; iter = iter, tol = tol, method = "lp",  returnall = returnall)
end

"""
    ot_dist(x::Union{PyObject, Array{Float64}}, y::Union{PyObject, Array{Float64}}, order::Union{Int64, PyObject}=2)

Computes the distance function with norm `order`. `dist` returns a ``n\\times m`` matrix, where ``x\\in \\mathbb{R}^{n\\times d}`` and
``y\\in \\mathbb{R}^{m\\times d}``, and the return ``M\\in \\mathbb{R}^{n\\times m}``
```math
M_{ij} = ||x_i - y_j||_{o}
```
"""
function ot_dist(x::Union{PyObject, Array{Float64}}, y::Union{PyObject, Array{Float64}}, order::Union{Int64, PyObject}=2)
    if !(isa(x, PyObject) || isa(y, PyObject) || isa(order, PyObject))
        x = reshape(x, size(x, 1), :)
        y = reshape(y, size(y, 1), :)
        M = zeros(size(x, 1), size(y, 1))
        for i = 1:size(x, 1)
            for j = 1:size(y, 1)
                M[i, j] = norm(x[i,:]-y[j,:], order)
            end
        end
        return M
    end
    x = convert_to_tensor(x)
    y = convert_to_tensor(y)
    order = convert_to_tensor(order, dtype=Float64)
    x = tf.expand_dims(x, axis=1)
    y = tf.expand_dims(y, axis=0)
    sum(abs(x-y)^order, dims=3)^(1.0/order)
end

"""
    dtw(s::Union{PyObject, Array{Float64}}, t::Union{PyObject, Array{Float64}},
        use_fast::Bool = false)

Computes the dynamic time wrapping (DTW) distance between two time series `s` and `t`.
Returns the distance and path. `use_fast` specifies whether fast algorithm is used. Note
fast algorithm may not be accurate.
"""
function dtw(s::Union{PyObject, Array{Float64}}, t::Union{PyObject, Array{Float64}},
            use_fast::Bool = false)
    use_fast = Int32(use_fast)
    pth = install("FastDTW")
    dtw_ = load_op_and_grad(pth, "dtw"; multiple=true)
    s,t, use_fast = convert_to_tensor([s,t,use_fast], [Float64,Float64, Int32])
    cost, path = dtw_(s,t,use_fast)
    return cost, path + 1
end


text:
include("read_expres_data_101501.jl")

#df_for_one_order_of_one_spectra = EXPRES.find_ranges_with_tellurics_in_order(all_spectra[5], 81)
#df_for_one_spectra = EXPRES.find_ranges_with_tellurics(all_spectra[5])

# Overwrite barycentric wavelengths with observatory-frame wavelengths
using FITSIO
map(obsid->all_spectra[obsid].λ .= FITSIO.read(FITS(df_files_use.Filename[obsid])["optimal"],"wavelength"), 1:length(all_spectra) )

wavelength_ranges_with_tellurics = EchelleInstruments.find_ranges_with_tellurics(all_spectra, min_Δv_clean=EXPRES.default_min_Δv_clean)

#using CSV
#CSV.write(joinpath(output_dir,"telluric_ranges.csv"),wavelength_ranges_with_tellurics)

min_λ_min = minimum(map(i->minimum(all_spectra[i].λ[.!isnan.(all_spectra[i].metadata[:tellurics]) .& (all_spectra[i].metadata[:tellurics] .== 1) ]), 1:length(all_spectra) ))
max_λ_max = maximum(map(i->maximum(all_spectra[i].λ[.!isnan.(all_spectra[i].metadata[:tellurics]) .& (all_spectra[i].metadata[:tellurics] .== 1) ]), 1:length(all_spectra) ))

wavelength_ranges_without_tellurics = EchelleInstruments.make_ranges_without_tellurics(wavelength_ranges_with_tellurics,
      λ_start = min_λ_min, λ_stop = max_λ_max, min_Δv= 1000)

wavelength_ranges_without_tellurics2 = EchelleInstruments.make_ranges_without_tellurics(wavelength_ranges_with_tellurics,
            λ_start = min_λ_min, λ_stop = max_λ_max, min_Δv= 1000, max_Δv = 1000000)

#CSV.write(joinpath(output_dir,"telluric_free_ranges_min_Δv=1000.csv"),wavelength_ranges_without_tellurics)

# Useful to compare how much extra information we'd get by reducing Δv
mapreduce(r->(r.lambda_hi-r.lambda_lo)*2/(r.lambda_hi+r.lambda_lo),+,eachrow(wavelength_ranges_without_tellurics))

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
using Random
using Test

using MultiAgentPOMDPs
using MultiAgentSysAdmin
using POMDPs

@testset "MultiAgentSysAdmin.jl" begin
    # Write your tests here.
    for sysfn in [UniSysAdmin, BiSysAdmin, RingofRingSysAdmin, RandomSysAdmin]
        for g in (true, false)
            sys = sysfn{g}()

            @test sys isa sysfn
            s = rand(initialstate(sys))
            a = ones(Int, n_agents(sys))
            @inferred @gen(:sp, :r)(sys, s, a, MersenneTwister(42))
        end
    end
end

max_stars_repo_name: JuliaPOMDP/MultiAgentSysAdmin

max_issues_repo_name: JuliaPOMDP/MultiAgentSysAdmin

max_forks_repo_name: JuliaPOMDP/MultiAgentSysAdmin


text:
module NiLang

using Reexport
@reexport using NiLangCore
import NiLangCore: ⊕, ⊖

using FixedPointNumbers: Q20f43, Fixed
export Fixed43
const Fixed43 = Q20f43

include("utils.jl")
include("vars.jl")
include("instructs.jl")
include("ulog.jl")
include("stack.jl")
include("complex.jl")
include("autobcast.jl")

include("autodiff/autodiff.jl")

include("stdlib/stdlib.jl")

include("deprecations.jl")

export AD
end # module


text:
## ----------------------------------------------------------------------------
# proj 2D
let
    method = ME_MAX_POL
    biom_ider = iJR.BIOMASS_IDER

    ps_pool = Dict()
    for exp in EXPS

        datfile = ME_INDEX[method, :DFILE, exp]
        dat = deserialize(datfile)

        model = dat[:model]

        for Hd_ider in FLX_IDERS

            # 2D Projection
            p = plot(;title = string("Heerden2013, exp", exp),
                xlabel = string(biom_ider), ylabel = string(Hd_ider),
                legend = :right
            )
            proj = DAT[method, :ep, :proj, Hd_ider, exp]
            ChP.plot_projection2D!(p, proj; l = 50)

            # bounds
            lb, ub = DAT[method, :bounds, :flx, Hd_ider, exp]
            hline!(p, [lb]; lw = 3,
                label = "fva lb",
                color = :blue, ls = :solid
            )
            hline!(p, [ub]; lw = 3,
                label = "fva ub",
                color = :red, ls = :solid
            )

            # EXPERIMENTAL FLXS
            exp_biom = DAT[method, :Hd, :flx, biom_ider, exp]
            exp_exch = DAT[method, :Hd, :flx, Hd_ider, exp]
            scatter!(p, [exp_biom], [exp_exch];
                m = 8, color = :red, label = "exp"
            )

            # MAXENT FLXS
            ep_biom = DAT[method, :ep, :flx, biom_ider, exp]
            ep_biom_err = DAT[method, :eperr, :flx, biom_ider, exp]
            ep_exch = DAT[method, :ep, :flx, Hd_ider, exp]
            ep_exch_err = DAT[method, :eperr, :flx, Hd_ider, exp]
            scatter!(p, [ep_biom], [ep_exch];
                xerr = [ep_biom_err], yerr = [ep_exch_err],
                m = 8, color = :blue, label = "maxent"
            )

            # mysavefig(p, "polytope"; Hd_ider, exp, method)
            ps_pool[(exp, Hd_ider)] = deepcopy(p)
        end
    end

    # collect
    for exp in EXPS
        ps = Plots.Plot[ps_pool[(exp, Hd_ider)] for Hd_ider in FLX_IDERS]
        mysavefig(ps, "polytope"; exp, method)
    end

    for Hd_ider in FLX_IDERS
        ps = Plots.Plot[ps_pool[(exp, Hd_ider)] for exp in EXPS]
        mysavefig(ps, "polytope"; Hd_ider, method)
    end
end


text:
using Combinatorics

mutable struct Moon
    pos::Vector{Int}
    v::Vector{Int}
end
Moon(pos) = Moon(pos, zeros(Int, 3))
Base.show(io::IO, m::Moon) = print("pos=<x = ",m.pos[1],", y= ",m.pos[2],", z= ",m.pos[3],">, vel=<x= ",m.v[1],", y= ",m.v[2],", z= ",m.v[3],">\n")
# Base.:(==)(a::Moon, b::Moon) = a.pos == b.pos && a.v == b.v
Base.:(==)(a::Moon, b::Moon) = hash(a) == hash(b)
Base.hash(m::Moon, h::UInt64) = hash((m.pos, m.v), h)

function energy(m::Moon)::Int
    return sum(abs.(m.pos)) * sum(abs.(m.v))
end

function create_moons(positions::Array{Vector{Int}})::Vector{Moon}
    moons = Moon[]
    for i in 1:length(positions)
        push!(moons, Moon(positions[i]))
    end
    return moons
end

function create_moons(positions::Array{})
    moons = Moon[]
    for i in 1:size(positions,1)
        push!(moons, Moon(positions[i,:]))
    end
    return moons
end

function simulate(moons::Vector{Moon})
    for pair in combinations(moons, 2)
        a, b = pair
        a.v = a.v + (a.pos .< b.pos) - (a.pos .> b.pos)
        b.v = b.v + (b.pos .< a.pos) - (b.pos .> a.pos)
    end
    for moon in moons
        moon.pos = moon.pos + moon.v
    end
end

function simulate(moons::Vector{Moon}, n::Int)
    for i = 1:n
        simulate(moons)
    end
end

function simulate(pos::Array{}, v::Array{})
    for pair in combinations(1:4, 2)
        a, b = pair
        v[a,:] = v[a,:] + (pos[a,:] .< pos[b,:]) - (pos[a,:] .> pos[b,:])
        v[b,:] = v[b,:] + (pos[b,:] .< pos[a,:]) - (pos[b,:] .> pos[a,:])
    end
    return pos+v, v
end

function simulate(pos::Array{}, v::Array{}, n::Int)
    for i in 1:n
        pos, v = simulate(pos, v)
    end
    return pos, v
end

function calculate_energy(positions::Array{Int,2}, n::Int)
    moons = create_moons(positions)
    # println(positions)
    simulate(moons, n)
    # println(moons)
    println("Total energy: ", sum(energy.(moons)))
end

# example 1
positions1 = [-1 0 2; 2 -10 -7; 4 -8 8; 3 5 -1]
# example 2
positions2 = [-8 -10 0; 5 5 10; 2 -7 3; 9 -8 -3]
# puzzle input
positions = [-16 15 -9; -14 5 4; 2 0 6; -3 18 9]

# velocities = zeros(Int,4,3)

calculate_energy(positions1, 10)
calculate_energy(positions2, 100)
calculate_energy(positions, 1000)
# function filter_history(x)
#     println(x)
# end

function find_periods(positions::Array{Int,2})
    # moons = create_moons(positions)
    velocities = zeros(Int,4,3)
    counter = 0
    history = Array{UInt64}(undef, 3, 1000000)
    #history[counter+1] = hash.(eachcol(vcat(positions, velocities)))
    history[:, counter+1] = hash.(eachcol(vcat(positions, velocities)))
    periods = [0,0,0]

    while true
        counter = counter + 1
        positions, velocities = simulate(positions, velocities)
        entry = hash.(eachcol(vcat(positions, velocities)))
        for axis in 1:3
            if periods[axis] == 0
                res = findfirst(x -> x == entry[axis], history[axis,1:counter])
                if res != nothing
                    println(counter, ": period found for axis ",axis, ": ", counter,"-",(res-1),"=", counter-(res-1))
                    periods[axis] = counter-(res-1)
                end
            end
        end
        if sum(periods .> 0) == 3
            break
        end
        if counter % 5000 == 0
            println(counter)
        end
        # history = hcat(history, entry)
        history[:,counter+1] = entry
    end

    return lcm(periods)
end

@time println(find_periods(positions1))
@time println(find_periods(positions2))
@time println(find_periods(positions))


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, "../.."))
Pkg.instantiate()

# ## visualizer
vis = Visualizer()
open(vis)

# ## setup
using Dojo
using IterativeLQR
using LinearAlgebra
using FiniteDiff
using DojoEnvironments
using JLD2

# ## scripts
include(joinpath(module_dir(), "examples/policy/methods/continuation.jl"))
include(joinpath(module_dir(), "examples/policy/methods/tvlqr.jl"))
include(joinpath(module_dir(), "DojoEnvironments/src",
    "quadruped/methods/template.jl"))

################################################################################
# ## system
################################################################################
gravity = -9.81
timestep = 0.02
friction_coefficient = 0.8
damper = 0.5
spring = 1.0
env = get_environment(:quadruped,
    representation=:minimal,
    timestep=timestep,
    contact_body=false,
    gravity=gravity,
    friction_coefficient=friction_coefficient,
    damper=damper,
    spring=spring,
    infeasible_control=true,
    vis=vis)

# ## dimensions
nx = env.num_states
nu = env.num_inputs
nu_infeasible = 6

################################################################################
# Load trajectory
################################################################################
file = JLD2.jldopen(joinpath(@__DIR__, "../data/trotting_forward.jld2"))
x_ref = file["x"]
u_ref = file["u"]
JLD2.close(file)

T = length(x_ref)

################################################################################
# TVLQR
################################################################################

K_tv, P_tv = tvlqr(x_ref, u_ref, env; q_tracking=q_tracking, r_tracking=r_tracking)

################################################################################
# tuning
################################################################################
function tracking_cost(x1, K; Q=Diagonal(ones(nx)))
    c = 0.0
    dcdx = [zeros(1, nx) for i = 1:T-1]
    dxdu = [zeros(nx, nu) for i = 1:T-1]
    dxpdx = [zeros(nx, nx) for i = 1:T-1]
    dudK = [spzeros(nu, nu*nx) for i = 1:T-1]
    dudx = [zeros(nu, nx) for i = 1:T-1]
    dcdK = [zeros(1, nu*nx) for i = 1:T-1]
    dxdK = [zeros(nx, nu*nx) for i = 1:T-1]
    dxdx = [zeros(nx, nx) for i = 1:T-1]
    dx = zeros(nx, nx)
    du = zeros(nx, nu)

    x = [zeros(nx) for i = 1:T]
    x[1] .= x1

    #Forward Pass
    for i = 1:T-1
        u_π = u_ref[i] + K[i] * (x_ref[i] - x[i])
        dynamics(x[i+1], env, x[i], u_π, zeros(0); gradients=true)
        dxpdx[i] .= env.dynamics_jacobian_state
        dxdu[i] .= env.dynamics_jacobian_input
        dudx[i] .= K[i]
        Δx = (x_ref[i] - x[i])
        for j = 1:nx
            dudK[i][:,(j-1)*nu .+ (1:nu)] = Δx[j] * I(nu)
        end
        dxdx[i] = dxdu[i] * dudx[i]
    end

    # Backward pass
    for i = T-1:-1:1
        c += 0.5 * (x[i+1] - x_ref[i+1])' * Q * (x[i+1] - x_ref[i+1])
        dcdx[i] .= (x[i+1] - x_ref[i+1])' * Q
        # for j = i:i
        for j = i:T-1
            X = I(nx)
            for k = j:-1:i+1
                X *= dxdx[k]
            end
            dcdK[i] .+= dcdx[j] * X *dxdu[j] * dudK[j]
        end
    end
    return x, c, dcdK
end

x, c, dcdK = tracking_cost(x_ref[1], K_tv, Q=Diagonal(ones(nx)))

# K_start = deepcopy(K_tv)
for i = 1:1
    c = 0.0
    dcdK = [zeros(1, nu*nx) for i=1:T-1]
    for j = 1:3
        x1 = deepcopy(x_ref[1])
        x1[j] += 0.05
        xi, ci, dcdKi = tracking_cost(x1, K_start, Q=Diagonal(q_tracking))
        c += ci / 3
        for k = 1:T-1
            dcdK[k] += dcdKi[k]/3
        end
    end
    # x, c, dcdK = tracking_cost(x_ref[1], K_start, Q=Diagonal(q_tracking))
    @show i, c
    for i = 1:T-1
        K_start[i] -= 0.001 * reshape(dcdK[i], (nu,nx))
    end
end

################################################################################
# Save
################################################################################
JLD2.jldsave(joinpath(@__DIR__, "../data/tvlqr_gains.jld2"), K=K_tv)
JLD2.jldsave(joinpath(@__DIR__, "../data/tuned_gains.jld2"), K=K_start)
file = JLD2.jldopen(joinpath(@__DIR__, "../data/tuned_gains.jld2"))
file["K"]
JLD2.close(file)

DojoEnvironments.visualize(env, x, build=false)


text:
export CholeskySolver
export cholesky_sym,
    cholesky_to!, solve_sym, solve, solve_to!, update!, downdate!

mutable struct CholeskySolver{T}
    list::Intlist

    ata::SparseMatrix{T}
    ld::SparseMatrix{T}
    y::Vector{T}

    perm::PermuteTable
    norm::Vector{T}

    row_count::Int64
    column_count::Int64

    tolerance::Float64
    use_permutation::Bool
    use_normalization::Bool

    function CholeskySolver{T}(row_count::Int64, column_count::Int64) where {T}
        local cs = new()

        cs.row_count = row_count
        cs.column_count = column_count

        cs.list = Intlist(column_count)

        cs.tolerance = 1e-10

        cs.use_permutation = true
        cs.use_normalization = true

        return cs
    end
end

function cholesky_sym(sym::SparseMatrix{T}) where {T}
    @assert lower_symmetric == sym.layout

    local n = sym.row_count

    local g = DirectedGraph{T}(n)

    local columns = Vector{Int64}()
    clear!(sym.list)

    for i = 1:n
        for j = sym.rows[i]:(sym.rows[i+1]-1)
            local jj = sym.rows_columns[j]
            push!(sym.list, jj)
            push!(columns, jj)
        end

        while length(columns) > 0
            local jj = pop!(columns)

            if !has_connections(g, jj)
                continue
            end

            local p = g.start[jj]
            while NIL != p
                local r = g.vertices[p]
                p = g.next[p]

                if jj == r
                    continue
                end

                if !contains(sym.list, r)
                    push!(sym.list, r)
                    push!(columns, r)
                end
            end
        end

        while !is_empty(sym.list)
            connect!(g, T(0), pop!(sym.list), i)
        end
    end

    local c = from_graph(g, n, n)
    local ct = transpose(c)
    ct.layout = lower_triangle

    return ct
end

function cholesky_to!(
    sym::SparseMatrix{T},
    ld::SparseMatrix{T},
    tolerance::Float64 = 1e-10,
) where {T}
    @assert lower_symmetric == sym.layout
    @assert lower_triangle == ld.layout

    local n = sym.row_count
    local zero = T(0)

    local acc = zeros(T, n)

    @inbounds begin
        for j = 1:n
            for i = sym.columns[j]:(sym.columns[j+1]-1)
                acc[sym.columns_rows[i]] = sym.values[i]
            end

            for k = ld.rows[j]:(ld.rows[j+1]-1)
                if ld.rows_columns[k] >= j
                    break
                end

                local r = ld.rows_columns[k]
                local a = ld.values[ld.positions[k]] * ld.values[ld.columns[r]]

                for i = ld.positions[k]:(ld.columns[r+1]-1)
                    acc[ld.columns_rows[i]] -= a * ld.values[i]
                end
            end

            local d = ld.values[ld.columns[j]] = acc[j]
            if d <= zero
                d = tolerance
            end

            for k = ld.columns[j]+1:(ld.columns[j+1]-1)
                ld.values[k] = acc[ld.columns_rows[k]] / d
                acc[ld.columns_rows[k]] = T(0)
            end
        end
    end
end

function solve_to!(
    ld::SparseMatrix{T},
    b::Vector{T},
    result::Vector{T},
) where {T}
    @assert lower_triangle == ld.layout
    @assert length(b) == ld.row_count

    local y = solve_lower(ld, b)

    local z = solve_diag(ld, y)

    local x = solve_upper(ld, z)

    result[:] = x[:]
end

function solve_lower(ld::SparseMatrix{T}, b::Vector{T}) where {T}
    @assert lower_triangle == ld.layout
    @assert length(b) == ld.row_count

    local n = ld.row_count
    local y = zeros(T, n)

    @inbounds begin
        for i = 1:n
            local sum = T(0)

            for j = ld.rows[i]:(ld.rows[i+1]-2)
                sum += ld.values[ld.positions[j]] * y[ld.rows_columns[j]]
            end

            y[i] = b[i] - sum
        end
    end

    return y
end

function solve_diag(ld::SparseMatrix{T}, y::Vector{T}) where {T}
    @assert lower_triangle == ld.layout
    @assert length(y) == ld.row_count

    local n = ld.row_count
    local x = zeros(T, n)

    @inbounds begin
        for j = n:-1:1
            x[j] = y[j] / ld.values[ld.columns[j]]
        end
    end

    return x
end

function solve_upper(ld::SparseMatrix{T}, z::Vector{T}) where {T}
    @assert lower_triangle == ld.layout
    @assert length(z) == ld.row_count

    local n = ld.row_count
    local x = zeros(T, n)

    @inbounds begin
        for j = n:-1:1
            local sum = T(0)

            for i = (ld.columns[j]+1):(ld.columns[j+1]-1)
                sum += ld.values[i] * x[ld.columns_rows[i]]
            end

            x[j] = z[j] - sum
        end
    end

    return x
end

function mul_transposed_to!(
    a::SparseMatrix{T},
    b::Vector{T},
    c::Vector{T},
) where {T}
    @assert a.row_count == length(b)
    @assert a.column_count == length(c)

    fill!(c, T(0))

    @inbounds begin
        for i = 1:a.row_count
            for j = a.rows[i]:(a.rows[i+1]-1)
                c[a.rows_columns[j]] += b[i] * a.values[a.positions[j]]
            end
        end
    end
end

function solve_sym(cs::CholeskySolver{T}, a::SparseMatrix{T}) where {T}
    if cs.use_permutation
        cs.perm = PermuteTable(a)
    else
        cs.perm = PermuteTable(a.column_count)
    end

    cs.ata = sqr_sym(a, cs.perm)
    cs.ld = cholesky_sym(cs.ata)
    cs.y = zeros(T, a.column_count)
end

function solve(
    cs::CholeskySolver{T},
    a::SparseMatrix{T},
    b::Vector{T},
) where {T}
    sqr_to!(a, cs.ata, cs.perm)

    if cs.use_normalization
        cs.norm = norm!(cs.ata)
    else
        cs.norm = ones(T, cs.ata.column_count)
    end

    cholesky_to!(cs.ata, cs.ld, cs.tolerance)

    mul_transposed_to!(a, b, cs.y)
    cs.y[:] = cs.y[cs.perm.permuted]

    if cs.use_normalization
        cs.y .*= cs.norm
    end

    local x = zeros(T, cs.ld.row_count)
    solve_to!(cs.ld, cs.y, x)

    if cs.use_normalization
        x .*= cs.norm
    end

    x[:] = x[cs.perm.primary]

    return x
end

function update!(
    ld::SparseMatrix{T},
    u::SparseArray{T},
    pt::PermuteTable,
) where {T}
    @assert lower_triangle == ld.layout
    @assert ld.row_count == u.size

    local zero = T(0)
    local a = T(1)
    local b = zero
    local c = zero

    local vals = zeros(T, u.size)
    for i = 1:u.nnz
        vals[pt.primary[u.indices[i]]] = u.values[i]
    end

    for j = 1:u.size
        if zero == vals[j]
            continue
        end

        local jj = ld.columns[j]

        local diag = ld.values[jj]
        local x = vals[j]
        b = a + x * x / diag
        ld.values[jj] = diag * b / a
        c = x / (diag * b)
        a = b

        for i = (jj+1):(ld.columns[j+1]-1)
            local ii = ld.columns_rows[i]

            vals[ii] -= x * ld.values[i]
            ld.values[i] += c * vals[ii]
        end
    end
end

function downdate!(
    ld::SparseMatrix{T},
    u::SparseArray{T},
    tolerance::T,
    pt::PermuteTable,
) where {T}
    @assert lower_triangle == ld.layout
    @assert ld.row_count == u.size

    local zero = T(0)
    local one = T(1)

    local vals = zeros(T, u.size)
    for i = 1:u.nnz
        vals[pt.primary[u.indices[i]]] = u.values[i]
    end

    local d = diag(ld)
    for i = 1:u.size
        set_columnwise!(ld, i, i, one)
    end

    local p = solve_lower(ld, vals)

    for i = 1:u.size
        set_columnwise!(ld, i, i, d[i])
    end

    local sum = zero
    for i = 1:u.size
        sum += p[i] / get_columnwise(ld, i, i) * p[i]
    end

    local a = one - sum

    if a <= tolerance
        a = tolerance
    end

    for j = u.size:-1:1
        local jj = ld.columns[j]
        local d = ld.values[jj]

        local b = a + p[j] * p[j] / d
        ld.values[jj] = d * a / b
        local c = -p[j] / (d * a)
        vals[j] = p[j]

        a = b

        for i = (jj+1):(ld.columns[j+1]-1)
            local ii = ld.columns_rows[i]

            local v = vals[ii]
            vals[ii] += p[j] * ld.values[i]
            ld.values[i] += c * v
        end
    end
end

function update!(cs::CholeskySolver{T}, u::SparseArray{T}, v::T) where {T}
    @assert cs.ld.row_count == u.size

    local zero = T(0)
    for i = 1:u.nnz
        if zero != u.values[i]
            local ii = cs.perm.primary[u.indices[i]]
            cs.y[ii] += v * u.values[i] * cs.norm[ii]
        end
    end

    update!(cs.ld, u, cs.perm)

    local x = zeros(T, cs.ld.row_count)
    solve_to!(cs.ld, cs.y, x)

    if cs.use_normalization
        x .*= cs.norm
    end

    if cs.use_permutation
        x[:] = x[cs.perm.primary]
    end

    return x
end

function downdate!(cs::CholeskySolver{T}, u::SparseArray{T}, v::T) where {T}
    @assert cs.ld.row_count == u.size

    local zero = T(0)
    for i = 1:u.nnz
        if zero != u.values[i]
            local ii = cs.perm.primary[u.indices[i]]
            cs.y[ii] -= v * u.values[i] * cs.norm[ii]
        end
    end

    downdate!(cs.ld, u, cs.tolerance, cs.perm)

    local x = zeros(T, cs.ld.row_count)
    solve_to!(cs.ld, cs.y, x)

    if cs.use_normalization
        x .*= cs.norm
    end

    if cs.use_permutation
        x[:] = x[cs.perm.primary]
    end

    return x
end


text:
#
# Advent of Code
# Day 6: Signals and Noise
#
println("Advent of Code")
println("Day 6: Signals and Noise")
println("")

NCOLS = 8

# Data file
# Input data has been saved to this file
inputDataFileName = "day06-input.txt"

# Read the input data, parse it down to a list of "commands"
inputDataFile = open(inputDataFileName, "r")
inputData = readstring(inputDataFile)
close(inputDataFile)

inputDataLines = split(strip(inputData), '\n')

# Set up a vector of dicts - character:frequency
signalDict = Vector{Dict{Char,Int}}(NCOLS)
# Can we initialize the dicts without doing this loop?
for counter in 1:NCOLS
  signalDict[counter] = Dict{Char,Int}()
end

tic()

# For every line, loop through each character,
# and count the occurances
for line in inputDataLines
  for (counter, letter) in enumerate(line)
    try
      signalDict[counter][letter] += 1
    catch error
      if isa(error, KeyError)
        signalDict[counter][letter] = 1
      end
    end
  end
end

signal1 = ""
signal2 = ""

# for every column, find the most frequent (part 1) and
# least frequent character (part 2) and concatenate them
# together

for subDict in signalDict
  maxCount = 0
  maxChar = 'a'

  minCount = 99999
  minChar = 'a'

  for (k, v) in subDict
    if v > maxCount
      maxCount = v
      maxChar = k
    elseif v < minCount
      minCount = v
      minChar = k
    end
  end
  signal1 = signal1 * string(maxChar)
  signal2 = signal2 * string(minChar)
end
toc()

println("The part 1 signal is: $signal1")
println("The part 2 signal is: $signal2")