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

using Markdown
using InteractiveUtils

# ╔═╡ f39f9c31-e702-4de9-81c8-07e43584ec96
# Files & Data processing
using CSV, DataFrames, Glob

# ╔═╡ 579cbf1f-be0a-4c80-b4a8-57e7a92a0fc7
begin
	# Plotting/Visualization
	using Plots
	theme(:ggplot2)
	gr()
end

# ╔═╡ 473dedb8-da32-4687-9b4c-ffb4976f2d54
# Math & Stats
using LinearAlgebra, Distributions

# ╔═╡ 330163b0-ad49-11ec-3128-afba43b7bdd0
md"""
# 2IFC Data Analysis

Goal: Find a range of participant σ's that best explain performance on a 2-Interval Forced-Choice Task
"""

# ╔═╡ 940b5333-3a5b-4fef-9c0e-356fd4f08354
md"""
# Data Preprocessing
"""

# ╔═╡ 5fa46984-a715-4c29-9373-5af0b62f0c70
md"""
Data Format & Legend

Condition A: "Block Order"
* 1 = Sides-Dots-Both
* 2 = Dots-Sides-Both

Condition B: "In BOTH block, which was asked first?"
* 1 = Sides asked first
* 2 = Dots asked first

Sides: "add 6 to variable"
* 0=6
* 1=7
* 2=8
* 3=9
* 4=10

Spacing: "add 4 to variable"
* 0=4
* 1=5
* 2=6
* 3=7
* 4=8

Interval Designation:
* Reference = chosen objects that aren't the middle one
* Comparison = Middle object (8 sides, 6mm dot spacing)
"""

# ╔═╡ d01b1239-57d4-4e64-b1c1-7b112387d6a1
md"""
## Load CSV Data as DataFrame

For each folder in 2IFC data, convert block csv files to DataFrames
"""

# ╔═╡ d3b0223d-5946-4480-a836-ffbc54c035d6
readdir("./2ifcdata")

# ╔═╡ 14bd5283-578a-4a08-872e-6e73abf9e2c5
function get_csv(dir_with_csvs::String)
    ispath(dir_with_csvs) || error("Could not find path: $dir_with_csvs")
    csv_files = glob("*csv", dir_with_csvs)
    return csv_files
end

# ╔═╡ 179c2ed9-664b-4f60-a5d8-e65cd7890062
md"""
Directory structure:

	2ifc_data
	|_ p00n
		|- b1.csv
		|- ...
		|- bk.csv

"""

# ╔═╡ bc5cbbe7-2f9d-45f8-bf65-830468516292
"""
	get_2ifc_data(ifc_dir::String)

For each participant directory ("p00n"), load csv block data into dataframe
"""
function get_2ifc_data(ifc_dir::String)
    ispath(ifc_dir) || error("Could not find path: $ifc_dir")
	# p = participant
	p_dirs = readdir(ifc_dir, join=true)
	p_data = [CSV.read(get_csv(p_dir), DataFrame) for p_dir in p_dirs]
	return p_data
end

# ╔═╡ 44cde6d2-584a-4d33-99ca-b056d8359225
# p = participant
p_dfs_init = get_2ifc_data("./2ifcdata")

# ╔═╡ 9646aea5-73c6-4eda-a9b0-e329f5a1cc9e
function filter_blocks(p_df, blocks)
	return filter(:Block => block_num -> block_num in blocks, p_df)
end

# ╔═╡ 75382891-f42b-4a5b-a49f-ca145150bd85
p_dfs = [filter_blocks(p_df, [3.0, 6.0]) for p_df in p_dfs_init]

# ╔═╡ 8a4433de-25cc-4eb9-9252-96f3f1e1917d
md"""
## Separate Feature Data

For each participant dataframe, separate the sides, dots, heavy sides, and heavy dots data
"""

# ╔═╡ 3558d50a-4972-4483-9e02-65f891bf1d2c
feats = (
	sides="Sides",
	dots="Spacing",
	sides_correct="Sides_Correct",
	dots_correct="Spacing_Correct",
	sides_first="Sides_First?",
	dots_first="Spacing_First?",
)

# ╔═╡ 0dd8340c-00fb-4e06-9b8d-ec30f8ee753e
"""
	get_feats_data(df, feats)

Separates sides, dots, heavy sides, and heavy dots data for given participant dataframe

Assumptions:
* `0.0` value in `dots_first` column means participant was not asked about dots (similar idea with `sides_first`)
"""
function get_feats_data(df, feats)
	sides_cols = [feats.sides, feats.sides_correct]
	dots_cols = [feats.dots, feats.dots_correct]

	sides_index = df[!, feats.dots_first] .== 0.0
	sides = df[sides_index, sides_cols]

	dots_index = df[!, feats.sides_first] .== 0.0
	dots = df[dots_index, dots_cols]

	# Heavy condition
	sides_dots_index = map(!, sides_index) .& map(!, dots_index)
	sides_dots = df[sides_dots_index, vcat(sides_cols, dots_cols)]

	heavy_sides = sides_dots[!, sides_cols]
	heavy_dots = sides_dots[!, dots_cols]

	nrow(sides) + nrow(dots) + nrow(sides_dots) == nrow(df) || error("DataFrames should be mutually exclusive")
	return (;sides, dots, heavy_sides, heavy_dots)
end

# ╔═╡ 01880206-93c4-49d6-a7d6-0501c568d7b8
"""
	get_2ifc_feats(p_dfs, feats)

Separates feature data (sides, dots, heavy sides, heavy dots) for each participant
"""
function get_2ifc_feats(p_dfs, feats)
	return [get_feats_data(p_df, feats) for p_df in p_dfs]
end

# ╔═╡ 257fb4e8-2b0a-4e73-9dad-cb778d144d40
p_feats = get_2ifc_feats(p_dfs, feats)

# ╔═╡ 3eaf9582-9392-4970-bed2-455cb5756fb0
md"""
### Calculate Proportion of Correct Results

For each set of participant responses, calculate the porportion of times the participant correctly said "Reference > Comparison" for given block type (sides, dots, or both)
"""

# ╔═╡ 03d82f2a-359f-4adb-8390-27bd144c71b7
begin
	"""
		get_feats_correct(df::DataFrame; MIDDLE, categories)

	Return proportion of correct participant responses for given feature
	"""
	function get_feats_correct(df::DataFrame; MIDDLE = 2, categories = (0.0, 1.0, 3.0, 4.0))
		feat, feat_correct = names(df)
	    feat_p_correct = Vector{Pair}(undef, length(categories))
		for (i, category) in enumerate(categories)
			feat_responses = df[df[!, feat] .== category, feat_correct]
			p_correct = sum(feat_responses) / length(feat_responses)
	    	# Does Reference object have more [feature] than Comparison object
			feat_p_correct[i] = category => (category < MIDDLE) ? 1.0 - p_correct : p_correct
	    end
	    return feat_p_correct
	end

	# Example input: Participant 1 sides
	get_feats_correct(p_feats[1].sides)
end

# ╔═╡ d436b0ea-fb7c-4e6f-98ff-3c0840099e4d
"""
	get_2ifc_correct(p_feats)

Applies `get_feats_correct()` for each feature, for each participant
"""
function get_2ifc_correct(p_feats)
	# For all participants
	p_feat_correct = Vector{NamedTuple}(undef, length(p_feats))
	for (i, p_feat) in enumerate(p_feats)
		# For each feature
		feats_correct = (get_feats_correct(df) for df in p_feat)
		p_feat_correct[i] = (; zip(keys(p_feat), feats_correct)...)
	end
	return p_feat_correct
end

# ╔═╡ 0c839d64-d766-45eb-95a5-db8e27eeb50b
p_feats_correct = get_2ifc_correct(p_feats)

# ╔═╡ dc3c304c-be64-4b45-bab8-f7079b699683
md"""
## Plot Participant Data
"""

# ╔═╡ e67ca332-e489-4f74-b6e3-1db35a672e23
function plot_feat_correct(p_feat_correct, title)
	p = scatter(first.(p_feat_correct), last.(p_feat_correct),
		title = title,
		xlabel = "Category",
		ylabel = "R > C",
		ylim = [0.0, 1.0],
		legend = false,
	)
	return p
end

# ╔═╡ 7424b7b8-b6d7-4107-9e0a-04889a5a512d
begin
	function plot_participant_data(p_data::NamedTuple)
	    feat_plots = [plot_feat_correct(feat_data, feat) for (feat, feat_data) in zip(keys(p_data), p_data)]

		if length(feat_plots) > 2
	    	layout = (length(feat_plots) ÷ 2, length(feat_plots) ÷ 2)
		else
			layout = length(feat_plots)
		end

	    p = plot(feat_plots...,
			layout = layout,
			legend = false,
	        plot_title = "2IFC Performance",
	        plot_titlefontsize= 14,
	    )
		return p
	end

	# Example participant 1
	plot_participant_data(p_feats_correct[1])
end

# ╔═╡ 5479abe9-8937-4de9-a31b-afa45d1a2a7e
Markdown.parse("""
# Bayesian Parameter Estimation
""")

# ╔═╡ 05e90ddc-2ed3-4432-b3ba-da72a4326222
Markdown.parse("""

From the csv data, we calculated the proportion of times the participant said the reference object was greater than the comparison. Consider the data for the sides feature:

```math
D_\\text{sides} = [0, 0.1, 0.85, 1.0]
```

The index of each value in d corresponds with the category (0-indexed). For example, consider ``d[\\text{Category }0] = 0.0``. The participant said Category 0 objects had more sides than the Comparison object 0.0% of the time. This isn't surprising since Category 0 objects have fewer sides than the Comparison.
""")

# ╔═╡ 8f769d7c-5776-4de2-a614-35c976f3d0ff
Markdown.parse("""
Each trial is a binary decision, so we can model the experiment with a Brenoulli process. The likelihood that a participant would say "Reference > Comparison" ``s`` times, over ``n`` trials, with a probability ``p``, for ``k`` data points, is:

```math
L(D | s, n, p) = \\prod_k {n \\choose s_k} p^{s_k} (1 - p)^{n - s_k}
```

But what is ``p``?
""")

# ╔═╡ 73f8e4b2-155b-4524-8163-780c48e33135
Markdown.parse("""
``p`` depends on the category Comparison and Reference objects, and how good the participant at comparing the objects. We'll assume the participant's probability of saying "Reference > Comparison" is normally distributed. So, when the participant is presented with an object of category ``x``, the chance that their ``\\text{response} = `` "Reference > Comparison" is modeled by a cumulative normal distribution:

```math
p = P(\\text{response} | x; μ=2, σ) = \\frac{1}{σ\\sqrt{2π}} \\int_{-\\infty}^x \\exp( - \\frac{(x - μ)^2}{2σ^2})
```

This makes sense because the more sides the Reference object has, the more likely the participant will say "Reference > Comparison".
""")

# ╔═╡ 3732a3e2-7760-4abf-b217-61ec019fa27f
Markdown.parse("""

Although we don't know exactly what the participant's sigma is, we can estimate it with Bayes theorem
(We can omit the prior because it is uniform across all hypotheses)

```math
P(μ, σ | D) = \\frac{P(D | μ, σ)}{\\sum_{μ, σ} P(D | μ, σ)}
= \\frac{L(D)}{\\sum_{μ, σ} L(D)}
```
""")

# ╔═╡ eca3fc1c-1841-4252-9b40-27e7b20c991a
Markdown.parse("""
### Computational Considerations

#### Simplifying Equations

When we use Bayes theorem, the ``{n \\choose s_k}`` will factor out, so we can just use: ``L(D) = \\prod_k p^{s_k} (1 - p)^{n - s_k}``

#### Log Likelihoods

The formulas presented above tend to lead to many numerical errors, like overflows and underflows.

```math
\\log \\big( L \\big) = \\sum_k s_k \\log(p(x_k; μ, σ)) + \\sum_k (n - s_k) \\log(1 - p(x_k; μ, σ))
```

Therefore, to find the best parameters, ``\\hat{θ}``, the log form lets us use the maximum likelihood estimation technique without running overflows or underflows:

```math
\\hat{θ} = \\text{argmax}_{μ, σ} ( \\log(L) )
```
""")

# ╔═╡ 90929ce1-042f-4045-854a-ded3d05cbe2e
"""
	LL(s, x, μ, σ)

Log Likelihood function

* `s` = performance for each category
* `x` = categories
* `μ` = mean
* `σ` = noise
"""
function LL(s::AbstractVector, x::AbstractVector, μ, σ)
	return sum(s .* logcdf(Normal(μ, σ), x) .+ (1.0 .- s) .* log.(1.0 .- cdf(Normal(μ, σ), x)))
end

# ╔═╡ aec0855f-5630-4620-972e-d6bddf98f95c
"""
	calc_lls(μ_s::Vector, σ_s::Vector, p_data::Vector{Pair})

Log Likelihood function

`p_data` = participant data
"""
function calc_lls(μ_hyp::Vector, σ_hyp::Vector, feat_data)
	categories = first.(feat_data)
    performance = last.(feat_data)
    lls = [LL(performance, categories, μ, σ) for μ in μ_hyp, σ in σ_hyp]
	return lls
end

# ╔═╡ 48ab0022-3dc5-4172-a0f1-16c7b36f0ea9
begin
	function calc_posteriors(likelihoods::Matrix)
	   return likelihoods ./ sum(likelihoods)
	end

	function calc_posteriors(likelihoods::Vector, priors::Vector)
	    return (likelihoods .* priors) ./ sum(likelihoods .* priors)
	end
end

# ╔═╡ 004a7f1d-e4a3-4727-ad44-adbd7de70f10
begin
	"""
		estimate_feat_σ()

	Estimates σ for the given feature (e.g. sides)
	"""
	function estimate_feat_σ(μ_hyp::Vector, σ_hyp::Vector, feat_data)
		loglikelihoods = calc_lls(μ_hyp, σ_hyp, feat_data)

		likelihoods = exp.(loglikelihoods .- maximum(loglikelihoods))

		posteriors = calc_posteriors(likelihoods)

		best_μ_index, best_σ_index = Tuple(argmax(posteriors))
		return (best_μ=μ_hyp[best_μ_index], best_σ=σ_hyp[best_σ_index])
	end
	# Example
	estimate_feat_σ([2], collect(0.5:0.0005:5.0), p_feats_correct[1].heavy_sides)
end

# ╔═╡ f0463c17-751d-4b4f-b510-4cc4b1934524
begin
	function test_estimation()
		# Example from teams:
		ex_feat_data = [0.0 => 0.0, 1.0 => 0.2, 3.0 => 1.0, 4 => 0.8]
		# Expected sigma = 1.2495
		true_ans = 1.2495
		ex_μs = [2]
		ex_σs = collect(0.5:0.0005:5.0)
		ans = estimate_feat_σ(ex_μs, ex_σs, ex_feat_data)
		return (; true_ans, ans)
	end

	test_estimation()
end

# ╔═╡ 97d046b8-0b96-4a7d-bb05-c401dbf225b7
function estimate_all_feat_σ(μ_hyp, σ_hyp, p_feats_correct)

	best_params = [estimate_feat_σ(μ_hyp, σ_hyp, feat_data) for feat_data in p_feats_correct]
	return (; zip(keys(p_feats_correct), best_params)...)
end

# ╔═╡ 4a05e74c-937b-471d-8a3d-9b3e59977e5f
estimate_all_feat_σ([2], collect(0.5:0.0005:5.0), p_feats_correct[1])

# ╔═╡ 806bede6-ad21-4ea4-aeb0-0bbea00f527d
1.7e-24

# ╔═╡ 2159804d-2b29-4807-9c69-ec77957ab2dd
estimate_all_feat_σ([2], collect(0.5:0.0005:5.0), p_feats_correct[2])

# ╔═╡ e21e7643-06cf-43e8-ade6-192d14aece79
function calc_avg_σ(μ_hyp, σ_hyp, p_feats_correct)
	all_p_params = [estimate_all_feat_σ(μ_hyp, σ_hyp, p_feat) for p_feat in p_feats_correct]

	all_sides_σs = [param.heavy_sides.best_σ for param in all_p_params]
	all_dots_σs = [param.heavy_dots.best_σ for param in all_p_params]
	mean_sides_σ = mean(all_sides_σs)
	mean_dots_σ = mean(all_dots_σs)
	return (; mean_sides_σ, mean_dots_σ)
end

# ╔═╡ 65001704-943e-4ce9-bb6b-c5ae8e3f45ee
estimate_all_feat_σ([2], collect(0.5:0.0005:5.0), p_feats_correct[2])

# ╔═╡ 837acbbb-db56-44c9-bdb4-4ee134fd53fe
# P2, B6 best_σ = 1.2495
# P2, B3 best_σ = 1.115
# Average = 1.18225

# P2, Combined B3 & B6 best_σ = 1.1835

# ╔═╡ 2e7ecc9c-d099-4c2d-a39c-ab3e60670f47
mean([1.2495, 1.115])

# ╔═╡ 32f407a8-1473-493b-9eb1-2ddeabba7c9c
avg_feat_σs = calc_avg_σ([2], collect(0.5:0.005:5.0), p_feats_correct)

# ╔═╡ 9affd7e3-e2b5-4e2a-8815-e10961cfa165
md"""
## Curve Fitting
"""

# ╔═╡ 17d16146-3092-4234-9400-fc62818f2b0c
function generate_model_cdf(μ, σ)
	sides_μ, dots_μ = μ
	sides_σ, dots_σ = σ
	xs = collect(0.0:0.1:4.0)

	sides = cdf(Normal(sides_μ, sides_σ), xs)

	dots = cdf(Normal(dots_μ, dots_σ), xs)
	return (; sides, dots, xs)
end

# ╔═╡ 1439a359-acc2-4fc7-8313-f4abeb92a028
function plot_model(p_data, model_data, avg_σs)
	p_plot = plot_participant_data(p_data[1])
	p_plot = plot_participant_data(p_data[2])

	plot!(p_plot[3], model_data.xs, model_data.sides, label=avg_σs.mean_sides_σ, legend=false)

	plot!(p_plot[4], model_data.xs, model_data.dots, label=avg_σs.mean_dots_σ, legend=false)

	return p_plot
end

# ╔═╡ 1e6843c9-9a98-4f66-87dd-5a653262ed7d
model_data = generate_model_cdf((2, 2), avg_feat_σs)

# ╔═╡ b8bedc21-ffbf-40cc-97ba-c5bb41dd764f
plot_model(p_feats_correct, model_data, avg_feat_σs)

# ╔═╡ 7b4362d9-38ab-48af-9949-6a2fe7ffcb2b
begin
	ex_xs = collect(-3.0:0.1:3.0)
	plot(ex_xs, pdf.(Normal(), ex_xs), lw=3, legend=false)
end

# ╔═╡ b4d25215-e94c-4c12-9614-f07da8901962
σ

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

[compat]
CSV = "~0.10.3"
DataFrames = "~1.3.2"
Distributions = "~0.25.53"
Glob = "~1.3.0"
Plots = "~1.27.3"
"""

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

julia_version = "1.7.2"
manifest_format = "2.0"

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

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

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

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

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

[[deps.CSV]]
deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings"]
git-tree-sha1 = "9310d9495c1eb2e4fa1955dd478660e2ecab1fbb"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.10.3"

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

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

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

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

[[deps.CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

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

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

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

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

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

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

[[deps.Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

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

[[deps.DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "ae02104e835f219b8930c7664b8012c93475c340"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.3.2"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.FilePathsBase]]
deps = ["Compat", "Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "129b104185df66e408edd6625d480b7f9e9823a0"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.18"

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.Glob]]
git-tree-sha1 = "4df9f7e06108728ebf00a0a11edee4b29a482bb2"
uuid = "c27321d9-0574-5035-807b-f59d2c89b15c"
version = "1.3.0"

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

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

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

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

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

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

[[deps.InlineStrings]]
deps = ["Parsers"]
git-tree-sha1 = "61feba885fac3a407465726d0c330b3055df897f"
uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
version = "1.1.2"

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

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

[[deps.InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.3"

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

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

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

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

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

[[deps.PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "28ef6c7ce353f0b35d0df0d5930e0d072c1f5b9b"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.4.1"

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

[[deps.PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "dfb54c4e414caa595a1f2ed759b160f5a3ddcba5"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.3.1"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

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

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

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

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

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

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

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

[[deps.WeakRefStrings]]
deps = ["DataAPI", "InlineStrings", "Parsers"]
git-tree-sha1 = "b1be2855ed9ed8eac54e5caff2afcdb442d52c23"
uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
version = "1.4.2"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# ╔═╡ Cell order:
# ╠═f39f9c31-e702-4de9-81c8-07e43584ec96
# ╠═579cbf1f-be0a-4c80-b4a8-57e7a92a0fc7
# ╠═473dedb8-da32-4687-9b4c-ffb4976f2d54
# ╟─330163b0-ad49-11ec-3128-afba43b7bdd0
# ╟─940b5333-3a5b-4fef-9c0e-356fd4f08354
# ╟─5fa46984-a715-4c29-9373-5af0b62f0c70
# ╟─d01b1239-57d4-4e64-b1c1-7b112387d6a1
# ╠═d3b0223d-5946-4480-a836-ffbc54c035d6
# ╠═14bd5283-578a-4a08-872e-6e73abf9e2c5
# ╟─179c2ed9-664b-4f60-a5d8-e65cd7890062
# ╠═bc5cbbe7-2f9d-45f8-bf65-830468516292
# ╠═44cde6d2-584a-4d33-99ca-b056d8359225
# ╠═9646aea5-73c6-4eda-a9b0-e329f5a1cc9e
# ╠═75382891-f42b-4a5b-a49f-ca145150bd85
# ╟─8a4433de-25cc-4eb9-9252-96f3f1e1917d
# ╠═3558d50a-4972-4483-9e02-65f891bf1d2c
# ╠═0dd8340c-00fb-4e06-9b8d-ec30f8ee753e
# ╠═01880206-93c4-49d6-a7d6-0501c568d7b8
# ╠═257fb4e8-2b0a-4e73-9dad-cb778d144d40
# ╟─3eaf9582-9392-4970-bed2-455cb5756fb0
# ╠═03d82f2a-359f-4adb-8390-27bd144c71b7
# ╠═d436b0ea-fb7c-4e6f-98ff-3c0840099e4d
# ╠═0c839d64-d766-45eb-95a5-db8e27eeb50b
# ╟─dc3c304c-be64-4b45-bab8-f7079b699683
# ╠═e67ca332-e489-4f74-b6e3-1db35a672e23
# ╠═7424b7b8-b6d7-4107-9e0a-04889a5a512d
# ╟─5479abe9-8937-4de9-a31b-afa45d1a2a7e
# ╟─05e90ddc-2ed3-4432-b3ba-da72a4326222
# ╟─8f769d7c-5776-4de2-a614-35c976f3d0ff
# ╟─73f8e4b2-155b-4524-8163-780c48e33135
# ╟─3732a3e2-7760-4abf-b217-61ec019fa27f
# ╠═eca3fc1c-1841-4252-9b40-27e7b20c991a
# ╠═90929ce1-042f-4045-854a-ded3d05cbe2e
# ╠═aec0855f-5630-4620-972e-d6bddf98f95c
# ╠═48ab0022-3dc5-4172-a0f1-16c7b36f0ea9
# ╠═004a7f1d-e4a3-4727-ad44-adbd7de70f10
# ╠═f0463c17-751d-4b4f-b510-4cc4b1934524
# ╠═97d046b8-0b96-4a7d-bb05-c401dbf225b7
# ╠═4a05e74c-937b-471d-8a3d-9b3e59977e5f
# ╠═806bede6-ad21-4ea4-aeb0-0bbea00f527d
# ╠═2159804d-2b29-4807-9c69-ec77957ab2dd
# ╠═e21e7643-06cf-43e8-ade6-192d14aece79
# ╠═65001704-943e-4ce9-bb6b-c5ae8e3f45ee
# ╠═837acbbb-db56-44c9-bdb4-4ee134fd53fe
# ╠═2e7ecc9c-d099-4c2d-a39c-ab3e60670f47
# ╠═32f407a8-1473-493b-9eb1-2ddeabba7c9c
# ╠═9affd7e3-e2b5-4e2a-8815-e10961cfa165
# ╠═17d16146-3092-4234-9400-fc62818f2b0c
# ╠═1439a359-acc2-4fc7-8313-f4abeb92a028
# ╠═1e6843c9-9a98-4f66-87dd-5a653262ed7d
# ╠═b8bedc21-ffbf-40cc-97ba-c5bb41dd764f
# ╠═7b4362d9-38ab-48af-9949-6a2fe7ffcb2b
# ╠═b4d25215-e94c-4c12-9614-f07da8901962
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# U.S. government securities/Treasury constant maturities/Nominal
function readUSTreasuries(dates::StepRange, )

    ## compose url
    url1 =
        "http://www.federalreserve.gov/datadownload/Output.aspx?rel=H15&series=bf17364827e38702b42a58cf8eaa3f78&lastObs=&from="
    url2 = "&to="
    url3 = "&filetype=csv&label=include&layout=seriescolumn"
    dat1 = EconDatasets.convertFedDate(dates[1])
    dat2 = EconDatasets.convertFedDate(dates[end])
    url = string(url1, dat1, url2, dat2, url3)

    ## download and format data
    filename = download(url)
    data = readtable(filename, skipstart = 5,
                     nastrings = ["", "ND"])
    dats = Date[Date(data[ii, 1]) for ii=1:size(data, 1)]
    intRates = Timenum(data[:, 2:end], dats)

    ## get metadata
    metaData = readtable(filename, nrows = 6, header = false,
                         nastrings = ["", "ND"])
    metaDf = EconDatasets.formatMetaData(metaData)

    return (intRates, metaDf)
end

function convertFedDate(dat::Date)
    yyyy = string(Dates.year(dat))
    m = string(Dates.month(dat))
    if length(m) == 1
        m = string("0", m)
    end
    d = string(Dates.day(dat))
    if length(d) == 1
        d = string("0", d)
    end
    return string(m, "/", d, "/", yyyy)
end

function formatMetaData(metaData::DataFrame)
    nVars = size(metaData, 2)-1
    nMeta = size(metaData, 1)
    varNames = Symbol[symbol(metaData[ii, 1]) for ii=1:nMeta]
    df = DataFrame()
    for ii=nMeta:-1:1
        vals = DataArray(typeof(metaData[ii, 2]), nVars)
        for jj=1:nVars
            vals[jj] = metaData[ii, jj+1]
        end
        df[varNames[ii]] = vals
    end
    return df
end

ext: jl


text:
abstract type GAMModel end

function fit!(gm :: Model, X :: Matrix{T},
    y :: Vector{T}) where {T <: Real, Model <: GAMModel} end
function predict(gm :: Model,
    X :: Matrix{T}) where {T <: Real, Model <: GAMModel} end


text:
export AdmittanceModel
"""
An abstract representation of a linear mapping from inputs `x` to outputs `y` of the form
`YΦ = Px`, `y = QᵀΦ`. Subtypes U <: AdmittanceModel are expected to implement:

    get_Y(am::U)
    get_P(am::U)
    get_Q(am::U)
    get_ports(am::U)
    partial_copy(am::U; Y, P, ports)
    compatible(AbstractVector{U})
"""
abstract type AdmittanceModel{T} end

export get_Y
"""
    get_Y(pso::PSOModel)
    get_Y(bbox::Blackbox)

Return a vector of admittance matrices.
"""
function get_Y end

export get_P
"""
    get_P(pso::PSOModel)
    get_P(bbox::Blackbox)

Return an input port matrix.
"""
function get_P end

export get_Q
"""
    get_Q(pso::PSOModel)
    get_Q(bbox::Blackbox)

Return an output port matrix.
"""
function get_Q end

export get_ports
"""
    get_ports(pso::PSOModel)
    get_ports(bbox::Blackbox)

Return a vector of port identifiers.
"""
function get_ports end

export partial_copy
"""
    partial_copy(pso::PSOModel{T, U};
        Y::Union{Vector{V}, Nothing}=nothing,
        P::Union{V, Nothing}=nothing,
        Q::Union{V, Nothing}=nothing,
        ports::Union{AbstractVector{W}, Nothing}=nothing) where {T, U, V, W}

    partial_copy(bbox::Blackbox{T, U};
        Y::Union{Vector{V}, Nothing}=nothing,
        P::Union{V, Nothing}=nothing,
        Q::Union{V, Nothing}=nothing,
        ports::Union{AbstractVector{W}, Nothing}=nothing) where {T, U, V, W}

Create a new model with the same fields except those given as keyword arguments.
"""
function partial_copy end

export compatible
"""
    compatible(psos::AbstractVector{PSOModel{T, U}}) where {T, U}
    compatible(bboxes::AbstractVector{Blackbox{T, U}}) where {T, U}

Check if the models can be cascaded. Always true for PSOModels and true for Blackboxes
that share the same value of `ω`.
"""
function compatible end

export canonical_gauge
"""
    canonical_gauge(pso::PSOModel)
    canonical_gauge(bbox::Blackbox)

Apply an invertible transformation that takes the model to coordinates in which
`P` is `[I ; 0]` (up to floating point errors). Note this will create a dense model.
"""
function canonical_gauge end

import Base: ==
function ==(am1::AdmittanceModel, am2::AdmittanceModel)
    t = typeof(am1)
    if typeof(am2) != t
        return false
    end
    return all([getfield(am1, name) == getfield(am2, name) for name in fieldnames(t)])
end

import Base: isapprox
function isapprox(am1::AdmittanceModel, am2::AdmittanceModel)
    t = typeof(am1)
    if typeof(am2) != t
        return false
    end
    return all([getfield(am1, name) ≈ getfield(am2, name) for name in fieldnames(t)])
end

export apply_transform
"""
    apply_transform(am::AdmittanceModel, transform::AbstractMatrix{<:Number})

Apply a linear transformation `transform` to the coordinates of the model.
"""
function apply_transform(am::AdmittanceModel, transform::AbstractMatrix{<:Number})
    Y = [transpose(transform) * m * transform for m in get_Y(am)]
    P = eltype(Y)(transpose(transform) * get_P(am))
    Q = eltype(Y)(transpose(transform) * get_Q(am))
    return partial_copy(am, Y=Y, P=P, Q=Q)
end

export cascade
"""
    cascade(ams::AbstractVector{U}) where {T, U <: AdmittanceModel{T}}
    cascade(ams::Vararg{U}) where {T, U <: AdmittanceModel{T}}

Cascade the models into one larger block diagonal model.
"""
function cascade(ams::AbstractVector{U}) where {T, U <: AdmittanceModel{T}}
    @assert length(ams) >= 1
    if length(ams) == 1
        return ams[1]
    end
    @assert compatible(ams)
    Y = [cat(m..., dims=(1,2)) for m in zip([get_Y(am) for am in ams]...)]
    P = cat([get_P(am) for am in ams]..., dims=(1,2))
    Q = cat([get_Q(am) for am in ams]..., dims=(1,2))
    ports = vcat([get_ports(am) for am in ams]...)
    return partial_copy(ams[1], Y=Y, P=P, Q=Q, ports=ports)
end

cascade(ams::Vararg{U}) where {T, U <: AdmittanceModel{T}} = cascade(collect(ams))

export ports_to_indices
"""
    ports_to_indices(am::AdmittanceModel{T}, ports::AbstractVector{T}) where T
    ports_to_indices(am::AdmittanceModel{T}, ports::Vararg{T}) where T

Find the indices corresponding to given ports.
"""
function ports_to_indices(am::AdmittanceModel{T}, ports::AbstractVector{T}) where T
    am_ports = get_ports(am)
    return [findfirst(isequal(p), am_ports) for p in ports]
end

ports_to_indices(am::AdmittanceModel{T}, ports::Vararg{T}) where T =
    ports_to_indices(am, collect(ports))

export unite_ports
"""
    unite_ports(am::AdmittanceModel{T}, ports::AbstractVector{T}) where T
    unite_ports(am::AdmittanceModel{T}, ports::Vararg{T}) where T

Unite the given ports into one port.
"""
function unite_ports(am::AdmittanceModel{T}, ports::AbstractVector{T}) where T
    if length(ports) <= 1
        return am
    end
    port_inds = ports_to_indices(am, ports)
    keep_inds = filter(x -> !(x in port_inds[2:end]),
        1:length(get_ports(am))) # keep the first port
    P = get_P(am)
    first_vector = P[:, port_inds[1]]
    constraint_mat = transpose(hcat([first_vector - P[:, i] for i in port_inds[2:end]]...))
    constrained_am = apply_transform(am, nullbasis(constraint_mat))
    return partial_copy(constrained_am, P=get_P(constrained_am)[:, keep_inds],
        Q=get_Q(constrained_am)[:, keep_inds], ports=get_ports(constrained_am)[keep_inds])
end

unite_ports(am::AdmittanceModel{T}, ports::Vararg{T}) where T =
    unite_ports(am, collect(ports))

export open_ports
"""
    open_ports(am::AdmittanceModel{T}, ports::AbstractVector{T}) where T
    open_ports(am::AdmittanceModel{T}, ports::Vararg{T}) where T

Remove the given ports.
"""
function open_ports(am::AdmittanceModel{T}, ports::AbstractVector{T}) where T
    if length(ports) == 0
        return am
    end
    port_inds = ports_to_indices(am, ports)
    keep_inds = filter(x -> !(x in port_inds), 1:length(get_ports(am)))
    return partial_copy(am, P=get_P(am)[:, keep_inds],
        Q=get_Q(am)[:, keep_inds], ports=get_ports(am)[keep_inds])
end

open_ports(am::AdmittanceModel{T}, ports::Vararg{T}) where T =
    open_ports(am, collect(ports))

export open_ports_except
"""
    open_ports_except(am::AdmittanceModel{T}, ports::AbstractVector{T}) where T
    open_ports_except(am::AdmittanceModel{T}, ports::Vararg{T}) where T

Remove all ports except those specified.
"""
function open_ports_except(am::AdmittanceModel{T}, ports::AbstractVector{T}) where T
    return open_ports(am, filter(x -> !(x in ports), get_ports(am)))
end

open_ports_except(am::AdmittanceModel{T}, ports::Vararg{T}) where T =
    open_ports_except(am, collect(ports))

export short_ports
"""
    short_ports(am::AdmittanceModel{T}, ports::AbstractVector{T}) where T
    short_ports(am::AdmittanceModel{T}, ports::Vararg{T}) where T

Replace the given ports by short circuits.
"""
function short_ports(am::AdmittanceModel{T}, ports::AbstractVector{T}) where T
    if length(ports) == 0
        return am
    end
    port_inds = ports_to_indices(am, ports)
    keep_inds = filter(x -> !(x in port_inds), 1:length(get_ports(am)))
    constraint_mat = transpose(hcat([get_P(am)[:, i] for i in port_inds]...))
    constrained_am = apply_transform(am, nullbasis(constraint_mat))
    return partial_copy(constrained_am, P=get_P(constrained_am)[:, keep_inds],
        Q=get_Q(constrained_am)[:, keep_inds], ports=get_ports(constrained_am)[keep_inds])
end

short_ports(am::AdmittanceModel{T}, ports::Vararg{T}) where T =
    short_ports(am, collect(ports))

export short_ports_except
"""
    short_ports_except(am::AdmittanceModel{T}, ports::AbstractVector{T}) where T
    short_ports_except(am::AdmittanceModel{T}, ports::Vararg{T}) where T

Replace all ports with short circuits, except those specified.
"""
function short_ports_except(am::AdmittanceModel{T}, ports::AbstractVector{T}) where T
    return short_ports(am, filter(x -> !(x in ports), get_ports(am)))
end

short_ports_except(am::AdmittanceModel{T}, ports::Vararg{T}) where T =
    short_ports_except(am, collect(ports))

export cascade_and_unite
"""
    cascade_and_unite(models::AbstractVector{U}) where {T, U <: AdmittanceModel{T}}
    cascade_and_unite(models::Vararg{U}) where {T, U <: AdmittanceModel{T}}

Cascade all models and unite ports with the same name.
"""
function cascade_and_unite(models::AbstractVector{U}) where {T, U <: AdmittanceModel{T}}
    @assert length(models) >= 1
    if length(models) == 1
        return models[1]
    end
    # number the ports so that the names are all distinct and then cascade
    port_number = 1
    function rename(model::U)
        ports = [(port_number + i - 1, port) for (i, port) in enumerate(get_ports(model))]
        port_number += length(ports)
        return partial_copy(model, ports=ports)
    end
    model = cascade(map(rename, models))
    # merge all ports with the same name
    original_ports = vcat([get_ports(m) for m in models]...)
    for port in unique(original_ports)
        inds = findall([p[2] == port for p in get_ports(model)])
        model = unite_ports(model, get_ports(model)[inds])
    end
    # remove numbering
    return partial_copy(model, ports=[p[2] for p in get_ports(model)])
end

cascade_and_unite(models::Vararg{U}) where {T, U <: AdmittanceModel{T}} =
    cascade_and_unite(collect(models))


text:
__precompile__(true)

module GlyphSheets

include("dsl.jl")

using .DSL

export create_library, create_typeface, create_glyph, outline_decompose
export close_library, close_typeface
export pt, px

end # module


text:
# # sparseProbVec.jl
#
# export SparseSBPrior, SparseSBPriorP, SparseSBPriorFull,
#   rpost_sparseStickBreak, logSBMmarginal;
#
# struct SparseSBPrior
#   η::Float64
#   p1::Float64
#   μ::Float64
#   M::Float64
# end
#
# mutable struct SparseSBPriorP
#   η::Float64
#   μ::Float64
#   M::Float64
#   a_p1::Float64
#   b_p1::Float64
#   p1_now::Float64
# end
#
# mutable struct SparseSBPriorFull
#   η::Float64
#   M::Float64
#   a_p1::Float64
#   b_p1::Float64
#   a_μ::Float64
#   b_μ::Float64
#   p1_now::Float64
#   μ_now::Float64
# end
#
#
# ### Sparse Stick Breaking mixture
# function h_slice_mu(z2::Vector{Float64}, n2::Int, M::Float64, μ::Float64,
#   logscale::Bool=false)
#
#   a = -n2*(loggamma(M*μ) + loggamma(M*(1.0 - μ)))
#   b = M*μ*sum( log(z2) .- log(1.0 .- z2))
#   out = a + b
#   if (!logscale)
#     out = exp(out)
#   end
#   out
# end
#
#
# function slice_mu(z::Vector{Float64}, ξ::Vector{Int}, μ_old::Float64, M::Float64,
#   a_μ::Float64, b_μ::Float64)
#
#   # assumes xi=1 for group 1, xi=2 for group 2
#
#   n2 = sum( ξ.==2 )
#   betadist = Beta(a_μ, b_μ)
#   if n2 == 0
#     μ_out = rand( betadist )
#   else
#     z2 = copy( z[ findall(ξ.==2) ] )
#     b_u = h_slice_mu(z2, n2, M, μ_old)
#     u = rand( Uniform(0.0, b_u) )
#     lu = log(u)
#
#     ii = 0
#     keepgoing = true
#     μ_cand = 0.0
#     while keepgoing
#       if ii > 1000
#         throw("too many slices!")
#       end
#
#       μ_cand = rand( betadist )
#       lh = h_slice_mu(z2, n2, M, μ_cand, true)
#
#       keepgoing = lu > lh
#       ii += 1
#     end
#
#     μ_out = copy(μ_cand)
#   end
#
#   μ_out
# end
#
#
#
# function rpost_sparseStickBreak(x::Vector{Int}, p1::Float64, η::Float64,
#   μ::Float64, M::Float64, logout::Bool=false)
#
#   K = length(x)
#   n = K - 1
#   rcrx = cumsum( x[range(K, step=-1, length=n)] )[range(n, step=-1, length=n)]  # ∑_{k+1}^K x_k
#   γ = M * μ
#   δ = M * (1.0 - μ)
#
#   ## mixture component 1 update
#   a1 = 1.0 .+ x[1:n]
#   b1 = η .+ rcrx
#
#   ## mixture component 2 update
#   a2 = γ .+ x[1:n]
#   b2 = δ .+ rcrx
#
#   ## calculate posterior mixture weights
#   lgwt1 = log(p1) .- logbeta.(1.0, η) .+ logbeta.(a1, b1)
#   lgwt2 = log(1.0 - p1) .- logbeta.(γ, δ) .+ logbeta.(a2, b2)
#   ldenom = [ logsumexp( [lgwt1[i], lgwt2[i]] ) for i in 1:n ]
#
#   post_p1 = exp.( lgwt1 .- ldenom )
#   is1 = rand(n) .< post_p1
#   ξ = -1*(is1 .* 1 .- 1) .+ 1
#
#   lz = Vector{Float64}(undef, n)
#   for i in 1:n
#     if is1[i]
#       lx1 = log( rand( Distributions.Gamma( a1[i] ) ) )
#       lx2 = log( rand( Distributions.Gamma( b1[i] ) ) )
#       lxm = max(lx1, lx2)
#       lxsum = lxm + log( exp(lx1 - lxm) + exp(lx2 - lxm) ) # logsumexp
#       lz[i] = lx1 - lxsum
#       # z[i] = rand( Beta(a1[i], b1[i]) )
#     else
#       lx1 = log( rand( Distributions.Gamma( a2[i] ) ) )
#       lx2 = log( rand( Distributions.Gamma( b2[i] ) ) )
#       lxm = max(lx1, lx2)
#       lxsum = lxm + log( exp(lx1 - lxm) + exp(lx2 - lxm) ) # logsumexp
#       lz[i] = lx1 - lxsum
#       # z[i] = rand( Beta(a2[i], b2[i]) )
#     end
#   end
#
#   ## break the Stick
#   lw = Vector{Float64}(undef, K)
#   lwhatsleft = 0.0
#
#   for i in 1:n
#     lw[i] = lz[i] + lwhatsleft
#     lwhatsleft += log( 1.0 - exp(lw[i] - lwhatsleft) ) # logsumexp
#   end
#   lw[K] = copy(lwhatsleft)
#
#   z = exp.(lz)
#   w = exp.(lw)
#
#   if logout
#       (lw, lz, ξ)
#     else
#       (w, z, ξ)
#   end
# end
# function rpost_sparseStickBreak(x::Vector{Int}, p1_old::Float64, η::Float64,
#   μ::Float64, M::Float64, a_p1::Float64, b_p1::Float64, logout::Bool=false)
#
#   ## inference for p1 only
#
#   w, z, ξ = rpost_sparseStickBreak(x, p1_old, η, μ, M, logout)
#   p1_now = rand( Beta(a_p1 + sum(ξ==1), b_p1 + sum(ξ==2) ) )
#
#   (w, z, ξ, p1_now) # w and z may be log-valued depending on logout
# end
# function rpost_sparseStickBreak(x::Vector{Int}, p1_old::Float64, η::Float64,
#   μ_old::Float64, M::Float64, a_p1::Float64, b_p1::Float64, a_μ::Float64, b_μ::Float64, logout::Bool=false)
#
#   ## inference for p1 and μ as well
#
#   w, z, ξ = rpost_sparseStickBreak(x, p1_old, η, μ_old, M, logout)
#   if logout
#       μ_now = slice_mu(exp.(z), ξ, μ_old, M, a_μ, b_μ)
#   else
#       μ_now = slice_mu(z, ξ, μ_old, M, a_μ, b_μ)
#   end
#
#   p1_now = rand( Beta(a_p1 + sum(ξ==1), b_p1 + sum(ξ==2) ) )
#
#   (w, z, ξ, μ_now, p1_now) # w and z may be log-valued depending on logout
# end
#
#
#
# """
#     logSBMmarginal(x, p1, η, μ, M)
#
#   Calculate the log of the SBM prior predictive probability mass function.
#
# """
# function logSBMmarginal(x::Vector{Int}, p1::Float64, η::Float64,
#   μ::Float64, M::Float64)
#
#       K = length(x)
#       n = K - 1
#       rcrx = cumsum( x[range(K, step=-1, length=n)] )[range(n, step=-1, length=n)]  # ∑_{k+1}^K x_k
#       γ = M * μ
#       δ = M * (1.0 - μ)
#
#       ## mixture component 1 update
#       a1 = 1.0 .+ x[1:n]
#       b1 = η .+ rcrx
#
#       ## mixture component 2 update
#       a2 = γ .+ x[1:n]
#       b2 = δ .+ rcrx
#
#       ## calculate posterior mixture weights
#       lgwt1 = log(p1) .- logbeta.(1.0, η) .+ logbeta.(a1, b1)
#       lgwt2 = log(1.0 - p1) .- logbeta.(γ, δ) .+ logbeta.(a2, b2)
#       lsum = [ logsumexp( [lgwt1[i], lgwt2[i]] ) for i in 1:n ]
#
#       sum( lsum )
# end


text:
"""
    Nucleotide_Essentials.writeFasta

readFasta(Path::String)
FastaRecord => write_fasta(input_fasta, out, compressed) => .fasta file/.fasta.gz file

Creates a single or multiple entry FastaRecord and outputs either a .fasta or compressed .fasta.gz file to the desired directory

Supported keyword arguments include:

* ```input_fasta::FastaRecord```: A FastaRecord with either a single entry or multiple entries
* ```out::String```: The full or relative path to the directory where files should be written to
* ```compressed::Bool```: Whether or not to write the .fasta files as compressed files or not.
    * If ```true```, files will written as .fasta.gz files
    * If ```false```, files will written as .fasta files

# Example:

```julia
# .fasta files can be written as from an already imported FastaRecord in Julia
myfasta = readFasta("myfasta.fasta");
writeFasta(input_fasta, "example/output/directory, false)

# .fasta files can be written as a .fasta.gz from an already imported FastaRecord in Julia
myfasta = readFasta("myfasta.fasta");
writeFasta(input_fasta, "example/output/directory", true)

# .fasta files with multiple sequences can be read and written as individual .fasta or .fasta.gz in the same step
myfasta = readFasta("myfasta.fasta");
writeFasta(readFasta("/myfasta.fasta"), "example/output/directory", true);
```
"""
function writeFasta(input_fasta::FastaRecord, out::String, compressed::Bool)
    if compressed == true
        for i in 1:length(input_fasta.ID)
            temp_record = string(">", input_fasta.ID[i], "\n", input_fasta.sequence[i])
            out_file = string(out, "/", input_fasta.ID[i], ".fasta")
            write(string(out, "/", input_fasta.ID[i], ".fasta"), temp_record)
            run(`gzip $out_file`); #this works
        end
    else
        for i in 1:length(input_fasta.ID)
            temp_record = string(">", input_fasta.ID[i], "\n", input_fasta.sequence[i])
            write(string(out, "/", input_fasta.ID[i], ".fasta"), temp_record)
        end
    end
end


text:
# BearingRangeTrackingServer

function parseProcTracking!{T <: AbstractString}(instSys::InSituSystem, d::Dict{T, Any})
    # parameters
    lsrnoise = []

    # unpack the data
    b1Dxb = Array{Float64,1}(d["dx"])
    lsrFeats = d["BRfeats"]
    len = length(lsrFeats)
    bfts = zeros(3,len)
    for i in 1:len, j in 1:3
      bfts[j,i] = Float64(lsrFeats[i][j])
    end
    # do the work
    propAllTrackers!(instSys.trackers, b1Dxb, [0.05;0.05;0.004])
    hardassc = assocMeasWFeats!(instSys.trackers, bfts)
    measUpdateTrackers!(instSys.trackers, hardassc, [0.4;0.02])

    # pack outgoing data
    od = Dict{Int64,Array{Float64,1}}()
    for t in instSys.trackers
      m = getKDEMean(t[2].bel) #getKDEMax
      od[t[1]] = cart2pol(m, [1.0;1.0])[1]
      # println("tracker $(t[1]), pol=$(od[t[1]])")
    end
    println("tracking $(length(instSys.trackers)) feats")
    rmsg = json(od)
    return rmsg # json(current tracked feature jsoned dict)
end

function parseTCP!(insitu::InSituSystem, line::AbstractString)
    # sp = split(line,' ')
    # cmd = sp[1]

    f = +
    goahead = true
    d = Dict{AbstractString, Any}()
    try
      d = JSON.parse(line)
      cmd = AbstractString(d["CMD"])
      # if cmd == "INIT"
      #   f = parseInit!
      if cmd == "PROCESS"
        f = parseProcTracking!
      # elseif cmd == "RESET"
      #   f = parseReset
      elseif cmd == "QUIT"
        println("parseTCP -- should quit now")
        return false, string("")
      else
        warn("parseTCP -- I don't know what $(cmd) means")
        goahead = false
      end
    catch
      warn("parseTCP! -- json parse and command failed: $(line)")
      goahead = false
    end
    retstr = nothing
    goahead ? retstr = f(insitu, d) : nothing
    return true, retstr
end

# function propAllTrackers!(trkrs::Dict{Int64,Feature}, bDxb1::Array{Float64,1}, s::Array{Float64,1})

function tcpStringBRTrackingServer(;port::Int=60002)
  println("Empty teacker object created")
  instSys = makeInSituSys(zeros(3), zeros(2,0))

  println("Listenting on $(port)")
  server = listen(port)
  loop = true
  while loop
   sock = accept(server)
     while isopen(sock)
         gotline = readline(sock)[1:(end-1)]
         loop, retstr = parseTCP!(instSys, gotline)
         loop ? ( isopen(sock) ? println(sock, retstr) : nothing ) : close(sock)
     end
     println("connection lost, loop=$(loop)")
  end
  !loop ? close(server) : nothing
  return instSys
end


text:
using SampleChainsDynamicHMC
using Test

@testset "SampleChainsDynamicHMC.jl" begin
    include("notebook.jl")
end


text:
"""
    Bounds.MultiEllipsoid([T=Float64], ndims)
    Bounds.MultiEllipsoid(::AbstractVector{Ellipsoid})

Use multiple [`Ellipsoid`](@ref)s in an optimal clustering to bound prior space. This implementation follows the MultiNest implementation outlined in Feroz et al. (2008,2009).[^3][^4] For more details about the bounding algorithm, see the extended help (`??Bounds.MultiEllipsoid`)

[^3]: Feroz and Hobson, 2008, MNRAS 384, 2 ["Multimodal nested sampling: an efficient and robust alternative to Markov Chain Monte Carlo methods for astronomical data analyses"](https://academic.oup.com/mnras/article/384/2/449/1023701)
[^4]: Feroz et al., 2009, MNRAS 398, 4 ["MultiNest: an efficient and robust Bayesian inference tool for cosmology and particle physics"](https://academic.oup.com/mnras/article/398/4/1601/981502)

## Extended help

The multiple-ellipsoidal implementation is defined as follows:

1. Fit a [`Bounds.Ellipsoid`](@ref) to the sample.
2. Perform K-means clustering (here using [Clustering.jl](https://github.com/JuliaStats/Clustering.jl)) centered at the endpoints of the bounding ellipsoid. This defines two clusters within the sample.
3. If either cluster has fewer than two points, consider it ill-defined and end any recursion.
4. Fit [`Bounds.Ellipsoid`](@ref) to each of the clusters assigned in (2).
5. If the volume of the parent ellipsoid is more than twice the volume of the two child ellipsoids, recurse (1-5) to each child.

To sample from this distribution, a random ellipsoid is selected, and a random sample is sampled from that ellipsoid. We then reverse this and find all of the ellipsoids which enclose the sampled point, and select one of those randomly for the enclosing bound.
"""
struct MultiEllipsoid{T} <: AbstractBoundingSpace{T}
    ellipsoids::Vector{Ellipsoid{T}}
end
MultiEllipsoid(ndims::Integer) = MultiEllipsoid(Float64, ndims)
MultiEllipsoid(T::Type, ndims::Integer) = MultiEllipsoid([Ellipsoid(T, ndims)])

Base.eltype(::MultiEllipsoid{T}) where {T} = T
Base.broadcastable(me::MultiEllipsoid) = (me,)
Base.length(me::MultiEllipsoid) = length(me.ellipsoids)
Base.ndims(me::MultiEllipsoid) = ndims(me.ellipsoids[1])

volume(me::MultiEllipsoid) = sum(volume, me.ellipsoids)

function scale!(me::MultiEllipsoid, factor)
    scale!.(me.ellipsoids, factor)
    return me
end

function fit(::Type{<:MultiEllipsoid}, x::AbstractMatrix; pointvol = 0)
    parent = fit(Ellipsoid, x, pointvol = pointvol)
    ells = fit(MultiEllipsoid, x, parent, pointvol = pointvol)
    return MultiEllipsoid(ells)
end

function fit(::Type{<:MultiEllipsoid}, x::AbstractMatrix, parent::Ellipsoid; pointvol = 0)
    ndim, npoints = size(x)

    # Clustering will fail with fewer than k=2 points
    npoints ≤ 2 && return [parent]

    p1, p2 = endpoints(parent)
    starting_points = hcat(p1, p2)
    R = kmeans!(x, starting_points; maxiter = 10)
    labels = assignments(R)
    x1 = x[:, labels .== 1]
    x2 = x[:, labels .== 2]

    # if either cluster has fewer than ndim points, it is ill-defined
    if size(x1, 2) < 2ndim || size(x2, 2) < 2ndim
        return [parent]
    end

    # Getting bounding ellipsoid for each cluster
    ell1, ell2 = fit.(Ellipsoid, (x1, x2), pointvol = pointvol)

    # If total volume decreased by over half, recurse
    if volume(ell1) + volume(ell2) < 0.5volume(parent)
        return vcat(fit(MultiEllipsoid, x1, ell1, pointvol = pointvol),
                    fit(MultiEllipsoid, x2, ell2, pointvol = pointvol))
    end

    # Otherwise see if total volume is much larger than expected
    # and split into more than 2 clusters
    if volume(parent) > 2npoints * pointvol
        out = vcat(fit(MultiEllipsoid, x1, ell1, pointvol = pointvol),
                    fit(MultiEllipsoid, x2, ell2, pointvol = pointvol))
        sum(volume, out) < 0.5volume(parent) && return out
    end

    # Otherwise, return single bounding ellipse
    return [parent]
end

Base.in(x::AbstractVector, me::MultiEllipsoid) = any(ell->x ∈ ell, me.ellipsoids)

function Base.rand(rng::AbstractRNG, me::MultiEllipsoid)
    length(me) == 1 && return rand(rng, me.ellipsoids[1])

    vols = volume.(me.ellipsoids)
    weights = vols ./ sum(vols)
    local x

    while true
        # Select random ellipsoid
        idx = rand(rng, Categorical(weights))
        ell = me.ellipsoids[idx]

        # Select point
        x = rand(rng, ell)

        # How many ellipsoids is the sample in
        n = count(ell -> x ∈ ell, me.ellipsoids)

        # Only accept with probability 1/n
        (n == 1 || rand(rng) < 1 / n) && break
    end

    return x
end

"""
Returns a random live point and a bounding ellipsoid, since MultiEllipsoid doesn't have
valid transformation axes
"""
function rand_live(rng::AbstractRNG, me::MultiEllipsoid, us)
    idx = rand(rng, Base.axes(us, 2))
    u = us[:, idx]

    # find which Ellipsoid/s it overlaps with
    idxs = findall(ell -> u ∈ ell, me.ellipsoids)
    # TODO if point isn't bounded, update bounds
    if isempty(idxs)
        return u, nothing
    end

    # pick random encompassing ellipsoid
    idx = rand(rng, idxs)

    return u, me.ellipsoids[idx]
end


text:
################################################################################
#
#  NfMaxOrd/LinearAlgebra.jl : Linear algebra over maximal orders
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
# (C) 2016 Tommy Hofmann
#
################################################################################

using Hecke

export pseudo_matrix, pseudo_hnf

function _det_bound(M::GenMat{NfOrdElem{NfMaxOrd}})
  n = rows(M)
  O = base_ring(M)
  d = degree(O)
  c1, c2 = Hecke.base_change_const(O)

  return fmpz(BigInt(ceil(sqrt(c2)*c1^(n/2)*BigInt(n)^n*BigInt(d)^n*BigInt(_max_max(M))^n)))
end

function _max_max(M::GenMat{NfOrdElem{NfMaxOrd}})
  d = FlintZZ(1)
  for i in 1:rows(M)
    for j in 1:cols(M)
      if !iszero(M[i, j])
        v = elem_in_basis(M[i, j])
        for k in degree(base_ring(M))
          d = max(d, abs(v[k]))
        end
      end
    end
  end
  return d
end

#function det(M::GenMat{NfOrdElem{NfMaxOrd}})
#  O = base_ring(M)::NfMaxOrd
#  B = _det_bound(M)
#  p = next_prime(2^60) # magic numbers
#  P = fmpz(1)
#  i = 1
#  res = zero(O)
#  t = 0.0
#  while P < 2*B
#    # reject some bad primes
#    if true
#      #println("$p");
#      Omodp, pi_p = quo(O, ideal(O, p))
#      Mmodp = MatrixSpace(Omodp, rows(M), cols(M))(M)
#      t += @elapsed detmodp = pi_p\Hecke.det(Mmodp)
#      if i == 1
#        res = detmodp
#      elseif i >= 2
#        g, e, f = gcdx(P, fmpz(p))
#        @assert isone(g)
#        res = f*p*res + e*P*detmodp
#        res = mod_sym(res, P*p)
#        #@assert _basis_coord_nonneg(res)
#      end
#      P = P*p
#      p = next_prime(p)
#      i += 1
#    end
#  end
#  #println("Modular determinant time: $t");
#  return res
#end

function _get_coeff_raw(x::nmod_poly, i::Int)
  u = ccall((:nmod_poly_get_coeff_ui, :libflint), UInt, (Ptr{nmod_poly}, Int), &x, i)
  return u
end

function det(M::GenMat{NfOrdElem{NfMaxOrd}}, b = 62)
  O = base_ring(M)::NfMaxOrd
  B = _det_bound(M)
  p = next_prime(2^b) # magic numbers
  P = fmpz(1)
  i = 1
  res = O()
  t_det = 0.0
  t_reducing = 0.0
  t_splitting = 0.0
  used_primes = Array{UInt, 1}()
  tmp_polys = Array{nmod_poly, 1}()

  Zx = PolynomialRing(FlintIntegerRing(), "x")[1]
  tmp_fmpz_poly = Zx()

  while P < 2*B
    # reject some bad primes
    if is_index_divisor(O, p)
      continue
    end

    ttt = @elapsed lp = prime_decomposition(O, p)

    ramified = false

    for j in 1:length(ttt)
      if lp[j][2] > 1 # || lp[j][1].splitting_type[2] > 3
        ramified = true
        break
      end
    end

    if ramified
      continue
    end

    push!(used_primes, UInt(p))

    t_splitting += ttt

    s = nmod_poly(UInt(p))
    t = nmod_poly(UInt(p))
    g = nmod_poly(UInt(p))
    z = nmod_poly(UInt(p))
    detmodp = nmod_poly(UInt(p))
    polyprod = nmod_poly(UInt(p))
    tmp_nmod_poly = nmod_poly(UInt(p))

    j = 1
    for (Q, e) in lp
      F, mF = ResidueField(O, Q)

      MmodQ = MatrixSpace(F, rows(M), cols(M))()
      t_reducing += @elapsed __helper!(MmodQ, mF, M.entries)
      t_det += @elapsed detmodQ = det(MmodQ)

      # now the CRT step

      if j == 1
        polyprod = mF.poly_of_the_field
        ccall((:nmod_poly_set, :libflint), Void, (Ptr{nmod_poly}, Ptr{fq_nmod}), &detmodp, &detmodQ)
      else
        crt!(tmp_nmod_poly, detmodp, detmodQ, polyprod, mF.poly_of_the_field, s, t)
        mul!(polyprod, polyprod, mF.poly_of_the_field)
        set!(detmodp, tmp_nmod_poly)
      end

      j += 1
    end

    push!(tmp_polys, detmodp)

    P = P*p
    p = next_prime(p)
    i += 1
  end

  # now the CRT on each coefficient

  fc = fmpz_comb(used_primes)
  fct = fmpz_comb_temp(fc)
  zz = fmpz()

  @assert length(used_primes) == length(tmp_polys)

  for i in 0:degree(O)
    fmpz_multi_crt_ui!(zz, [ _get_coeff_raw(tmp_polys[j], i) for j in 1:length(tmp_polys) ], fc, fct)
    setcoeff!(tmp_fmpz_poly, i, zz)
  end

  tut = fmpq_poly(tmp_fmpz_poly)
  tut.parent = parent(nf(O).pol)
  res = mod_sym(O(nf(O)(tut)), P)

  #println("Modular determinant time: $t_det");
  #println("Time for reducing: $t_reducing");
  #println("Time for splitting: $t_splitting");
  #println("Used $i primes")

  return res
end

# s, t are auxillary variables, r1, r2 are the residues, m1, m2 are the moduli
# aliasing is not allowed (?)
function crt!(z::nmod_poly, r1::nmod_poly, r2::Union{nmod_poly, fq_nmod}, m1::nmod_poly, m2::nmod_poly, s::nmod_poly, t::nmod_poly)
  ccall((:nmod_poly_xgcd, :libflint), Void, (Ptr{nmod_poly}, Ptr{nmod_poly}, Ptr{nmod_poly}, Ptr{nmod_poly}, Ptr{nmod_poly}), &z, &s, &t, &m1, &m2)
  @assert Bool(ccall((:nmod_poly_is_one, :libflint), Cint, (Ptr{nmod_poly}, ), &z))
  # z = s*m1*r2 + t*m2*r1
  mul!(z, s, m1)
  ccall((:nmod_poly_mul, :libflint), Void, (Ptr{nmod_poly}, Ptr{nmod_poly}, Ptr{fq_nmod}), &z, &z, &r2)
  mul!(t, t, m2)
  mul!(t, t, r1)
  add!(z, z, t)
  mul!(s, m1, m2)
  rem!(z, z, s)
end

function set!(z::nmod_poly, x::nmod_poly)
  ccall((:nmod_poly_set, :libflint), Void, (Ptr{nmod_poly}, Ptr{nmod_poly}), &z, &x)
end

function __helper!(z, mF, entries)
  s = size(entries)
  for i in 1:s[2]
    for j in 1:s[1]
      z[j, i] = mF(entries[j, i])
    end
  end
  return z
end

function mod_sym(x, m)
  z = elem_in_basis(x)
  for i in 1:length(z)
    z[i] = mod(z[i], m)
    if 2*z[i] > m
      z[i] = z[i] - m
    end
  end
  return parent(x)(z)
end

function _basis_coord_nonneg(x)
  b = elem_in_basis(x)
  for i in 1:length(b)
    if b[i] < 0
      return false
    end
  end
  return true
end

function rand(M::GenMatSpace{NfOrdElem{NfMaxOrd}}, B::Union{Int, fmpz})
  z = M()
  for i in 1:rows(z)
    for j in 1:cols(z)
      z[i, j] = Hecke.rand(Hecke.base_ring(M), B)
    end
  end
  return z
end

type PMat
  parent
  matrix::GenMat{nf_elem}
  coeffs::Array{NfMaxOrdFracIdl, 1}

  function PMat(m::GenMat{nf_elem}, c::Array{NfMaxOrdFracIdl, 1})
    z = new()
    z.matrix = m
    z.coeffs = c
    return z
  end
end

function show(io::IO, P::PMat)
  print(io, "Pseudo-matrix over $(parent(P.matrix[1, 1]))\n")
  for i in 1:rows(P.matrix)
    print(io, "$(P.coeffs[i]) with row vector $(sub(P.matrix, i:i, 1:cols(P.matrix)))\n")
  end
end

function PseudoMatrix(m::GenMat{nf_elem}, c::Array{NfMaxOrdFracIdl, 1})
  # sanity checks
  return PMat(m ,c)
end

function PseudoMatrix(m::GenMat{NfOrdElem{NfMaxOrd}}, c::Array{NfMaxOrdIdl, 1})
  mm = change_ring(m, nf(base_ring(m)))
  cc = map(z -> NfMaxOrdFracIdl(z, fmpz(1)), c)
  return PMat(mm, cc)
end

function rows(m::PMat)
  return rows(m.matrix)
end

function cols(m::PMat)
  return cols(m.matrix)
end

function change_ring(m::GenMat{NfOrdElem{NfMaxOrd}}, K::AnticNumberField)
  return MatrixSpace(K, rows(m), cols(m))(map(z -> K(z), m.entries))
end

function det(m::PMat)
  z = m.coeffs[1]
  for i in 2:rows(m)
    if isdefined(z.num, :gen_two) && isdefined(m.coeffs[i].num, :gen_two)
    end
    z = z*m.coeffs[i]
  end
  return det(m.matrix)*z
end

# this is slow
function _coprime_integral_ideal_class(x::NfMaxOrdFracIdl, y::NfMaxOrdIdl)
  O = order(y)
  c = conjugates_init(nf(O).pol)
  #num_x_inv = inv(num(x))
  x_inv = inv(x)
  check = true
  z = ideal(O, O(1))
  a = nf(O)()
  i = 0
  while check
    i += 1
    a = rand(x_inv, 10)
    b = x*a
    z = divexact(num(b), den(b))
    norm(z + y) == 1 ? (check = false) : (check = true)
  end
  return z, a
end

# this is slow
function _coprime_norm_integral_ideal_class(x::NfMaxOrdFracIdl, y::NfMaxOrdIdl)
  O = order(y)
  c = conjugates_init(nf(O).pol)
  #num_x_inv = inv(num(x))
  x_inv = inv(x)
  check = true
  z = ideal(O, O(1))
  a = nf(O)()
  i = 0
  while check
    i += 1
    a = rand(x_inv, 10)
    b = x*a
    z = divexact(num(b), den(b))
    gcd(norm(z), norm(y)) == 1 ? (check = false) : (check = true)
  end
  @assert gcd(norm(z), norm(y)) == 1
  return z, a
end

function rand(I::NfMaxOrdIdl, B::Int)
  r = rand(1:B, degree(order(I)))
  b = basis(I)
  z = r[1]*b[1]
  for i in 2:degree(order(I))
    z = z + r[i]*b[i]
  end
  return z
end

function rand(I::NfMaxOrdFracIdl, B::Int)
  z = rand(num(I), B)
  return divexact(elem_in_nf(z), den(I))
end

function pseudo_hnf(P::PMat, m::NfMaxOrdIdl, shape::Symbol = :upperright)
  O = order(m)

  t_comp_red = 0.0
  t_mod_comp = 0.0
  t_sum = 0.0
  t_div = 0.0
  t_idem = 0.0

  t_comp_red += @elapsed z = _matrix_for_reduced_span(P, m)
  t_mod_comp += @elapsed zz = strong_echelon_form(z, shape)

  res_mat = MatrixSpace(nf(O), rows(P), cols(P))()
  for i in 1:rows(P)
    for j in 1:cols(P)
      res_mat[i, j] = elem_in_nf(zz[i, j].elem)
    end
  end

  res = PMat(res_mat, [ deepcopy(x)::NfMaxOrdFracIdl for x in P.coeffs])

  for i in 1:rows(P)
    if iszero(zz[i, i].elem)
      res.matrix[i, i] = one(nf(O))
      res.coeffs[i] = NfMaxOrdFracIdl(deepcopy(m), fmpz(1))
    else
      o = ideal(O, zz[i, i].elem)
      t_sum += @elapsed g = o + m
      t_div += @elapsed oo = divexact(o, g)
      t_div += @elapsed mm = divexact(m, g)
      t_idem += @elapsed e, f = idempotents(oo, mm)
      res.coeffs[i] = NfMaxOrdFracIdl(deepcopy(g), fmpz(1))
      mul_row!(res.matrix, i, elem_in_nf(e))
      divide_row!(res.matrix, i, elem_in_nf(zz[i, i].elem))
      res.matrix[i, i] = one(nf(O))
    end
  end

  #println("computation of reduction : $t_comp_red")
  #println("modular computation      : $t_mod_comp")
  #println("computation of ideal sum : $t_sum")
  #println("computation of ideal div : $t_div")
  #println("computation of idems     : $t_idem")

  return res
end

#this is Algorithm 4 of FH2014
# we assume that span(P) \subseteq O^r
function _matrix_for_reduced_span(P::PMat, m::NfMaxOrdIdl)
  O = order(m)
  c = Array(NfMaxOrdIdl, rows(P))
  mat = deepcopy(P.matrix)
  for i in 1:rows(P)
    I, a = _coprime_norm_integral_ideal_class(P.coeffs[i], m)
    divide_row!(mat, i, a)
    c[i] = I
  end
  Om, OtoOm = quo(O, m)
  z = MatrixSpace(Om, rows(P), cols(P))()
  for i in 1:rows(z)
    for j in 1:cols(z)
      @assert norm(c[i])*mat[i, j] in O
      @assert euclid(OtoOm(O(norm(c[i])))) == 1
      q = OtoOm(O(norm(c[i])*mat[i,j]))
      qq = inv(OtoOm(O(norm(c[i]))))
      z[i, j] = q*qq
    end
  end
  return z
end

_check(a) = a.has_coord ? dot(a.elem_in_basis, basis(parent(a))) == a : true

function _check(m::GenMat{NfOrdElem{NfMaxOrd}})
  for i in 1:rows(m)
    for j in 1:cols(m)
      if !_check(m[i, j].elem)
        println("$i $j not consistent: $(m[i, j]) : $(m[i, j].elem_in_basis)")
        error("dasdsad")
      end
    end
  end
end

function _check(m::GenMat{NfMaxOrdQuoRingElem})
  for i in 1:rows(m)
    for j in 1:cols(m)
      if !_check(m[i, j].elem)
        println("$i $j not consistent: $(m[i, j].elem) : $(m[i, j].elem.elem_in_basis)")
        error("dasdsad")
      end
    end
  end
end

function divide_row!{T}(M::GenMat{T}, i::Int, r::T)
  for j in 1:cols(M)
    M[i, j] = divexact(M[i, j], r)
  end
end

function mul_row!{T}(M::GenMat{T}, i::Int, r::T)
  for j in 1:cols(M)
    M[i, j] = M[i, j] * r
  end
end

function _contained_in_span_of_pseudohnf(v::GenMat{nf_elem}, P::PMat)
  # assumes that P is upper right pseudo-HNF
  w = deepcopy(v)
  for i in 1:rows(P)
    if !(w[1, i] in P.coeffs[i])
      return false
    end
    e = w[1, i]
    for j in 1:cols(P)
      w[1, j] = w[1, j] - e*P.matrix[i, j]
    end
  end
  if !iszero(w)
    return false
  end
  return true
end

function _contained_in_span_of_pseudohnf_lowerleft(v::GenMat{nf_elem}, P::PMat)
  # assumes that P is lowerleft pseudo-HNF
  w = deepcopy(v)
  for i in rows(P):-1:1
    if !(w[1, i] in P.coeffs[i])
      return false
    end
    e = w[1, i]
    for j in 1:cols(P)
      w[1, j] = w[1, j] - e*P.matrix[i, j]
    end
  end
  if !iszero(w)
    return false
  end
  return true
end

function _spans_subset_of_pseudohnf(M::PMat, P::PMat)
# P upperright
  for i in rows(M)
    A = M.coeffs[i]
    v = sub(M.matrix, i:i, 1:cols(P))
    for b in basis(A.num)
      bb = divexact(elem_in_nf(b), A.den)
      if !Hecke._contained_in_span_of_pseudohnf(bb*v, P)
        return false
      end
    end
  end
  return true
end

function sub(M::GenMat, rows::UnitRange{Int}, cols::UnitRange{Int})
  @assert step(rows) == 1 && step(cols) == 1
  z = MatrixSpace(base_ring(M), length(rows), length(cols))()
  for i in rows
    for j in cols
      z[i - start(rows) + 1, j - start(cols) + 1] = M[i, j]
    end
  end
  return z
end

function in(x::nf_elem, y::NfMaxOrdFracIdl)
  B = inv(basis_mat(y))
  O = order(y)
  M = MatrixSpace(FlintZZ, 1, degree(O))()
  t = FakeFmpqMat(M)
  elem_to_mat_row!(t.num, 1, t.den, x)
  v = t*basis_mat_inv(O)
  v = v*B

  return v.den == 1
end

pseudo_matrix(x...) = PseudoMatrix(x...)


text:
# Eckstein, Bertsekas, "On the Douglas-Rachford Splitting Method and the
# Proximal Point Algorithm for Maximal Monotone Operators",
# Mathematical Programming, vol. 55, no. 1, pp. 293-318 (1989).

# this is the plug and play version of DR
# Definition of the DR recursion from Eckstein, Bertsekas is:
# z^{k+1} = J_{\gamma A}( (2J_{\gamma B} -I)z^k ) + (I - J_{\gamma B})z^k

# This unravels to
# x^{k+1} = J_{\gamma B} z^k
# y^{k+1} = J_{\gamma A}( 2x^{k+1} -z^k )
# z^{k+1} = z^k + y^{k+1} - x^{k+1}

using Base.Iterators
using ProximalAlgorithms.IterationTools
using ProximalOperators: Zero
using LinearAlgebra
using Printf

Base.@kwdef struct PnpDrsIteration{R,C<:Union{R,Complex{R}},Tx<:AbstractArray{C},F<:Function,G<:Function}
    J_A!::F
    J_B!::G
    u0::Tx #initial point
    #gamma::R #don't think gamma is needed
end

Base.IteratorSize(::Type{<:PnpDrsIteration}) = Base.IsInfinite()

Base.@kwdef struct PnpDrsState{Tx}
    u::Tx
    x::Tx = similar(u) #empty
    y::Tx = similar(u)
    r::Tx = similar(u)
    res::Tx = similar(u)
end

function Base.iterate(iter::PnpDrsIteration, state::PnpDrsState = PnpDrsState(u=copy(iter.u0)))
    iter.J_B!(state.x,state.u)
    @. state.r = 2 *state.x - state.u
    iter.J_A!(state.y,state.r)
    state.res .= state.y .- state.x
    state.u .+= state.res
    return state, state
end


text:
include("parseCpp.jl")

function  main()
    code = """
   /* test etdthad;lkfjtashd
atdsavapsodfiuase */

if ( y < 60)
    // do stuff
    printf(test);

else if (y > 49)
    // do more stuff
    test = 32 + 3;
else
    // do most stuff
    yadda ydadda yadda;

if ( y < 32){
    // do stuff

    test text for the block
}
else if (y > 42) {
    // do more stuff

    more teset text for this stuff
}
else {
    // do most stuff
    wohooo
    test text for this
}
"""

    len = length(code)
    # test find pattern functions
    println("newline is at index ",findPattern("\n",code,1, len))

    # filter out comments
    code_wo_comments = filterComments(code)

    i = 1
    str = ""
    while str != "done"
        # test block finding function
        i,str = findBlock(code_wo_comments,i+1)
        println(i," ",str)
    end

    output = parseCpp(code)

    for block in output

        println(block[2].true_code)

    end

end

main()

max_stars_repo_name: aklapatch/CodeVisualized

max_issues_repo_name: aklapatch/CodeVisualized

max_forks_repo_name: aklapatch/CodeVisualized