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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ f4fda666-7b9c-11eb-0304-716c5e710462
begin
    using Symbolics, ForwardDiff, Plots, PlutoUI, LaTeXStrings
	using ForwardDiff: jacobian
end

# ╔═╡ c09f092a-887e-11eb-1ccc-4367bf05602d
html"""

Section 1.6

 The Newton Method

Lecture Video

"""

# ╔═╡ d82f1eae-7b9c-11eb-24d8-e1dcb2eef71a
md"""
## Solving equations and finding inverse transformations using the Newton method
"""

# ╔═╡ e410c1d0-7ba1-11eb-394f-71dac89756b7
md"""
In science and engineering we often need to *solve systems of equations*.

If the equations are *linear* then linear algebra tells us a general method to solve them; these are now routinely applied to solve systems of millions of linear equations.

If the equations are *non*linear then things are less obvious. The main solution methods we know work by... reducing the nonlinear equations to a sequence of linear equations! They do this by *approximating* the function by a linear function and solving that to get a better solution, then repeating this operation as many times as necessary to get a *sequence* of increasingly better solutions. This is an example of an **iterative algorithm**.

A well-known and elegant method, which can be used in many different contexts, is the **Newton method**. It does, however, have the disadvantage that it requires derivatives of the function. This can be overcome using **automatic differentiation** techniques.

We will illustrate the Newton method using the `ForwardDiff.jl` package to carry out automatic differentiation, but we will also try to understand what's going on "under the hood".
"""

# ╔═╡ 5ea7344c-7ba2-11eb-2cc5-0bbdca218c82
md"""
## The Newton method in 1D

We would like to solve equations like $f(x) = g(x)$.
We rewrite that by moving all the terms to one side of the equation so that we can write $h(x) = 0$, with $h(x) := f(x) - g(x)$.

A point $x^*$ such that $h(x^*) = 0$ is called a **root** or **zero** of $h$.

The Newton method finds zeros, and hence solves the original equation.
"""

# ╔═╡ 8c0c412e-7c2f-11eb-1880-4f6c45d77597
md"""

The idea of the Newton method is to *follow the direction in which the function is pointing*! We do this by building a **tangent line** at the current position and following that instead, until it hits the $x$-axis.

Let's look at that visually first:

"""

# ╔═╡ ce44554e-847f-4129-8841-1a729dfa7a2e
md"""
n = $(@bind n2 Slider(0:10, show_value=true, default=0))
"""

# ╔═╡ 77ef0cfb-60db-4599-bec2-b65e99e5b246
md"""
x₀ = $(@bind x02 Slider(-10:10, show_value=true, default=6))
"""

# ╔═╡ 2445da24-7b9d-11eb-02bd-eb99a3d95a2e
md"""
n = $(@bind n Slider(0:10, show_value=true, default=0))
"""

# ╔═╡ 9addbcbe-7b9e-11eb-3e8c-fbab3be40e05
md"""
x₀ = $(@bind x0 Slider(-10:10, show_value=true, default=6))
"""

# ╔═╡ c0b4defe-7c2f-11eb-1913-bdb01d28a4a8
md"""
## Using symbolic calculations to understand derivatives and nonlinear maps
"""

# ╔═╡ 615aff3c-7c30-11eb-2ca8-9d2fdf299017
md"""
We can use Julia's new symbolic capabilities to understand what's going on with a nonlinear (polynomial) function.

Let's see what happens if we perturb a function $f$ around a point $z$ by a small amount $\eta$.
"""

# ╔═╡ 71efd6b0-7c30-11eb-0da7-0d4a5ab8f8ff
@variables z, η

# ╔═╡ a869e6c6-7c31-11eb-13c8-155d08be02eb
md"""
m = $(@bind m Slider(1:6, show_value=true))
"""

# ╔═╡ 6dc89964-7c30-11eb-0a41-8d97b210ed34
f(x) = x^m - 2;

# ╔═╡ d35e0cc8-7c30-11eb-28d3-17c9e221ea62
f′(x) = ForwardDiff.derivative(f, x);

# ╔═╡ 63dbf052-7c32-11eb-1062-5b3581d38f70
f(z)

# ╔═╡ 9371f930-7c30-11eb-1f77-c7f31b97ea26
f(z + η)

# ╔═╡ 9d778e36-7c30-11eb-1f4b-894af86a8f5d
md"""
When $\eta$	is small, $\eta^2$ is *very* small, so we can ignore it. We are left with terms that either don't contain $\eta$ (constants), or multiply $\eta$ (linear). The part that multiplies $\eta$ is the derivative:
"""

# ╔═╡ db26375a-7c30-11eb-066e-ab9e8ded3356
f′(z)

# ╔═╡ ea741018-7c30-11eb-3912-a50475e6ec49
f(z) + η*f′(z)

# ╔═╡ 389e990e-7c40-11eb-37c4-5ba0f59173b3
md"""
The derivative gives the "*linear* part" of the function. `ForwardDiff.jl`, and forward-mode automatic differentiation in general, effectively uses this (although not symbolically in this sense) to just propagate the linear part of each function through a calculation.
"""

# ╔═╡ 5123c038-7ba2-11eb-1be2-19f789b02c1f
md"""
## Mathematics of the Newton method
"""

# ╔═╡ 9bfafcc0-7ba2-11eb-1b67-e3a3803ead08
md"""
We can convert the idea of "following the tangent line" into equations as follows.
(You can also do so by just looking at the geometry in 1D, but that does not help in 2D.)
"""

# ╔═╡ f153b4b8-7ba0-11eb-37ec-4f1a3dbe20e8
md"""
Suppose we have a guess $x_0$ for the root and we want to find a (hopefully) better guess $x_1$.

Let's set $x_1 = x_0 + \delta$, where $x_1$ and $\delta$ are still unknown.

We want $x_1$ to be a root, so
"""

# ╔═╡ 9cfa9062-7ba0-11eb-3a93-197ac0287ab4
md"""
$$f(x_1) = f(x_0 + \delta) \simeq 0$$
"""

# ╔═╡ 1ba1ae44-7ba1-11eb-21ff-558c95446435
md"""
If we are already "quite close" to the root then $\delta$ should be small, so we can approximate $f$ using the tangent line:

$$f(x_0) + \delta \, f'(x_0) \simeq 0$$

and hence

$$\delta \simeq \frac{-f(x_0)}{f'(x_0)}$$

so that

$$x_1 = x_0 - \frac{f(x_0)}{f'(x_0)}$$

Now we can repeat so that

$$x_2 = x_1 - \frac{f(x_1)}{f'(x_1)}$$

and in general

$$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}.$$

This is the Newton method in 1D.
"""

# ╔═╡ ba570c4c-7ba2-11eb-2125-9f23e415a1dc
md"""
## Implementation in 1D
"""

# ╔═╡ d690f83a-7c2e-11eb-14d7-79a250deb473
function newton1D(f, x0)

	f′(x) = ForwardDiff.derivative(f, x)   # \prime<TAB>

	x = x0

	for i in 1:10
		x -= f(x) / f′(x)
	end

	return x

end

# ╔═╡ 2fb40dc6-7c2f-11eb-2469-8deb4db59b5c
newton1D(x -> x^2 - 2, 37.0)

# ╔═╡ 35791bca-7c2f-11eb-1cfb-8d5ebd0208cb
sqrt(2)

# ╔═╡ 1d7dd328-7c2d-11eb-2b35-bdbf5df686f0
md"""
## Symbolic derivative in 2D

"""

# ╔═╡ d44c73b4-7c3e-11eb-1302-8ba9039ae789
md"""
Let's see what happens when we perturb by small amounts $\delta$ in the $x$ direction and $\epsilon$ in the $y$ direction around the point $(a, b)$:
"""

# ╔═╡ fe742fec-7c3e-11eb-1f54-55cdf02a1574
md"""
p = $(@bind p Slider(0:0.01:1, show_value=true))
"""

# ╔═╡ 23536420-7c2d-11eb-20b0-9523f7a5f9d7
@variables a, b, δ, ϵ

# ╔═╡ 4dd2322c-7ba0-11eb-2b3b-af7c6c1d60a0
md"""
## Newton for transformations in 2 dimensions

$$T: \mathbb{R}^2 \to \mathbb{R}^2$$

"""

# ╔═╡ 5c9edb2c-7ba0-11eb-14f6-3d5e52123bc7
md"""
We want to find the inverse $T^{-1}(y)$, i.e. to solve the equation $T(x) = y$ for $x$.

We use the same idea as in 1D, but now in 2D:
"""

# ╔═╡ 80917990-7ba0-11eb-029a-dba981c52b58
md"""
$$T(x_0 + \delta) \simeq 0$$

$$T(x_0) + J \cdot \delta \simeq 0,$$

where $J := DT_{x_0}$ is the Jacobian matrix of $T$ at $x_0$, i.e. the best linear approximation of $T$ near to $x_0$.
"""

# ╔═╡ af887dea-7ba1-11eb-3b0d-6925756382a7
md"""
Hence $\delta$ is the solution of the system of linear equations
"""

# ╔═╡ b7dc4666-7ba1-11eb-32eb-fd3d720c2960
md"""
$$J \cdot \delta = -T(x_0)$$

Then we again construct the new approximation $x_1$ as $x_1 := x_0 + \delta$.
"""

# ╔═╡ c519704c-7ba1-11eb-12da-8b9b176daa0d
md"""
In 2D we have an explicit formula for the inverse of the matrix.
"""

# ╔═╡ e1afc6ca-7ba1-11eb-3fb9-ef3a7f82d750
md"""
## Implementation in 2D
"""

# ╔═╡ 1db66b0e-7ba4-11eb-2157-d5a399a73b1f
function newton2D_step(T, x)

	J = ForwardDiff.jacobian(T, x)   # should use StaticVectors

	δ = J \ T(x)   # J^(-1) * T(x)

	return x - δ
end

# ╔═╡ 923bde64-7ba4-11eb-21e9-a11993aaab2e
"Looks for x such that T(x) = 0"
function newton2D(T, x0, n=10)

	x = x0

	for i in 1:n
		x = newton2D_step(T, x)
	end

	return x
end

# ╔═╡ 61905ae0-7ba6-11eb-0773-17e9aa4e9991
md"""
Remember that Newton is designed to look for *roots*, i.e. places where $T(x) = 0$.
We want $T(x) = y$, so we need another layer:
"""

# ╔═╡ ff8b6aec-7ba5-11eb-0d83-19803b1bdda7
"Looks for x such that f(x) = y, i.e. f(x) - y = 0"
function inverse(f, y, x0=[0, 0])
	return newton2D(x -> f(x) - y, x0)
end

# ╔═╡ 2e2e5f0e-7c31-11eb-0da7-770b07ee6202
inverse(f) = y -> inverse(f, y)

# ╔═╡ 1b77fada-7b9d-11eb-3266-ebb3895cb76a
straight(x0, y0, x, m) = y0 + m * (x - x0)

# ╔═╡ f25af026-7b9c-11eb-1f11-77a8b06b2d71
function standard_Newton(f, n, x_range, x0, ymin=-10, ymax=10)

    f′ = x -> ForwardDiff.derivative(f, x)

	p = plot(f, x_range, lw=3, ylim=(ymin, ymax), legend=:false, size=(400, 300))

	hline!([0.0], c="magenta", lw=3, ls=:dash)
	scatter!([x0], [0], c="green", ann=(x0, -5, L"x_0", 10))

	for i in 1:n

		plot!([x0, x0], [0, f(x0)], c=:gray, alpha=0.5)
		scatter!([x0], [f(x0)], c=:red)
		m = f′(x0)

		plot!(x_range, [straight(x0, f(x0), x, m) for x in x_range],
			  c=:blue, alpha=0.5, ls=:dash, lw=2)

		x1 = x0 - f(x0) / m

		scatter!([x1], [0], c="green", ann=(x1, -5, L"x_%$i", 10))

		x0 = x1

	end

	p |> as_svg

end

# ╔═╡ ecb40aea-7b9c-11eb-1476-e54faf32d91c
let
	f(x) = x^2 - 2

	standard_Newton(f, n2, -1:0.01:10, x02, -10, 70)
end

# ╔═╡ ec6c6328-7b9c-11eb-1c69-dba12ae522ad
let
	f(x) = 0.2x^3 - 4x + 1

	standard_Newton(f, n, -10:0.01:10, x0, -10, 70)
end

# ╔═╡ 515c23b6-7c2d-11eb-28c9-1b1d92eb4ba0
T(α) = ((x, y),) -> [x + α*y^2, y + α*x^2]

# ╔═╡ 09b97be8-7c2e-11eb-05fd-65bbd097afb8
jacobian(T(p), [a, b]) .|> Text

# ╔═╡ 18ce2fac-7c2e-11eb-03d2-b3a674621662
jacobian(T(p), [a, b]) * [δ, ϵ]

# ╔═╡ 395fd8e2-7c31-11eb-1933-dd719fa0cd22
md"""
α = $(@bind α Slider(0.0:0.01:1.0, show_value=true))
"""

# ╔═╡ 02b1b470-7c31-11eb-28f4-411956f73f12
T(α)( [0.3, 0.4] )

# ╔═╡ 07a754da-7c31-11eb-0394-4bef4d79fc30
inverse(T(α))( [0.3, 0.4] )

# ╔═╡ 5faa2784-7c31-11eb-34f1-3f8224dbdbde
( T(α) ∘ inverse(T(α)) )( [0.3, 0.4] )

# ╔═╡ ee91563e-7c3e-11eb-3f65-1f336073869a
md"""
## Appendix
"""

# ╔═╡ 786f8e78-7c2d-11eb-1bb8-c5cb2e349f45
expand(ex) = simplify(ex, expand=true)

# ╔═╡ 98158a38-7c30-11eb-0796-2335e97ec6d0
expand( f(z + η) )

# ╔═╡ e18f2470-7c31-11eb-2b74-d59d00d20ba4
expand( f(z + η) ) - ( f(z) + η*f′(z) )

# ╔═╡ 3828b94c-7c2d-11eb-2e01-79038b0f5226
image = expand.(T(p)( [ (a + δ), (b + ϵ) ] ))

# ╔═╡ ed605b90-7c3e-11eb-34e9-776a05a177dd
image - T(p)([a, b])

# ╔═╡ 35b5c5c6-7c3f-11eb-2723-4b406a809114
simplify.(expand.(image - T(p)([a, b]) - jacobian(T(p), [a, b]) * [δ, ϵ]))

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"

[compat]
ForwardDiff = "~0.10.19"
LaTeXStrings = "~1.2.1"
Plots = "~1.21.3"
PlutoUI = "~0.7.9"
Symbolics = "~3.3.1"
"""

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

[[AbstractTrees]]
git-tree-sha1 = "03e0550477d86222521d254b741d470ba17ea0b5"
uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
version = "0.3.4"

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

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

[[ArrayInterface]]
deps = ["IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "d84c956c4c0548b4caf0e4e96cf5b6494b5b1529"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "3.1.32"

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

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

[[Bijections]]
git-tree-sha1 = "705e7822597b432ebe152baa844b49f8026df090"
uuid = "e2ed5e7c-b2de-5872-ae92-c73ca462fb04"
version = "0.1.3"

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

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "4ce9393e871aca86cc457d9f66976c3da6902ea7"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.4.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"

[[Combinatorics]]
git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860"
uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
version = "1.0.2"

[[CommonSolve]]
git-tree-sha1 = "68a0743f578349ada8bc911a5cbd5a2ef6ed6d1f"
uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
version = "0.2.0"

[[CommonSubexpressions]]
deps = ["MacroTools", "Test"]
git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
version = "0.3.0"

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

[[CompositeTypes]]
git-tree-sha1 = "d5b014b216dc891e81fea299638e4c10c657b582"
uuid = "b152e2b5-7a66-4b01-a709-34e65c35f657"
version = "0.1.2"

[[ConstructionBase]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f74e9d5388b8620b4cee35d4c5a618dd4dc547f4"
uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
version = "1.3.0"

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

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

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

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

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

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

[[DiffResults]]
deps = ["StaticArrays"]
git-tree-sha1 = "c18e98cba888c6c25d1c3b048e4b3380ca956805"
uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
version = "1.0.3"

[[DiffRules]]
deps = ["NaNMath", "Random", "SpecialFunctions"]
git-tree-sha1 = "3ed8fa7178a10d1cd0f1ca524f249ba6937490c0"
uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
version = "1.3.0"

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

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

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

[[DomainSets]]
deps = ["CompositeTypes", "IntervalSets", "LinearAlgebra", "StaticArrays", "Statistics"]
git-tree-sha1 = "2655d0dd8bec4e01cbe903e5faa4617de40be779"
uuid = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
version = "0.5.5"

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

[[DynamicPolynomials]]
deps = ["DataStructures", "Future", "LinearAlgebra", "MultivariatePolynomials", "MutableArithmetics", "Pkg", "Reexport", "Test"]
git-tree-sha1 = "05b68e727a192783be0b34bd8fee8f678505c0bf"
uuid = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
version = "0.3.20"

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

[[EllipsisNotation]]
deps = ["ArrayInterface"]
git-tree-sha1 = "8041575f021cba5a099a456b4163c9a08b566a02"
uuid = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
version = "1.1.0"

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

[[ExprTools]]
git-tree-sha1 = "b7e3d17636b348f005f11040025ae8c6f645fe92"
uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
version = "0.1.6"

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

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

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

[[ForwardDiff]]
deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "NaNMath", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
git-tree-sha1 = "b5e930ac60b613ef3406da6d4f42c35d8dc51419"
uuid = "f6369f11-7733-5829-9624-2563aa707210"
version = "0.10.19"

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

[[Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

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

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

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

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

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

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

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

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

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

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

[[IfElse]]
git-tree-sha1 = "28e837ff3e7a6c3cdb252ce49fb412c8eb3caeef"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.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"

[[IntervalSets]]
deps = ["Dates", "EllipsisNotation", "Statistics"]
git-tree-sha1 = "3cc368af3f110a767ac786560045dceddfc16758"
uuid = "8197267c-284f-5f27-9208-e0e47529a953"
version = "0.5.3"

[[IrrationalConstants]]
git-tree-sha1 = "f76424439413893a832026ca355fe273e93bce94"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.0"

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

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

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

[[LabelledArrays]]
deps = ["ArrayInterface", "LinearAlgebra", "MacroTools", "StaticArrays"]
git-tree-sha1 = "bdde43e002847c34c206735b1cf860bc3abd35e7"
uuid = "2ee39098-c373-598a-b85f-a56591580800"
version = "1.6.4"

[[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 = "0b4a5d71f3e5200a7dff793393e09dfc2d874290"
uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
version = "3.2.2+1"

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

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

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

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

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

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

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

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

[[LogExpFunctions]]
deps = ["ChainRulesCore", "DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "34dc30f868e368f8a17b728a1238f3fcda43931a"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.3"

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

[[MultivariatePolynomials]]
deps = ["DataStructures", "LinearAlgebra", "MutableArithmetics"]
git-tree-sha1 = "45c9940cec79dedcdccc73cc6dd09ea8b8ab142c"
uuid = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
version = "0.3.18"

[[MutableArithmetics]]
deps = ["LinearAlgebra", "SparseArrays", "Test"]
git-tree-sha1 = "3927848ccebcc165952dc0d9ac9aa274a87bfe01"
uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
version = "0.2.20"

[[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 = "887579a3eb005446d514ab7aeac5d1d027658b8f"
uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
version = "1.3.5+1"

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

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

[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+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"

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

[[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 = "2dbafeadadcf7dadff20cd60046bba416b4912be"
uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
version = "1.21.3"

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

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

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

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

[[QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "12fbe86da16df6679be7521dfb39fbc861e1dc7b"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.1"

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

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

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

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

[[RecursiveArrayTools]]
deps = ["ArrayInterface", "ChainRulesCore", "DocStringExtensions", "LinearAlgebra", "RecipesBase", "Requires", "StaticArrays", "Statistics", "ZygoteRules"]
git-tree-sha1 = "00bede2eb099dcc1ddc3f9ec02180c326b420ee2"
uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
version = "2.17.2"

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

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

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

[[RuntimeGeneratedFunctions]]
deps = ["ExprTools", "SHA", "Serialization"]
git-tree-sha1 = "cdc1e4278e91a6ad530770ebb327f9ed83cf10c4"
uuid = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
version = "0.5.3"

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

[[SciMLBase]]
deps = ["ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "RecipesBase", "RecursiveArrayTools", "StaticArrays", "Statistics", "Tables", "TreeViews"]
git-tree-sha1 = "ff686e0c79dbe91767f4c1e44257621a5455b1c6"
uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
version = "1.18.7"

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

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

[[Setfield]]
deps = ["ConstructionBase", "Future", "MacroTools", "Requires"]
git-tree-sha1 = "fca29e68c5062722b5b4435594c3d1ba557072a3"
uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
version = "0.7.1"

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

[[SpecialFunctions]]
deps = ["ChainRulesCore", "LogExpFunctions", "OpenSpecFun_jll"]
git-tree-sha1 = "a322a9493e49c5f3a10b50df3aedaf1cdb3244b7"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.6.1"

[[Static]]
deps = ["IfElse"]
git-tree-sha1 = "a8f30abc7c64a39d389680b74e749cf33f872a70"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.3.3"

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

[[StatsFuns]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "46d7ccc7104860c38b11966dd1f72ff042f382e4"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.10"

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

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

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

[[SymbolicUtils]]
deps = ["AbstractTrees", "Bijections", "ChainRulesCore", "Combinatorics", "ConstructionBase", "DataStructures", "DocStringExtensions", "DynamicPolynomials", "IfElse", "LabelledArrays", "LinearAlgebra", "MultivariatePolynomials", "NaNMath", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArrays", "TermInterface", "TimerOutputs"]
git-tree-sha1 = "7db169562da2c70a476810ec6a5d84392b0f1127"
uuid = "d1185830-fcd6-423d-90d6-eec64667417b"
version = "0.15.4"

[[Symbolics]]
deps = ["ConstructionBase", "DiffRules", "Distributions", "DocStringExtensions", "DomainSets", "IfElse", "Latexify", "Libdl", "LinearAlgebra", "MacroTools", "NaNMath", "RecipesBase", "Reexport", "Requires", "RuntimeGeneratedFunctions", "SciMLBase", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArrays", "SymbolicUtils", "TreeViews"]
git-tree-sha1 = "0e4d48d9c416563e3b067759e3cd72066d4f491b"
uuid = "0c5d862f-8b57-4792-8d23-62f2024744c7"
version = "3.3.1"

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

[[TermInterface]]
git-tree-sha1 = "02a620218eaaa1c1914d228d0e75da122224a502"
uuid = "8ea1fca8-c5ef-4a55-8b96-4e9afe9c9a3c"
version = "0.1.8"

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

[[TimerOutputs]]
deps = ["ExprTools", "Printf"]
git-tree-sha1 = "209a8326c4f955e2442c07b56029e88bb48299c7"
uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
version = "0.5.12"

[[TreeViews]]
deps = ["Test"]
git-tree-sha1 = "8d0d7a3fe2f30d6a7f833a5f19f7c7a5b396eae6"
uuid = "a2a6695c-b41b-5b7d-aed9-dbfdeacea5d7"
version = "0.3.0"

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

[[ZygoteRules]]
deps = ["MacroTools"]
git-tree-sha1 = "9e7a1e8ca60b742e508a315c17eef5211e7fbfd7"
uuid = "700de1a5-db45-46bc-99cf-38207098b444"
version = "0.2.1"

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

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

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

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

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

[[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:
# ╟─c09f092a-887e-11eb-1ccc-4367bf05602d
# ╠═f4fda666-7b9c-11eb-0304-716c5e710462
# ╟─d82f1eae-7b9c-11eb-24d8-e1dcb2eef71a
# ╟─e410c1d0-7ba1-11eb-394f-71dac89756b7
# ╟─5ea7344c-7ba2-11eb-2cc5-0bbdca218c82
# ╟─8c0c412e-7c2f-11eb-1880-4f6c45d77597
# ╟─ce44554e-847f-4129-8841-1a729dfa7a2e
# ╟─77ef0cfb-60db-4599-bec2-b65e99e5b246
# ╠═ecb40aea-7b9c-11eb-1476-e54faf32d91c
# ╟─2445da24-7b9d-11eb-02bd-eb99a3d95a2e
# ╟─9addbcbe-7b9e-11eb-3e8c-fbab3be40e05
# ╠═ec6c6328-7b9c-11eb-1c69-dba12ae522ad
# ╟─c0b4defe-7c2f-11eb-1913-bdb01d28a4a8
# ╟─615aff3c-7c30-11eb-2ca8-9d2fdf299017
# ╠═71efd6b0-7c30-11eb-0da7-0d4a5ab8f8ff
# ╠═6dc89964-7c30-11eb-0a41-8d97b210ed34
# ╠═d35e0cc8-7c30-11eb-28d3-17c9e221ea62
# ╟─a869e6c6-7c31-11eb-13c8-155d08be02eb
# ╠═63dbf052-7c32-11eb-1062-5b3581d38f70
# ╠═9371f930-7c30-11eb-1f77-c7f31b97ea26
# ╠═98158a38-7c30-11eb-0796-2335e97ec6d0
# ╟─9d778e36-7c30-11eb-1f4b-894af86a8f5d
# ╠═db26375a-7c30-11eb-066e-ab9e8ded3356
# ╠═ea741018-7c30-11eb-3912-a50475e6ec49
# ╠═e18f2470-7c31-11eb-2b74-d59d00d20ba4
# ╟─389e990e-7c40-11eb-37c4-5ba0f59173b3
# ╟─5123c038-7ba2-11eb-1be2-19f789b02c1f
# ╟─9bfafcc0-7ba2-11eb-1b67-e3a3803ead08
# ╟─f153b4b8-7ba0-11eb-37ec-4f1a3dbe20e8
# ╟─9cfa9062-7ba0-11eb-3a93-197ac0287ab4
# ╟─1ba1ae44-7ba1-11eb-21ff-558c95446435
# ╟─ba570c4c-7ba2-11eb-2125-9f23e415a1dc
# ╠═d690f83a-7c2e-11eb-14d7-79a250deb473
# ╠═2fb40dc6-7c2f-11eb-2469-8deb4db59b5c
# ╠═35791bca-7c2f-11eb-1cfb-8d5ebd0208cb
# ╟─1d7dd328-7c2d-11eb-2b35-bdbf5df686f0
# ╟─d44c73b4-7c3e-11eb-1302-8ba9039ae789
# ╟─fe742fec-7c3e-11eb-1f54-55cdf02a1574
# ╠═23536420-7c2d-11eb-20b0-9523f7a5f9d7
# ╠═3828b94c-7c2d-11eb-2e01-79038b0f5226
# ╠═09b97be8-7c2e-11eb-05fd-65bbd097afb8
# ╠═18ce2fac-7c2e-11eb-03d2-b3a674621662
# ╠═ed605b90-7c3e-11eb-34e9-776a05a177dd
# ╠═35b5c5c6-7c3f-11eb-2723-4b406a809114
# ╟─4dd2322c-7ba0-11eb-2b3b-af7c6c1d60a0
# ╟─5c9edb2c-7ba0-11eb-14f6-3d5e52123bc7
# ╟─80917990-7ba0-11eb-029a-dba981c52b58
# ╟─af887dea-7ba1-11eb-3b0d-6925756382a7
# ╟─b7dc4666-7ba1-11eb-32eb-fd3d720c2960
# ╟─c519704c-7ba1-11eb-12da-8b9b176daa0d
# ╟─e1afc6ca-7ba1-11eb-3fb9-ef3a7f82d750
# ╠═1db66b0e-7ba4-11eb-2157-d5a399a73b1f
# ╠═923bde64-7ba4-11eb-21e9-a11993aaab2e
# ╟─61905ae0-7ba6-11eb-0773-17e9aa4e9991
# ╠═ff8b6aec-7ba5-11eb-0d83-19803b1bdda7
# ╠═2e2e5f0e-7c31-11eb-0da7-770b07ee6202
# ╟─1b77fada-7b9d-11eb-3266-ebb3895cb76a
# ╠═f25af026-7b9c-11eb-1f11-77a8b06b2d71
# ╠═515c23b6-7c2d-11eb-28c9-1b1d92eb4ba0
# ╟─395fd8e2-7c31-11eb-1933-dd719fa0cd22
# ╠═02b1b470-7c31-11eb-28f4-411956f73f12
# ╠═07a754da-7c31-11eb-0394-4bef4d79fc30
# ╠═5faa2784-7c31-11eb-34f1-3f8224dbdbde
# ╟─ee91563e-7c3e-11eb-3f65-1f336073869a
# ╠═786f8e78-7c2d-11eb-1bb8-c5cb2e349f45
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
module NumericalIntegration

using LinearAlgebra
using Logging
using Base.Iterators
using Interpolations

export integrate, cumul_integrate
export Trapezoidal, TrapezoidalEven, TrapezoidalFast, TrapezoidalEvenFast
export SimpsonEven, SimpsonEvenFast
export RombergEven
export IntegrationMethod

abstract type IntegrationMethod end

struct Trapezoidal         <: IntegrationMethod end
struct TrapezoidalEven     <: IntegrationMethod end
struct TrapezoidalFast     <: IntegrationMethod end
struct TrapezoidalEvenFast <: IntegrationMethod end
struct SimpsonEven         <: IntegrationMethod end # https://en.wikipedia.org/wiki/Simpson%27s_rule#Alternative_extended_Simpson.27s_rule
struct SimpsonEvenFast     <: IntegrationMethod end
struct RombergEven{T<:AbstractFloat} <: IntegrationMethod
    acc::T
end # https://en.wikipedia.org/wiki/Romberg%27s_method
RombergEven() = RombergEven(1e-12)

const HALF = 1//2

#documentation

"""
    integrate(x,y...)

Compute numerical integral of y(x) from x=x[1] to x=x[end]. Return a scalar of the same type as the input. If no method is supplied, use Trapezdoial.
"""
function integrate(x,y...) end

"""
    cumul_integrate(x,y...)

Compute cumulative numerical integral of y(x) from x=x[1] to x=x[end]. Return a vector with elements of the same type as the input. If no method is supplied, use Trapezdoial.
"""
function cumul_integrate(x,y...) end

#implementation

function _zero(x,y)
    ret = zero(eltype(x)) + zero(eltype(y))
end

function _zeros(x::AbstractVector,y::AbstractVector)
    ret = zeros(eltype(x),size(x)) + zeros(eltype(y),size(y))
end

"""
    integrate(x::AbstractVector, y::AbstractVector, ::Trapezoidal)

Use Trapezoidal rule. This is the default when no method is supplied.
"""
function integrate(x::AbstractVector, y::AbstractVector, ::Trapezoidal)
    @assert length(x) == length(y) "x and y vectors must be of the same length!"
    retval = _zero(x,y)
    for i in 1 : length(y)-1
        retval += (x[i+1] - x[i]) * (y[i] + y[i+1])
    end
    return HALF * retval
end

"""
    integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEven)

Use Trapezoidal rule, assuming evenly spaced vector x.
"""
function integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEven)
    @assert length(x) == length(y) "x and y vectors must be of the same length!"
    return (x[2] - x[1]) * (HALF * (y[1] + y[end]) + sum(y[2:end-1]))
end

"""
    integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalFast)

Use Trapezoidal rule. Unsafe method: no bound checking.
"""
function integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalFast)
    retval = _zero(x,y)
    @fastmath @simd for i in 1 : length(y)-1
        @inbounds retval += (x[i+1] - x[i]) * (y[i] + y[i+1])
    end
    return HALF * retval
end

"""
    integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEvenFast)

Use Trapezoidal rule, assuming evenly spaced vector x. Unsafe method: no bound checking.
"""
function integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEvenFast)
    retval = _zero(x,y)
    N = length(y) - 1
    @fastmath @simd for i in 2 : N
        @inbounds retval += y[i]
    end
    @inbounds return (x[2] - x[1]) * (retval + HALF*y[1] + HALF*y[end])
end

"""
    integrate(x::AbstractVector, y::AbstractVector, ::SimpsonEven)

Use Simpson's rule, assuming evenly spaced vector x.
"""
function integrate(x::AbstractVector, y::AbstractVector, ::SimpsonEven)
    @assert length(x) == length(y) "x and y vectors must be of the same length!"
    retval = (17*y[1] + 59*y[2] + 43*y[3] + 49*y[4] + 49*y[end-3] + 43*y[end-2] + 59*y[end-1] + 17*y[end]) / 48
    for i in 5 : length(y) - 1
        retval += y[i]
    end
    return (x[2] - x[1]) * retval
end

"""
    integrate(x::AbstractVector, y::AbstractVector, ::SimpsonEven)

Use Simpson's rule, assuming evenly spaced vector x.  Unsafe method: no bound checking.
"""
function integrate(x::AbstractVector, y::AbstractVector, ::SimpsonEvenFast)
    @inbounds retval = 17*y[1] + 59*y[2] + 43*y[3] + 49*y[4]
    @inbounds retval += 49*y[end-3] + 43*y[end-2] + 59*y[end-1] + 17*y[end]
    retval /= 48
    @fastmath @inbounds for i in 5 : length(y)-1
        retval += y[i]
    end
    @inbounds return (x[2] - x[1]) * retval
end

"""
    integrate(x::AbstractVector, y::AbstractMatrix, method; dims=2)

When y is an array, compute integral along dimension specified by dims (default 2: columns).
"""
function integrate(x::AbstractVector, y::AbstractMatrix, M::IntegrationMethod; dims=2)
    out = [integrate(x,selectdim(y,dims,j),M) for j=1:size(y,dims)]
    return out
end

function integrate(x::AbstractVector, y::AbstractVector, m::RombergEven)
    @assert length(x) == length(y) "x and y vectors must be of the same length!"
    @assert ((length(x) - 1) & (length(x) - 2)) == 0 "Need length of vector to be 2^n + 1"
    maxsteps::Integer = Int(log2(length(x)-1))
    rombaux = zeros(eltype(y), maxsteps, 2)
    prevrow = 1
    currrow = 2
    @inbounds h = x[end] - x[1]
    @inbounds rombaux[prevrow, 1] = (y[1] + y[end])*h*HALF
    @inbounds for i in 1 : (maxsteps-1)
        h *= HALF
        npoints = 1 << (i-1)
        jumpsize = div(length(x)-1, 2*npoints)
        c = 0.0
        for j in 1 : npoints
            c += y[1 + (2*j-1)*jumpsize]
        end
        rombaux[1, currrow] = h*c + HALF*rombaux[1, prevrow]
        for j in 2 : (i+1)
            n_k = 4^(j-1)
            rombaux[j, currrow] = (n_k*rombaux[j-1, currrow] - rombaux[j-1, prevrow])/(n_k - 1)
        end

        if i > maxsteps//3 && norm(rombaux[i, prevrow] - rombaux[i+1, currrow], Inf) < m.acc
            return rombaux[i+1, currrow]
        end

        prevrow, currrow = currrow, prevrow
    end
    finalerr = norm(rombaux[maxsteps-1, prevrow] - rombaux[maxsteps, currrow], Inf)
    @warn "RombergEven :: final step reached, but accuracy not: $finalerr > $(m.acc)"
    @inbounds return rombaux[maxsteps, prevrow]
end

@inline function _midpoints(x::AbstractVector{T}) where T
    length(x) == 1 && return x
    retval = Vector{T}(undef, length(x)-1)
    @simd for i in eachindex(retval)
        retval[i] = HALF * (x[i] + x[i+1])
    end
    return retval
end
@inline function _midpoints(x::AbstractRange)
    length(x) == 1 && return x
    Δx = HALF*step(x)
    return range(first(x)+Δx, stop=last(x)-Δx, length=length(x)-1)
end

function integrate(X::NTuple{N,AbstractVector}, Y::AbstractArray{T,N}, ::Trapezoidal) where {T,N}
    @assert length.(X) == size(Y)
    return integrate(X, Y, TrapezoidalFast())
end

function integrate(X::NTuple{N,AbstractVector}, Y::AbstractArray{T,N}, ::TrapezoidalFast) where {T,N}
    midpnts = map(_midpoints, X)
    Δ(x::AbstractVector) = length(x) > 1 ? diff(x) : 1
    Δxs = map(Δ, X)
    interp = LinearInterpolation(X,Y)
    f((Δx,x)) = prod(Δx)*interp(x...)
    return sum(f, zip(product(Δxs...), product(midpnts...)))
end

function integrate(X::NTuple{N,AbstractVector}, Y::AbstractArray{T,N}, ::TrapezoidalEvenFast) where {T,N}
    midpnts = map(_midpoints, X)
    Δ(x::AbstractVector) = length(x) > 1 ? x[2] - x[1] : 1
    Δx = prod(Δ, X)
    interp = LinearInterpolation(X,Y)
    f(x) = interp(x...)
    return Δx*sum(f, product(midpnts...))
end

function integrate(X::NTuple{N,AbstractRange}, Y::AbstractArray{T,N}, ::TrapezoidalEvenFast) where {T,N}
    midpnts = map(_midpoints, X)
    Δ(x::AbstractVector) = length(x) > 1 ? step(x) : 1
    Δx = prod(Δ, X)
    interp = LinearInterpolation(X,Y)
    f(x) = interp(x...)
    return Δx*sum(f, product(midpnts...))
end

#function integrate(X::Tuple{AbstractVector}, Y::AbstractVector{T}, M::IntegrationMethod) :: T where {T}
#    return integrate(X[1], Y, M)
#end
#"""
#    integrate(X::NTuple{N,AbstractVector}, Y::AbstractArray{T,N}, method, cache=nothing)
#
#Given an n-dimensional grid of values, compute the total integral along each dim
#"""
#function integrate(X::NTuple{N,AbstractVector}, Y::AbstractArray{T,N}, M::IntegrationMethod) :: T where {T,N}
#    n = last(size(Y))
#    cache = Vector{T}(undef, n)
#    x = X[1:N-1]
#    @inbounds for i in 1:n
#        cache[i] = integrate(x, selectdim(Y,N,i), M)
#    end
#    return integrate(X[end], cache, M)
#end

# cumulative integrals

"""
    cumul_integrate(x::AbstractVector, y::AbstractVector, ::Trapezoidal)

Use Trapezoidal rule. This is the default when no method is supplied.
"""
function cumul_integrate(x::AbstractVector, y::AbstractVector, ::Trapezoidal)
    @assert length(x) == length(y) "x and y vectors must be of the same length!"
    retarr = _zeros(x,y)
    for i in 2 : length(y)
        retarr[i] = retarr[i-1] + (x[i] - x[i-1]) * (y[i] + y[i-1])
    end
    return HALF * retarr
end

"""
    cumul_integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEven)

Use Trapezoidal rule, assuming evenly spaced vector x.
"""
function cumul_integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEven)
    @assert length(x) == length(y) "x and y vectors must be of the same length!"
    retarr = _zeros(x,y)
    for i in 2 : length(y)
        retarr[i] = retarr[i-1] + (y[i-1] + y[i])
    end
    return (x[2] - x[1]) * HALF * retarr
end

"""
    cumul_integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalFast)

Use Trapezoidal rule. Unsafe method: no bound checking.
"""
function cumul_integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalFast)
    retarr = _zeros(x,y)
    @fastmath for i in 2 : length(y) #not sure if @simd can do anything here
        @inbounds retarr[i] = retarr[i-1] + (x[i] - x[i-1]) * (y[i] + y[i-1])
    end
    return HALF * retarr
end

"""
    cumul_integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEvenFast)

Use Trapezoidal rule, assuming evenly spaced vector x. Unsafe method: no bound checking.
"""
function cumul_integrate(x::AbstractVector, y::AbstractVector, ::TrapezoidalEvenFast)
    retarr = _zeros(x,y)
    @fastmath for i in 2 : length(y)
        @inbounds retarr[i] = retarr[i-1] + (y[i] + y[i-1])
    end
    @inbounds return (x[2] - x[1]) * HALF * retarr
end

"""
    cumul_integrate(x::AbstractVector, y::AbstractMatrix, method; dims=2)

When y is an array, compute integral along each dimension specified by dims (default 2: columns)
"""
function cumul_integrate(x::AbstractVector, y::AbstractMatrix, M::IntegrationMethod; dims=2)
    return hcat([cumul_integrate(x,selectdim(y,dims,j),M) for j=1:size(y,dims)]...)
end

#default behaviour
integrate(x::AbstractVector, y::AbstractVector) = integrate(x, y, Trapezoidal())

integrate(x::AbstractVector, y::AbstractMatrix; dims=2) = integrate(x, y, Trapezoidal(); dims=dims)

integrate(X::NTuple, Y::AbstractArray) = integrate(X, Y, Trapezoidal())

cumul_integrate(x::AbstractVector, y::AbstractVector) = cumul_integrate(x, y, Trapezoidal())

cumul_integrate(x::AbstractVector, y::AbstractMatrix; dims=2) = cumul_integrate(x, y, Trapezoidal(); dims=dims)

end # module

ext: jl


text:
using Test
using QuantumOptics, CollectiveSpins

@testset "rotate" begin

cs = CollectiveSpins
e_dipole = [0, 0, 1]

function compare_rotate(;Tend=0., phi=0., theta=0., N=1., α=0., axis=[1.,0.,0.])
    system = SpinCollection(geometry.chain(0.9, N), e_dipole)
    T = Float64[0, Tend]
    state0_mf = cs.meanfield.blochstate(phi, theta, N)
    if abs(Tend)<1e-12
        stateT_mf = state0_mf
    else
        stateT_mf = CollectiveSpins.meanfield.timeevolution(T, system, state0_mf)[2][end]
    end
    state_mf = cs.meanfield.rotate(axis, α, stateT_mf)
    ρmf = cs.mpc.densityoperator(state_mf)

    state0_mpc = cs.mpc.blochstate(phi, theta, N)
    if abs(Tend)<1e-12
        stateT_mpc = state0_mpc
    else
        stateT_mpc = CollectiveSpins.mpc.timeevolution(T, system, state0_mpc)[2][end]
    end
    state_mpc = cs.mpc.rotate(axis, α, stateT_mpc)
    ρmpc = cs.mpc.densityoperator(state_mpc)

    Ψ₀ = CollectiveSpins.quantum.blochstate(phi,theta,N)
    ρ₀ = Ψ₀⊗dagger(Ψ₀)
    if abs(Tend)<1e-12
        ρT = ρ₀
    else
        ρT = CollectiveSpins.quantum.timeevolution(T, system, ρ₀)[2][end]
    end
    ρ = cs.quantum.rotate(axis, α, ρT)

    td_mpc_mf = QuantumOptics.tracedistance(ρmpc, ρmf)
    td_master_mf = QuantumOptics.tracedistance(ρ, ρmf)
    td_master_mpc = QuantumOptics.tracedistance(ρ, ρmpc)
    return td_mpc_mf, td_master_mf, td_master_mpc
end

td = compare_rotate(Tend=0., phi=0., theta=0., N=2, α=[1.0,0.3], axis=[1.,1.,0])
@test maximum(td) < 1e-12

td = compare_rotate(Tend=0., phi=0.3, theta=1., N=5, α=1.6, axis=[1.,1.,5.])
@test maximum(td) < 1e-12

angles = [1.2, 2.7]
td_mpc_mf, td_master_mf, td_master_mpc = compare_rotate(Tend=1., phi=0.3, theta=1., N=2, α=angles, axis=[2.,1.,5.])
@test td_mpc_mf < 0.05
@test td_master_mf < 0.05
@test td_master_mpc < 1e-6

angles = [1.2, 2.7, 2.1, 0.3, -0.5]
td_mpc_mf, td_master_mf, td_master_mpc = compare_rotate(Tend=1., phi=0.3, theta=1., N=5, α=angles, axis=[2.,1.,5.])
@test td_mpc_mf < 0.1
@test td_master_mf < 0.1
@test td_master_mpc < 0.01

end # testset


text:
using Test
using Forecast
using Random
using Logging

const tests = [
    "acf",
    "ccf",
    "datasets",
    "d",
    "loess",
    "sma",
    "stl",
    "hma"
]

printstyled("\nTest Summary List:\n", color=:underline)

Random.seed!(36)
Base.disable_logging(Base.CoreLogging.Error) # disable warnings

for t in tests
    @testset "Test $t" begin
        Random.seed!(36)
        include("$t.jl")
        println()
    end
end

println()


text:
function optimise(
        model::Model,
        target = XXX,
        num_params=length(model.params),
        max_steps=Int(2e4),
        loss=:integral,
        silent=false,
        cost_kwargs = Dict{Symbol, Any}(),
        error_on_abort=false,
    )

    ### Set up for optimisation
    cost_kwargs[:loss] = loss

    opt_kwargs = Dict(
        :SearchRange => (0.,1.),
        :NumDimensions => num_params,
        :Method => :adaptive_de_rand_1_bin_radiuslimited,
        :TraceInterval => 5.,
        :MaxSteps => max_steps,
        :TraceMode => silent ? :silent : :compact,
    )

    ### Do the optimisation
    res = bboptimize(
        p -> cost(SimpleParam(model, transform_param(model, p)), target; cost_kwargs...);
        opt_kwargs...)

    error_on_abort && res.stop_reason == "InterruptException()" && error("InterruptException from blackboxoptim.")

    ### Format info for storage in parameter type.
    lower_param_bound = transform_param(model, fill(0., length(model.params)))
    upper_param_bound = transform_param(model, fill(1., length(model.params)))

    parameter_set = transform_param(model, best_candidate(res))
    cost = best_fitness(res)

    cost_kwargs[:target] = target

    metadata = Dict(
        :cost => cost_kwargs,
        :opt => opt_kwargs,
        :lower_param_bound => lower_param_bound,
        :upper_param_bound => upper_param_bound,
        )

    ### Create parameter type instance
    p = Param(model, parameter_set, cost, metadata)

    return p
end


text:
#=##############################################################################
# DESCRIPTION
    Miscellaneous functions.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Feb 2021
  * License   : MIT
=###############################################################################

function remove_particles_lowstrength(sim, PFIELD, T, DT; crit=0.00001^2)

    for i in vpm.get_np(PFIELD):-1:1
        P = vpm.get_particle(PFIELD, i)

        if P.Gamma[1]*P.Gamma[1] + P.Gamma[2]*P.Gamma[2] + P.Gamma[3]*P.Gamma[3] < crit
            vpm.remove_particle(PFIELD, i)
        end
    end

    return false
end

function remove_particles_sphere(sim, PFIELD, T, DT; Rsphere2=(1.25*5.86)^2, Xoff=[4.0, 0, 0])

    Xvehicle = sim.vehicle.system.O

    for i in vpm.get_np(PFIELD):-1:1
        P = vpm.get_particle(PFIELD, i)
        X1 = P.X[1] - (Xvehicle[1] + Xoff[1])
        X2 = P.X[2] - (Xvehicle[2] + Xoff[2])
        X3 = P.X[3] - (Xvehicle[3] + Xoff[3])

        if X1*X1 + X2*X2 + X3*X3 > Rsphere2
            vpm.remove_particle(PFIELD, i)
        end
    end

    return false
end


text:
using Revise
using ADCME
using NNFEM
using JLD2
using PyPlot

include("nnutil.jl")
nntype= "ae_scaled"
stress_scale = 100.0

loss = constant(0.0)
for i = [1,2,4,5]
    @load "Data/domain$i.jld2" domain
    X, Y = prepare_strain_stress_data1D(domain)
    y = squeeze(nn(constant(X[:,1]), constant(X[:,2]), constant(X[:,3])))
    global loss += sum((y-Y)^2)
end
sess = Session(); init(sess)
@show run(sess,loss)
BFGS!(sess, loss, 2000)
ADCME.save(sess, "Data/learned_nn.mat")
# error()

@load "Data/domain3.jld2" domain
X, Y = prepare_strain_stress_data1D(domain)
x = constant(X)
y = squeeze(nn(constant(X[:,1]), constant(X[:,2]), constant(X[:,3])))
sess = Session(); init(sess)
close("all")
ADCME.load(sess, "Data/learned_nn.mat")
out = run(sess, y)
plot(X[:,1], out,"+", label="NN")
plot(X[:,1], Y, ".", label="Exact")
legend()

# @load "Data/domain.jld2" domain
# X, Y = prepare_strain_stress_data1D(domain)
# x = constant(X)
# y = squeeze(nn(constant(X[:,1]), constant(X[:,2]), constant(X[:,3])))
# sess = Session(); init(sess)
# close("all")
# ADCME.load(sess, "Data/learned_nn.mat")
# out = run(sess, y)
# plot(X[:,1], out,"+", label="NN")
# plot(X[:,1], Y, ".", label="Exact")
# ADCME.load(sess, "Data/trained_nn_fem.mat")
# out = run(sess, y)
# plot(X[:,1], out,">", label="End-to-end")
# legend()


text:
#!/usr/bin/julia

# pangenome_make_table.jl - Reads a pangenome profile and generates core, new, and shared gene hits for the pangenome
# By: Shaun Adkins (sadkins@som.umaryland.edu)

using ArgParse

#############
# Type defs #
#############

type HitsProfile
    profile::Dict
    genome_names::Vector{String}
    num_genomes::Integer
end

HitsProfile() = HitsProfile(Dict(), Vector{String}(), 0)

###############
# HitsProfile #
###############

function read_in_profile(hp::HitsProfile, profile)
    # Read in tab-delimited file, and grab genome names
    (data, headers) = readdlm(profile, '\t', header=true, comments=false)
    hp.num_genomes = length(headers) - 2
    hp.genome_names = headers[3:end]
    hp.profile = parse_data(hp.genome_names, data)
end

function parse_data(headers::Vector{String}, data)
    # Iterate through each line of data
    profile = Dict{String, Any}()
    for i = 1:size(data,1)
        genome = data[i,1]
        gene = data[i,2]
        # Hits to each subject genome
        hits = data[i,3:end]
        g = Dict{String, Bool}()
        for j = 1:length(hits)
            if hits[j] == 1
                g[headers[j]] = true
            end
        end

        # Create typed dictionary for genome if one does not exist
        try
			profile[genome][gene] = g
		catch KeyError e
			profile[genome] = Dict{String, Any}(gene => g)
		end
    end
    return profile
end

########
# Main #
########

function determine_sampling_rate( hp::HitsProfile, args)
    if args["comparisons"] > 0
        return estimate_multiplicity(hp.num_genomes, args["comparisons"])
    else
        return args["multiplicity"]
    end
end

function estimate_multiplicity( genomes::Integer, comparisons::Integer)
    # Determine multiplex sampling rate based on specified comparisons
    lower_mult = 0
    lower_comp = 0
    lower_theor = 0
    upper_mult = 0
    upper_comp = 0
    upper_theor = 0
    ldiff = 0
    udiff = 0

	# Bound limits for estimating multiplicity rate
	lower_bounds = 5
	upper_bounds = 5000

    # Establish multiplex bounds
    for i = lower_bounds:upper_bounds
        (est_comps, theor_comps) = estimate_comparisons(genomes, i)
        if est_comps < comparisons
            lower_mult = i
            lower_comp = est_comps
            lower_theor = theor_comps
        else
            upper_mult = i
            upper_comp = est_comps
            upper_theor = theor_comps
            break
        end
    end
    if upper_mult == lower_bounds
        (est_comps, theor_comps) = estimate_comparisons(genomes, lower_bounds)
    elseif lower_mult == upper_bounds
        return upper_bounds
    else
        lower_diff = comparisons - lower_comp
        upper_diff = upper_comp - comparisons
        if ldiff <= udiff
            return lower_mult
        else
            return upper_mult
        end
    end
end

function estimate_comparisons(genomes::Integer, mult::Integer)
    # Estimate the number of comparisons necessary for good bootstrapping
    total_comp = 0
    theor_comp = 0
    for i = range(2,genomes+1)
        theor = get_theoretical_comps(genomes, i)
        real = mult * genomes
        theor_comp += theor
        total_comp += min(theor, real)
    end
    return (total_comp, theor_comp)
end

function get_theoretical_comps(genomes::Integer, i::Integer)
    # Get number of combinations for the given size pangenome
    factorial(genomes) / (factorial(genomes - i) * factorial(i-1))
end

function create_reference_string(reference_set)
    # Create a string from a random set of genomes
    seen_vec = [ reference_set[i] for i = 1:length(reference_set) ]
    sort!(seen_vec)
    join(seen_vec, ":")
end

function categorize_shared_gene(genes_by_category, comp_genome::String, gene::String)
    # Increment shared count for this gene if we have a hit
    get!(genes_by_category[comp_genome]["shared"], gene, 0)
    genes_by_category[comp_genome]["shared"][gene] += 1
    return genes_by_category
end

function categorize_core_gene(genes_by_category, comp_genome::String, gene::String)
    # Determine if the given gene is present in all other genomes up to this point.
    genes_by_category[comp_genome]["core"][gene] = 1
    return genes_by_category
end

function categorize_new_gene(genes_by_category, comp_genome::String, gene::String)
    # Determine if the given gene is present in no other genomes up to this point.
    genes_by_category[comp_genome]["new"][gene] = 1
    return genes_by_category
end

function write_single_genomes(hp::HitsProfile, fh)
    # Write the number of genes for the individual genomes
    genome_index = Dict{String, Integer}(hp.genome_names[i] => i for i in 1:hp.num_genomes)

    # For each genome, get number of genes and write to file (also write assigned index)
    for genome in hp.genome_names
        num_genes = length(keys(hp.profile[genome]))
        index = genome_index[genome]
        write(fh, "1\t0\t0\t$num_genes\t$index\t$genome\n")
    end
end

function write_results(genes_by_category, comp_genome::String, size::Integer, hp::HitsProfile, fh)
    index = findfirst(hp.genome_names,comp_genome)
    p_size = size + 1 # Pangenome size = Ref genomes + seed genome

    # Write core/shared/new results to file
    shared_count = length(keys(genes_by_category[comp_genome]["shared"]))
    new_count = length(keys(genes_by_category[comp_genome]["new"]))
    core_count = length(keys(genes_by_category[comp_genome]["core"]))

    write(fh, "$p_size\t$core_count\t$shared_count\t$new_count\t$index\t$comp_genome\n")
end

function do_analysis_with_sampling(args, hp::HitsProfile, fh)
    multiplex = determine_sampling_rate(hp, args)

    # Iterate through each possible pangenome size
    for i = 1:(hp.num_genomes-1)
        tic()
        print("Running analysis with $i genomes...")
        sampling_counter = Dict{String, Integer}()

        fac_numer = factorial(big(hp.num_genomes-1))    # Must convert to BigInt first
        fac_denom1 = factorial(big((i+1)-1))            # ... otherwise OverflowError
        fac_denom2 = factorial(big(hp.num_genomes-(i+1)))
        true_max = fac_numer / (fac_denom1 * fac_denom2)

        # Iterate through the all genomes in set
        for j = 1:hp.num_genomes
            seen_strings = Dict{String, Bool}()
            seed_genome = hp.genome_names[j]
            # Create group of reference genomes minus the seed genome
            ref_genomes = copy(hp.genome_names) # copying the variable creates a pointer
            splice!(ref_genomes, j)
            point_count = 0
            while point_count < min(multiplex, true_max)
                # Create of string of N-reference genomes for an N+1 sized pangenome
                reference_set = rand(ref_genomes, i)
                sort_string = create_reference_string(reference_set)

                # Proceed if this is a new string
                if haskey(seen_strings, sort_string) == false
                    get!(sampling_counter, seed_genome, 1)
                    genes_by_category = Dict()
                    genes_by_category[seed_genome] = Dict("shared" => Dict(), "new" => Dict(), "core" => Dict())

                    # Iterate through all the genes in the seed genome
                    for gene in keys(hp.profile[seed_genome])
                        for ref_g in reference_set
                            # If reference genome is a key in the HitsProfile, add to "shared" genes count
                            if haskey(hp.profile[seed_genome][gene], ref_g)
                                genes_by_category = categorize_shared_gene(genes_by_category, seed_genome, gene)
                            end
                        end
                        # Determine if gene is shared by all or none of the referece genomes
                        count = get(genes_by_category[seed_genome]["shared"], gene, 0)
                        if count == length(reference_set)
                            genes_by_category = categorize_core_gene(genes_by_category, seed_genome, gene)
                        elseif count == 0
                            genes_by_category = categorize_new_gene(genes_by_category, seed_genome, gene)
                        end
                    end
                    write_results(genes_by_category, seed_genome, i, hp, fh)
                    sampling_counter[seed_genome] += 1
                    seen_strings[sort_string] = true
                    point_count += 1
                end
            end
        end

        elapsed = toq()
        println("time elapsed: $elapsed secs")
    end

end

function do_analysis_without_sampling()
    println("Doing analysis without sampling is not implemented at this time")
end

function parse_commandline()
    s = ArgParseSettings(description = "Reads a pangenome profile matrix and generates total genes in the pangenome",
    prog = "pangenome_make_pangenome.jl",
	version = "1.0",
	add_version = true,
	add_help = true)

    # The macro to add a table of arguments and options to the ArgParseSettings object
    @add_arg_table s begin
        "--profile", "-p"
            help = "Path to profile matrix"
            metavar = "/path/to/pangenome/profile.txt"
            required = true
        "--output_dir", "-o"
            help = "Directory to write the output"
            metavar = "/path/to/pangenome/dir/"
            required = true
        "--comparisons", "-c"
            help = "The number of comparisons to make for any 1 value of N (sampling)"
            metavar = "100000"
			default = 0
			range_tester = (x->x>=0)
            arg_type = Int
        "--multiplicity", "-m"
            help = "Another option for sampling based on a multiplicity factor ((sum(m*n) for n=[2..n]) number of comparisons)"
            metavar = "20"
			default = 0
			range_tester = (x->x>=0)
            arg_type = Int
    # Will not add log_file or debug options for now
    end

    # Converts the ArgParseSettings object into key/value pairs
    return parse_args(s)
end

function main()
    args = parse_commandline()
    println("Parsed args:")
    for (arg,val) in args
        println("  $arg  =>  $val")
    end

    hp = HitsProfile()
    read_in_profile(hp, args["profile"])

    output_dir = args["output_dir"]
    output_file = "$output_dir/core_shared_new.tsv"
    out_fh = open(output_file, "w")

    ## our output table will take the form:
    ## genome_count  core    shared    new  root_genome
    write_single_genomes(hp, out_fh)

    if args["comparisons"] == 0 && args["multiplicity"] == 0
        do_analysis_without_sampling()
    else
        do_analysis_with_sampling(args, hp, out_fh)
    end

    close(out_fh)
end

main()
quit()

max_stars_repo_name: IGS/ergatis

max_issues_repo_name: IGS/ergatis

max_forks_repo_name: IGS/ergatis


text:
########################################################################
# Expected item info -  for 1PL pars and 1-D latent                    #
########################################################################
"""
```julia
_item_expected_information(parameters::Parameters1PL, latent::Latent1D)
```

# Description

It computes the expected information (-second derivative of the likelihood) with respect to the difficulty parameter of the 1PL model.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`parameters::Parameters1PL`** : Required. A 1-parameter logistic parameters object.
- **`latent::Latent1D`** : Required. A 1-dimensional `Latent1D` latent variable.

# Output

A `Float64` scalar.
"""
function _item_expected_information(parameters::Parameters1PL, latent::Latent1D)
    p = _probability(latent, parameters)
    return p * (1 - p)::Float64
end

########################################################################
# Expected info -  2PL pars and 1-D latent                             #
########################################################################
"""
```julia
_item_expected_information(parameters::Parameters2PL, latent::Latent1D)
```

# Description

It computes the expected information (-second derivative of the likelihood) with respect to the 2 parameters of the 2PL model.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`parameters::Parameters1PL`** : Required. A 2-parameter logistic parameters object.
- **`latent::Latent1D`** : Required. A 1-dimensional `Latent1D` latent variable.

# Output

A ``2 \times 2`` matrix of the expected informations.
"""
function _item_expected_information(parameters::Parameters2PL, latent::Latent1D)
    p = _probability(latent, parameters)
    p_1_p = (1 - p) * p
    i_aa = p_1_p * (latent.val - parameters.b)^2
    i_ab = -parameters.a * p_1_p * (latent.val - parameters.b)
    i_bb = parameters.a^2 * p_1_p
    return [i_aa i_ab; i_ab i_bb]::Matrix{Float64}
end

########################################################################
# Expected info -  for 3PL pars and 1-D latent                         #
########################################################################
"""
```julia
_item_expected_information(parameters::Parameters3PL, latent::Latent1D)
```

# Description

It computes the expected information (-second derivative of the likelihood) with respect to the 3 parameters of the 3PL model.
It follows the parametrization ``a(θ - b)``.
It should be always positive.

# Arguments
- **`parameters::Parameters1PL`** : Required. A 3-parameter logistic parameters object.
- **`latent::Latent1D`** : Required. A 1-dimensional `Latent1D` latent variable.

# Output
A ``3 \times 3`` matrix of the expected informations.
"""
function _item_expected_information(parameters::Parameters3PL, latent::Latent1D)
    p = _probability(latent, parameters)
    den = p * (1 - parameters.c)^2
    theta_b = latent.val - parameters.b
    p_c = p - parameters.c
    i_cc = (1 - p) / den # vv p. 51 Kim, Baker, 2008 : - E(L_33)
    i_aa = i_cc * p_c^2 * theta_b^2  # v p. 51 Kim, Baker, 2008
    i_ab = - i_aa / theta_b * parameters.a
    i_ac = i_cc * p_c * theta_b
    i_bc = i_cc * p_c * parameters.a
    i_bb = (i_bc ^ 2) / i_cc #v p. 51 Kim, Baker, 2008
    return [i_aa i_ab i_ac; i_ab i_bb i_bc; i_ac i_bc i_cc]::Matrix{Float64}
end


text:
using CodecBase
using Test
import TranscodingStreams:
    TranscodingStream,
    test_roundtrip_read,
    test_roundtrip_write,
    test_roundtrip_lines,
    test_roundtrip_transcode

@testset "DecodeError" begin
    error = CodecBase.DecodeError("basexx: invalid data")
    @test error isa CodecBase.DecodeError
    @test error isa Exception
    @test sprint(showerror, error) == "DecodeError: basexx: invalid data"
    try
        throw(error)
        @test false
    catch ex
        @test ex isa CodecBase.DecodeError
    end
end

@testset "Base16" begin
    CodeTable16 = CodecBase.CodeTable16
    @test_throws ArgumentError CodeTable16("Sigur Rós")
    @test_throws ArgumentError CodeTable16("Takk")
    table = copy(CodecBase.BASE16_UPPER)
    @test_throws ArgumentError CodecBase.ignorechars!(table, "Jónsi")

    @test transcode(Base16Encoder(), b"") == b""
    @test transcode(Base16Encoder(), b"f") == b"66"
    @test transcode(Base16Encoder(), b"fo") == b"666F"
    @test transcode(Base16Encoder(), b"foo") == b"666F6F"
    @test transcode(Base16Encoder(), b"foob") == b"666F6F62"
    @test transcode(Base16Encoder(), b"fooba") == b"666F6F6261"
    @test transcode(Base16Encoder(), b"foobar") == b"666F6F626172"

    @test transcode(Base16Encoder(lowercase=false), b"自問自答") == b"E887AAE5958FE887AAE7AD94"
    @test transcode(Base16Encoder(lowercase=true),  b"自問自答") == b"e887aae5958fe887aae7ad94"

    @test transcode(Base16Decoder(), b"") == b""
    @test transcode(Base16Decoder(), b"66") == b"f"
    @test transcode(Base16Decoder(), b"666F") == b"fo"
    @test transcode(Base16Decoder(), b"666F6F") == b"foo"
    @test transcode(Base16Decoder(), b"666F6F62") == b"foob"
    @test transcode(Base16Decoder(), b"666F6F6261") == b"fooba"
    @test transcode(Base16Decoder(), b"666F6F626172") == b"foobar"
    @test transcode(Base16Decoder(), b"666f6F626172") == b"foobar"

    @test transcode(Base16Decoder(), b"   ") == b""
    @test transcode(Base16Decoder(), b"666\nF6F") == b"foo"
    @test transcode(Base16Decoder(), b"666\r\nF6F") == b"foo"
    @test transcode(Base16Decoder(), b"6  66\r\nF6F   ") == b"foo"
    @test transcode(Base16Decoder(), b"  66\t6F\t6F\n") == b"foo"

    DecodeError = CodecBase.DecodeError
    @test_throws DecodeError transcode(Base16Decoder(), b"a")
    @test_throws DecodeError transcode(Base16Decoder(), b"aaa")
    @test_throws DecodeError transcode(Base16Decoder(), b"aaaaa")
    @test_throws DecodeError transcode(Base16Decoder(), b"\0")
    @test_throws DecodeError transcode(Base16Decoder(), b"a\0")
    @test_throws DecodeError transcode(Base16Decoder(), b"aa\0")
    @test_throws DecodeError transcode(Base16Decoder(), b"aaa\0")

    test_roundtrip_read(Base16EncoderStream, Base16DecoderStream)
    test_roundtrip_write(Base16EncoderStream, Base16DecoderStream)
    test_roundtrip_lines(Base16EncoderStream, Base16DecoderStream)
    test_roundtrip_transcode(Base16Encoder, Base16Decoder)
end

@testset "Base32" begin
    CodeTable32 = CodecBase.CodeTable32
    @test_throws ArgumentError CodeTable32("Sigur Rós", '=')
    @test_throws ArgumentError CodeTable32("Takk", '=')

    # standard
    @test transcode(Base32Encoder(), b"") == b""
    @test transcode(Base32Encoder(), b"f") == b"MY======"
    @test transcode(Base32Encoder(), b"fo") == b"MZXQ===="
    @test transcode(Base32Encoder(), b"foo") == b"MZXW6==="
    @test transcode(Base32Encoder(), b"foob") == b"MZXW6YQ="
    @test transcode(Base32Encoder(), b"fooba") == b"MZXW6YTB"
    @test transcode(Base32Encoder(), b"foobar") == b"MZXW6YTBOI======"

    @test transcode(Base32Decoder(), b"") == b""
    @test transcode(Base32Decoder(), b"MY======") == b"f"
    @test transcode(Base32Decoder(), b"MZXQ====") == b"fo"
    @test transcode(Base32Decoder(), b"MZXW6===") == b"foo"
    @test transcode(Base32Decoder(), b"MZXW6YQ=") == b"foob"
    @test transcode(Base32Decoder(), b"MZXW6YTB") == b"fooba"
    @test transcode(Base32Decoder(), b"MZXW6YTBOI======") == b"foobar"

    DecodeError = CodecBase.DecodeError
    @test_throws DecodeError transcode(Base32Decoder(), b"MZXW6=")
    @test_throws DecodeError transcode(Base32Decoder(), b"MZXW6==")
    @test_throws DecodeError transcode(Base32Decoder(), b"MZX\0W6===")
    @test_throws DecodeError transcode(Base32Decoder(), b"MZXW6===\0")

    # extended hex
    @test transcode(Base32Encoder(hex=true), b"") == b""
    @test transcode(Base32Encoder(hex=true), b"f") == b"CO======"
    @test transcode(Base32Encoder(hex=true), b"fo") == b"CPNG===="
    @test transcode(Base32Encoder(hex=true), b"foo") == b"CPNMU==="
    @test transcode(Base32Encoder(hex=true), b"foob") == b"CPNMUOG="
    @test transcode(Base32Encoder(hex=true), b"fooba") == b"CPNMUOJ1"
    @test transcode(Base32Encoder(hex=true), b"foobar") == b"CPNMUOJ1E8======"

    @test transcode(Base32Decoder(hex=true), b"") == b""
    @test transcode(Base32Decoder(hex=true), b"CO======") == b"f"
    @test transcode(Base32Decoder(hex=true), b"CPNG====") == b"fo"
    @test transcode(Base32Decoder(hex=true), b"CPNMU===") == b"foo"
    @test transcode(Base32Decoder(hex=true), b"CPNMUOG=") == b"foob"
    @test transcode(Base32Decoder(hex=true), b"CPNMUOJ1") == b"fooba"
    @test transcode(Base32Decoder(hex=true), b"CPNMUOJ1E8======") == b"foobar"

    test_roundtrip_read(Base32EncoderStream, Base32DecoderStream)
    test_roundtrip_write(Base32EncoderStream, Base32DecoderStream)
    test_roundtrip_lines(Base32EncoderStream, Base32DecoderStream)
    test_roundtrip_transcode(Base32Encoder, Base32Decoder)
end

@testset "Base64" begin
    CodeTable64 = CodecBase.CodeTable64
    @test_throws ArgumentError CodeTable64("Sigur Rós", '=')
    @test_throws ArgumentError CodeTable64("Takk", '=')

    @test transcode(Base64Encoder(), b"") == b""
    @test transcode(Base64Encoder(), b"f") == b"Zg=="
    @test transcode(Base64Encoder(), b"fo") == b"Zm8="
    @test transcode(Base64Encoder(), b"foo") == b"Zm9v"
    @test transcode(Base64Encoder(), b"foob") == b"Zm9vYg=="
    @test transcode(Base64Encoder(), b"fooba") == b"Zm9vYmE="
    @test transcode(Base64Encoder(), b"foobar") == b"Zm9vYmFy"

    @test transcode(Base64Encoder(urlsafe=false), b"響き") == b"6Z+/44GN"
    @test transcode(Base64Encoder(urlsafe=true),  b"響き") == b"6Z-_44GN"

    @test transcode(Base64Decoder(), b"") == b""
    @test transcode(Base64Decoder(), b"Zg==") == b"f"
    @test transcode(Base64Decoder(), b"Zm8=") == b"fo"
    @test transcode(Base64Decoder(), b"Zm9v") == b"foo"
    @test transcode(Base64Decoder(), b"Zm9vYg==") == b"foob"
    @test transcode(Base64Decoder(), b"Zm9vYmE=") == b"fooba"
    @test transcode(Base64Decoder(), b"Zm9vYmFy") == b"foobar"

    @test transcode(Base64Decoder(), b"Zg=\n=") == b"f"
    @test transcode(Base64Decoder(), b"Zm9v\nYmFy") == b"foobar"
    @test transcode(Base64Decoder(), b"Zg=\r\n=") == b"f"
    @test transcode(Base64Decoder(), b"Z m  9\rv\r\nYm\nF\t\ty") == b"foobar"
    @test transcode(Base64Decoder(), b"  Zm9vYmFy  ") == b"foobar"
    @test transcode(Base64Decoder(), b"      Zm     9v      Ym    Fy      ") == b"foobar"

    @test transcode(Base64Decoder(urlsafe=false), b"6Z+/44GN") == b"響き"
    @test transcode(Base64Decoder(urlsafe=true),  b"6Z-_44GN") == b"響き"

    DecodeError = CodecBase.DecodeError
    @test_throws DecodeError transcode(Base64Decoder(), b"a")
    @test_throws DecodeError transcode(Base64Decoder(), b"aa")
    @test_throws DecodeError transcode(Base64Decoder(), b"aaa")
    @test_throws DecodeError transcode(Base64Decoder(), b"aaaaa")
    @test_throws DecodeError transcode(Base64Decoder(), b"\0")
    @test_throws DecodeError transcode(Base64Decoder(), b"a\0")
    @test_throws DecodeError transcode(Base64Decoder(), b"aa\0")
    @test_throws DecodeError transcode(Base64Decoder(), b"aaa\0")

    test_roundtrip_read(Base64EncoderStream, Base64DecoderStream)
    test_roundtrip_write(Base64EncoderStream, Base64DecoderStream)
    test_roundtrip_lines(Base64EncoderStream, Base64DecoderStream)
    test_roundtrip_transcode(Base64Encoder, Base64Decoder)
end


text:
module DSL

export create_window, show_window, setup_grid, swap_window_buffers, save_framebuffer, isclosed
export create_chart_2d, create_chart_3d
# chart
export set_axes_limits, set_axes_titles
export create_image
export create_plot, create_surface, create_vector_field
export set_color
export set_legend, set_marker_size
export draw_chart
export update_vertex_buffer, update_direction_buffer

import ..Forge
using ColorTypes

struct Window
    ref::Ref{Forge.fg_window}
end

struct Image
    ref::Ref{Forge.fg_image}
end

struct Chart
    ref::Ref{Forge.fg_chart}
end

struct Plot
    ref::Ref{Forge.fg_plot}
end

struct Surface
    ref::Ref{Forge.fg_surface}
end

struct VectorField
    ref::Ref{Forge.fg_vector_field}
end

struct Histogram
    ref::Ref{Forge.fg_histogram}
end

# window
function create_window(width::Real, height::Real, title::String, shared_window=C_NULL, visible::Bool=true)::Window
    refwindow = Ref{Forge.fg_window}()
    Forge.fg_create_window(refwindow, Cint(width), Cint(height), title, shared_window, visible)
    Window(refwindow)
end

function show_window(window::Window)
    Forge.fg_show_window(window.ref[])
end

function setup_grid(window::Window, rows::Real, cols::Real)
    Forge.fg_setup_window_grid(rows, cols, window.ref[])
end

function create_image(width::Real, height::Real, format::Forge.fg_channel_format=Forge.FG_RGBA, dtype::Forge.fg_dtype=Forge.u8)::Image
    refimage = Ref{Forge.fg_image}()
    Forge.fg_create_image(refimage[], width, height, format, dtype)
    Image(refimage)
end

function swap_window_buffers(window::Window)
    Forge.fg_swap_window_buffers(window.ref[])
end

function save_framebuffer(window::Window, path::String)
    Forge.fg_save_window_framebuffer(path, window.ref[])
end

function isclosed(window::Window)::Bool
    refisclosed = Ref{Bool}()
    Forge.fg_close_window(refisclosed, window.ref[])
    refisclosed[]
end

# chart
function create_chart_2d()::Chart
    refchart = Ref{Forge.fg_chart}()
    Forge.fg_create_chart(refchart, Forge.FG_CHART_2D)
    Chart(refchart)
end

function create_chart_3d()::Chart
    refchart = Ref{Forge.fg_chart}()
    Forge.fg_create_chart(refchart, Forge.FG_CHART_3D)
    Chart(refchart)
end

function set_axes_limits(chart::Chart, xmin::Real, xmax::Real, ymin::Real, ymax::Real, zmin::Real=0, zmax::Real=0)
    Forge.fg_set_chart_axes_limits(chart.ref[], Cfloat.([xmin, xmax, ymin, ymax, zmin, zmax])...)
end

function set_axes_titles(chart::Chart, x::String, y::String, z::String="")
    Forge.fg_set_chart_axes_titles(chart.ref[], x, y, z)
end

function create_surface(chart::Chart, xsize::Real, ysize::Real, dtype::UInt32=Forge.FG_FLOAT32, plot_type::Forge.fg_plot_type=Forge.FG_PLOT_LINE, marker_type::Forge.fg_marker_type=Forge.FG_MARKER_NONE)::Surface
    refsurface = Ref{Forge.fg_surface}()
    Forge.fg_add_surface_to_chart(refsurface, chart.ref[], Cuint(xsize), Cuint(ysize), dtype, plot_type, marker_type)
    Surface(refsurface)
end

function create_plot(chart::Chart, npoints::Real)::Plot
    create_plot(chart, npoints, Forge.FG_FLOAT32, Forge.FG_PLOT_LINE, Forge.FG_MARKER_NONE)
end
function create_plot(chart::Chart, npoints::Real, plot_type::Forge.fg_plot_type, marker_type::Forge.fg_marker_type)::Plot
    create_plot(chart, npoints, Forge.FG_FLOAT32, plot_type, marker_type)
end
function create_plot(chart::Chart, npoints::Real, dtype::UInt32, plot_type::Forge.fg_plot_type, marker_type::Forge.fg_marker_type)::Plot
    refplot = Ref{Forge.fg_plot}()
    Forge.fg_add_plot_to_chart(refplot, chart.ref[], Cuint(npoints), dtype, plot_type, marker_type)
    Plot(refplot)
end

function create_vector_field(chart::Chart, npoints::Real, dtype::UInt32=Forge.FG_FLOAT32)::VectorField
    reffield = Ref{Forge.fg_vector_field}()
    Forge.fg_add_vector_field_to_chart(reffield, chart.ref[], Cuint(npoints), dtype)
    VectorField(reffield)
end

function create_histogram(chart::Chart, nbins::Real, dtype::Forge.fg_dtype=Forge.s32)::Histogram
    refhisto = Ref{Forge.fg_histogram}()
    Forge.fg_add_histogram_to_chart(refhisto, nbinsp, dtype)
    Histogram(refhisto)
end

function draw_chart(f::Function, window::Window, chart::Chart)
    while f(window)
        Forge.fg_draw_chart(window.ref[], chart.ref[])
    end
end

# set_color
function (::Type{RGBA})(color::Unsigned)
    RGBA(map(n -> (color >> n) & 0xFF / 255, (24, 16, 8, 0))...)
end

set_color(plot::Plot, color::Unsigned) = set_color(plot, RGBA(color))
set_color(plot::Plot, color::RGB) = set_color(plot, convert(RGBA, color))
function set_color(plot::Plot, color::RGBA)
    Forge.fg_set_plot_color(plot.ref[], Cfloat.([color.r, color.g, color.b, color.alpha])...)
end

set_color(field::VectorField, color::Unsigned) = set_color(field, RGBA(color))
set_color(field::VectorField, color::RGB) = set_color(field, convert(RGBA, color))
function set_color(field::VectorField, color::RGBA)
    Forge.fg_set_vector_field_color(field.ref[], Cfloat.([color.r, color.g, color.b, color.alpha])...)
end

set_color(histogram::Histogram, color::Unsigned) = set_color(histogram, RGBA(color))
set_color(histogram::Histogram, color::RGB) = set_color(histogram, convert(RGBA, color))
function set_color(histogram::Histogram, color::RGBA)
    Forge.fg_set_histogram_color(histogram.ref[], Cfloat.([color.r, color.g, color.b, color.alpha])...)
end

function set_legend(plot::Plot, legend::String)
    Forge.fg_set_plot_legend(plot.ref[], legend)
end

function set_marker_size(plot::Plot, marker_size::Real)
    Forge.fg_set_plot_marker_size(plot.ref[], Cfloat(marker_size))
end

# update vertex buffer
function update_vertex_buffer(plot::Plot, data::Array)
    refvertices = Ref{Cuint}()
    refvertexbuffersize = Ref{Cuint}()
    Forge.fg_get_plot_vertex_buffer(refvertices, plot.ref[])
    Forge.fg_get_plot_vertex_buffer_size(refvertexbuffersize, plot.ref[])
    Forge.fg_update_vertex_buffer(refvertices[], refvertexbuffersize[], Array{Cfloat}(data))
end

function update_vertex_buffer(surface::Surface, data::Array)
    refvertices = Ref{Cuint}()
    refvertexbuffersize = Ref{Cuint}()
    Forge.fg_get_surface_vertex_buffer(refvertices, surface.ref[])
    Forge.fg_get_surface_vertex_buffer_size(refvertexbuffersize, surface.ref[])
    Forge.fg_update_vertex_buffer(refvertices[], refvertexbuffersize[], Array{Cfloat}(data))
end

function update_vertex_buffer(field::VectorField, data::Array)
    refvertices = Ref{Cuint}()
    refvertexbuffersize = Ref{Cuint}()
    Forge.fg_get_vector_field_vertex_buffer(refvertices, field.ref[])
    Forge.fg_get_vector_field_vertex_buffer_size(refvertexbuffersize, field.ref[])
    Forge.fg_update_vertex_buffer(refvertices[], refvertexbuffersize[], Array{Cfloat}(data))
end

function update_direction_buffer(field::VectorField, data::Array)
    refdirections = Ref{Cuint}()
    refdirectionbuffersize = Ref{Cuint}()
    Forge.fg_get_vector_field_direction_buffer(refdirections, field.ref[])
    Forge.fg_get_vector_field_direction_buffer_size(refdirectionbuffersize, field.ref[])
    Forge.fg_update_vertex_buffer(refdirections[], refdirectionbuffersize[], Array{Cfloat}(data))
end

function update_pixel_buffer(image::Image, data::Array)
    bufferid
    buffersize
    Forge.fg_update_pixel_buffer(bufferid, buffersize, Array{Cfloat}(data))
end

function update_vertex_buffer(histogram::Histogram, data::Array)
    refvertices = Ref{Cuint}()
    refvertexbuffersize = Ref{Cuint}()
    Forge.fg_get_histogram_vertex_buffer(refvertices, field.ref[])
    Forge.fg_get_histogram_vertex_buffer_size(refvertexbuffersize, field.ref[])
    Forge.fg_update_vertex_buffer(refvertices[], refvertexbuffersize[], Array{Cfloat}(data))
end

function update_color_buffer(histogram::Histogram, data::Array)
    refcolorbuffer = Ref{Cuint}()
    refcolorbuffersize = Ref{Cuint}()
    Forge.fg_get_histogram_color_buffer(refcolorbuffer, histogram.ref[])
    Forge.fg_get_histogram_color_buffer_size(refcolorbuffersize, histogram.ref[])
    Forge.fg_update_vertex_buffer(refcolorbuffer[], refcolorbuffersize[], Array{Cfloat}(data))
end

end # module DSL


text:
using AdFem

for degree in [1, 2]
    m = 100; n = 100; h = 1/n
    mesh = Mesh(m, n, h, degree = degree)

    function f(x, y)
        2 * x^2 + y
    end

    f_fem = eval_f_on_fem_pts(f, mesh)
    f_gauss = fem_to_gauss_points(constant(f_fem), mesh)
    sess = Session(); init(sess)
    f_gauss = run(sess, f_gauss)

    f_gauss2 = eval_f_on_gauss_pts(f, mesh)

    @show maximum(abs.(f_gauss .- f_gauss2))
end