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

using Markdown
using InteractiveUtils

# ╔═╡ 7ac0c689-419a-428b-9006-69d3e6c1ea7e
using PlutoUI; TableOfContents(title = "目录", aside = false, depth = 4)

# ╔═╡ d0dfdb8d-2b7a-4956-802d-73922db32330
begin
	using DataFrames, CSV
	using FreqTables
	using Statistics: mean, std
	using GLM
	using CategoricalArrays
end

# ╔═╡ b36a8b22-992c-4b50-a8ef-23852f73e26c
let
	using CairoMakie #GLMakie
	set_theme!(theme_minimal())
	update_theme!(
		backgroundcolor = RGBf(0.98, 0.98, 0.98),
		resolution = (1024, 768)
	 )
end

# ╔═╡ e8a9ae4d-438a-4496-b75c-ee9a5e91d740
md"""
### 载入包
"""

# ╔═╡ b7a11e97-fba8-4656-b529-ff6f21a75a0f
md"""
### 实验一
"""

# ╔═╡ ab814eb8-b384-40d0-95b2-139031baaa35
md"""
#### 读入数据
"""

# ╔═╡ 3e149cf5-4387-4a05-a13c-4b00167965b5
md"""
- trialType: complete表示完全反馈条件； partial表示部分反馈条件
- x1，y1：左边转盘的两个分数
- p，p2：x1的概率为p1；y1的概率为p2
- x2，y2：右边转盘的两个分数
- q，q2：x2的概率为q1；y2的概率为q2
- left points：选择左边转盘得到的分数
- right points：选择右边转盘得到的分数
- LeftYourUn：选择左边转盘没有得到的分数
- RightYouUn：选择右边转盘没有得到的分数
- `key_resp_4.keys`：被试选择了左边转盘（left）还是右边转盘（right）。因此从key_resp_4.keys和left  points、right points就能看出被试这轮得到的分数。
- rating：被试的情绪评分
- total score：被试在转盘上累计的总分，也可以通过下一轮的总分减去上一轮的总分得出被试上一轮得到的分数。
"""

# ╔═╡ 2a9d311c-a766-4f32-ac8b-64a9b63cedcc
csv_list1 = filter(endswith(".csv"), readdir("data/Experiment_1", join = true));

# ╔═╡ 7235b483-8484-4ca2-a700-2dbea497ac26
df1 = mapreduce(vcat, csv_list1[1:41]) do csv
	CSV.read(csv, DataFrame, stringtype = String,
			 drop = (i, name) -> startswith(string(name), "Column")
	)
end;

# ╔═╡ 91006353-4ea7-4fd0-98c2-f23e6a38a351
dropmissing!(df1, [Symbol("key_resp_4.keys"), :leftPoints, :rightPoints]);

# ╔═╡ 06e153e3-f076-4e16-8ab6-141cc2b0d7bb
function neg2pos2(key, left, right)
	if key == "left"
		if left + right < 0
			return (OwnScore = left + 200, OtherScore = right + 200, Gain = "Lose")
		elseif left + right >= 0
			return (OwnScore = left,       OtherScore = right,       Gain = "Gain")
		end
	elseif key == "right"
		if left + right < 0
			return (OwnScore = right + 200, OtherScore = left + 200, Gain = "Lose")
		elseif left + right >= 0
			return (OwnScore = right,       OtherScore = left,       Gain = "Gain")
		end
	end
end

# ╔═╡ 30ad6818-714b-4362-b8b4-af193dcd70cf
transform!(df1, [Symbol("key_resp_4.keys"), :leftPoints, :rightPoints] =>
		ByRow(neg2pos2) => AsTable
);

# ╔═╡ 401c49d8-4096-480f-8271-5e37e9a307f5
select!(df1,
	["participant", "trialType", "OwnScore", "OtherScore", "Gain", "rating"]
);

# ╔═╡ 80ad622e-0783-4b7c-942d-73adb7eff2e8
rename!(df1, "participant" => "Participant",
	"trialType" => "Feedback", "rating" => "Emotion");

# ╔═╡ 1ba1c644-a1d9-4233-9085-3ac871d95a9e
md"""
#### 总体结果
"""

# ╔═╡ 9c8b0137-d0af-4a6e-b2ef-28d30bd80c28
df1c = combine(
	groupby(df1, [:Participant, :Feedback, :Gain]),
	:Emotion => mean => :Emotion
);

# ╔═╡ 436f5420-3ef0-4a9c-8fc7-b20ea714a889
e1fm1 = lm(@formula(Emotion ~ Feedback * Gain), df1c);

# ╔═╡ e5d58b40-7c37-43c8-a83a-5127b9b0ea38
e1fm2 = lm(@formula(Emotion ~ Feedback + Gain), df1c);

# ╔═╡ a1e0d71c-01a4-42bd-af08-993369cd74b9
e1fm3 = lm(@formula(Emotion ~ Feedback), df1c);

# ╔═╡ 83dfce03-65c2-4d0e-989d-5da22471a024
e1fm4 = lm(@formula(Emotion ~ Gain), df1c);

# ╔═╡ b14a6b03-768d-4bac-8dab-96db011a2acd
ftest(e1fm1.model, e1fm2.model)

# ╔═╡ a61067b2-83b9-4dbb-aa5f-94aaa399ac43
coeftable(e1fm2)

# ╔═╡ 12d7c32f-4a62-4cae-a9b2-39ef597bd084
let
	replace!(df1c.Gain,     "Gain"     => "获益框架", "Lose"    => "损失框架");
	replace!(df1c.Feedback, "complete" => "完全反馈", "partial" => "不完全反馈");
	transform!(df1c,
		:Feedback => (x -> categorical(x, ordered = true)) => :Feedback
	)
	transform!(df1c,
		:Gain => (x -> categorical(x, ordered = true)) => :Gain
	)
	Gain = unique(df1c.Gain)
	Feedback = unique(df1c.Feedback)

	fig = Figure()
	for (id1, gg) in enumerate(Gain)
		ax = Axis(fig[1, id1],
			title = gg, xlabel = "反馈类型", ylabel = "被试情绪自评")
		ylims!(-20, 40)
		dt = df1c[df1c.Gain .== gg, :]
		for (id2, fd) in enumerate(Feedback)
			dtt = dt[dt.Feedback .== fd, :]
			xx = Int.(dtt.Feedback.refs)
			yy = dtt.Emotion
			boxplot!(ax, xx, yy, width = 1, label = fd, whiskerwidth = 0.5)
		end
		ax.xticks = (1:length(levels(Feedback)), levels(Feedback))
		id1 > 1 && hideydecorations!(ax)
	end
	fig
end

# ╔═╡ 76d4e45d-225e-4c4c-9f5e-3f30db7c53d9
md"""
#### 被试得分
"""

# ╔═╡ 1abbbc4d-5213-4632-9f3f-f45ad638b015
df1d = combine(
	groupby(df1, [:Participant, :Feedback, :Gain, :OwnScore]),
	:Emotion => mean => :Emotion
);

# ╔═╡ 5bd3d6c1-2187-4d24-b354-244ffcb9d41b
e1fm51 = lm(@formula(Emotion ~ Feedback * Gain * OwnScore), df1d);

# ╔═╡ 6b4a960b-205c-4e92-bc7e-ac4d36728316
e1fm52 = lm(@formula(Emotion ~ Feedback + Gain + OwnScore), df1d);

# ╔═╡ bda05cb8-6e31-4c5d-a0df-7b0eadefb1ce
ftest(e1fm51.model, e1fm52.model)

# ╔═╡ d0178ff1-4562-4be2-84d1-5f73ff561c91
coeftable(e1fm52)

# ╔═╡ 4b7d5041-7657-410b-abd6-2cbcb3ae7eb8
let
	replace!(df1d.Gain,     "Gain"     => "获益框架", "Lose"    => "损失框架");
	replace!(df1d.Feedback, "complete" => "完全反馈", "partial" => "不完全反馈");
	fig = Figure()
	for (id1, ff) in enumerate(unique(df1d.Feedback))
		dtt = df1d[df1d.Feedback .== ff, :]
		ax = Axis(fig[1, id1], title = ff,
			xlabel = "被试得分", ylabel = "被试情绪自评",
			xticks = -200:50:200, yticks = -40:20:40)
		for (id2, gg) in enumerate(unique(df1d.Gain))
			dt = dtt[dtt.Gain .== gg, :]
			xx = dt.OwnScore .+ (iseven(id2) ? 1 : -1) * 6.5
			yy = dt.Emotion
			boxplot!(ax, xx, yy, width = 15, label = gg, whiskerwidth = 0.5)
		end
		id1 > 1 && hideydecorations!(ax)
	end
	axislegend("框架", position = :rb)
	fig
end

# ╔═╡ 396cbce1-4ab7-44f6-ae14-5ae041bf487a
md"""
### 实验二
"""

# ╔═╡ 29c07409-5d62-42aa-a002-b14dfd7c1685
md"""
#### 读入数据
"""

# ╔═╡ 3223d964-3424-4a91-87cd-690876844b4d
md"""
相关说明
- `mychoice`: 被试本试次得分
- `otherchoice`： 电脑本试次得分
- `stay`: 剩余箱子得分
- `key_resp_4.keys`： 电脑先选， 被试按键
- `key_resp_6.keys`： 被试先选， 被试按键
- `rating`: 情绪评分
- `choicepos`: 电脑先选条件下， 电脑选择的位置
"""

# ╔═╡ af32e1d2-b575-11ec-34f6-25b42c355184
cd(@__DIR__) # pwd()

# ╔═╡ 0f5bd14c-638a-4eb1-85e5-729767c5fc4f
md"""
#### 总体结果
"""

# ╔═╡ d5967d18-62d7-4870-b3b7-af7416acb2b4
md"""
#### 被试得分
"""

# ╔═╡ ebc5cfbd-dd1d-4089-8798-e35357b38cac
md"""
### 实验三
"""

# ╔═╡ 4e5c0d80-2871-4f57-a0a0-acbffb4ec199
md"""
#### 读入数据
"""

# ╔═╡ f5d03c3a-3854-4c08-9cf7-074dd628393b
md"""
#### 总体结果
"""

# ╔═╡ dec0718e-799d-48ce-b444-cf07daf803c6
md"""
#### 被试得分
"""

# ╔═╡ af57ef9e-0030-4b2d-8b1b-127f50d9754a
md"""
### 附录
"""

# ╔═╡ ef982930-6388-466d-b0e1-8d715316a884
function list2df(experiment)

	# 1. List csv files
	csv_list = filter(endswith(".csv"), readdir(joinpath("data", experiment), join = true))

	# 2. Read data to dataframe
	df = mapreduce(vcat, csv_list) do csv
		CSV.read(csv, DataFrame, stringtype = String,
				 drop = (i, name) -> startswith(string(name), "Column")
		)
	end

	# 3. Delete rows
	dropmissing!(df, [:mychoice, :otherchoice, :stay, :rating])

	# 4. Create new columns
	## a. Define a function to convert negative values to postive ones
		function neg2pos(x, y, z)
			if (x + y + z) < 0
				mn = -200
				return (Human = x-mn, Computer = y-mn, Remain = z-mn, Gain = "Lose")
			else
				return (Human = x,    Computer = y,    Remain = z,    Gain = "Gain")
			end
		end

	## b. Do the converting
	transform!(df,
		[:mychoice, :otherchoice, :stay] => ByRow(neg2pos)    => AsTable,
		:choicepos => ByRow(x -> ismissing(x) ? "HC" : "CH")  => :Order,
		[:mychoice, :stay] => ByRow( (x, y) -> x - y)         => :HumanRemainDiff,
		[:mychoice, :otherchoice] => ByRow( (x, y) -> x - y)  => :HumanComputerDiff
	)

	# 5. Rename column names
	rename!(df, "rating" => "Emotion", "participant"  => "Participant")

	# 6. Select usefule columns, "mychoice", "otherchoice", "stay"
	select!(df,
		["Participant", "Gain", "Order", "Emotion", "Human", "Computer", "Remain",
		"HumanRemainDiff", "HumanComputerDiff"]
	)

	return df
end

# ╔═╡ e2d72e3d-b16e-4782-92fb-743fd332f202
ex2 = list2df("Experiment_2");

# ╔═╡ 43f57dbc-3d7e-48e0-bce1-23e058938bd1
size(ex2, 1) == 192 * 35 # 192 trial * 35 participant

# ╔═╡ 218d3696-9ee4-4841-ad53-447849a238f9
freqtable(ex2, :Gain, :Order)

# ╔═╡ ad4a9897-84fa-4e42-a653-fb36bb2e7657
ex2t = groupby(ex2, [:Participant, :Gain, :Order]);

# ╔═╡ a90937a1-1961-465a-bd7f-f8f35da9badc
ex2ts = combine(ex2t, :Emotion => mean => :Emotion);

# ╔═╡ 2eaf3273-0ac0-4cc8-8879-fd15903c75de
fmt1 = lm(@formula(Emotion ~ Gain * Order), ex2ts);

# ╔═╡ 699dd367-d4a3-431d-87a1-dd6ff276c84b
fmt2 = lm(@formula(Emotion ~ Gain), ex2ts);

# ╔═╡ a0691d8b-d80c-4b5c-b78c-163ee13faf85
ftest(fmt1.model, fmt2.model)

# ╔═╡ 71e9a4b4-5dbb-4a67-b898-0cfffd90998d
coeftable(fmt2)

# ╔═╡ cc35db1a-715b-419a-b5a7-3f4a3f8faab4
let
	replace!(ex2ts.Gain, "Gain" => "获益框架", "Lose" => "损失框架")
	transform!(ex2ts,
		:Gain => (x -> categorical(x, ordered = true)) => :Gain
	)
	xx = ex2ts.Gain
	yy = ex2ts.Emotion
	colors = Makie.default_palettes.patchcolor.val
	boxplot(Int.(xx.refs), yy, whiskerwidth = 0.5,
		color = colors[Int.(xx.refs)],
		axis = (xticks = (1:length(levels(xx)), levels(xx)),
			xlabel = "框架", ylabel = "被试情绪自评"
		)
	)
end

# ╔═╡ 1ebe900d-2301-44ea-86b3-0fd7f03769b1
HM = combine(
	groupby(ex2, [:Participant, :Gain, :Human]),
	:Emotion => mean => :Emotion
);

# ╔═╡ a6d31357-bb6d-4fb2-a642-112f39d32832
fHM1 = lm(@formula(Emotion ~ Gain * Human), HM);

# ╔═╡ ac88d890-6514-4e82-b776-ff98467d8708
fHM2 = lm(@formula(Emotion ~ Gain + Human), HM);

# ╔═╡ 241d663c-dc41-4039-b9a2-0f48d792b8a1
ftest(fHM1.model, fHM2.model)

# ╔═╡ 2a111cc4-a80c-43d0-83b3-7c827761a4fa
coeftable(fHM2)

# ╔═╡ 5d2801db-a987-48e5-b71d-6674c0f7af58
let
	replace!(HM.Gain, "Gain" => "获益框架", "Lose" => "损失框架")
	transform!(HM,
		:Gain => (x -> categorical(x, ordered = true)) => :Gain
	)
	fig = Figure()
	ax = Axis(fig[1, 1],
		xlabel = "被试得分", ylabel = "被试情绪自评", xticks = -200:50:200)
	for (id, gg) in enumerate(unique(HM.Gain))
		dt = HM[HM.Gain .== gg, :]
		xx = dt.Human .+ (iseven(id) ? 1 : -1) * 6.5
		yy = dt.Emotion
		boxplot!(ax, xx, yy, width = 15, label = gg, whiskerwidth = 0.5)
	end
	axislegend("框架", position = :rb)
	fig
end

# ╔═╡ af3c116c-510d-4c5c-9b64-fb56a106e5eb
ex2n = list2df("Experiment_2_N");

# ╔═╡ fe698e99-939b-44d4-aaed-1bdfdc5c58f0
ex2nts = combine(
	groupby(ex2n, [:Participant, :Gain, :Order]),
	:Emotion => mean => :Emotion);

# ╔═╡ b4f7a35f-0f74-465b-82a0-31c404b49a10
ex2n1 = lm(@formula(Emotion ~ Gain * Order), ex2nts);

# ╔═╡ 7516d09f-7356-41ac-a565-2c4e2c21408d
ex2n2 = lm(@formula(Emotion ~ Gain + Order), ex2nts);

# ╔═╡ 3e754434-5934-42df-861f-614bfa19643f
ftest(ex2n1.model, ex2n2.model)

# ╔═╡ 00762ff1-c860-4d2c-b5ef-4a27b28b5d7f
coeftable(ex2n2)

# ╔═╡ dccfe48e-4068-4f35-8446-8476033f90e6
let
	replace!(ex2nts.Gain, "Gain" => "获益框架", "Lose" => "损失框架")
	transform!(ex2nts,
		:Gain => (x -> categorical(x, ordered = true)) => :Gain
	)
	xx = ex2nts.Gain
	yy = ex2nts.Emotion
	colors = Makie.default_palettes.patchcolor.val
	boxplot(Int.(xx.refs), yy, whiskerwidth = 0.5, color = colors[Int.(xx.refs)],
		axis = (xticks = (1:length(levels(xx)), levels(xx)),
			xlabel = "框架", ylabel = "被试情绪自评")
	)
end

# ╔═╡ c36cb78e-09ce-4521-803f-3e94d5ebd7cd
HMn = combine(
	groupby(ex2n, [:Participant, :Gain, :Human]),
	:Emotion => mean => :Emotion
);

# ╔═╡ fa2bd5ca-1fc5-4706-86b1-c7baf5bd2a47
fHMn1 = lm(@formula(Emotion ~ Gain * Human), HMn);

# ╔═╡ 4a90cefa-9793-481d-ae1b-863a379369ec
fHMn2 = lm(@formula(Emotion ~ Gain + Human), HMn);

# ╔═╡ dc676a84-430c-4316-9e42-e935d3dc9e74
ftest(fHMn1.model, fHMn2.model)

# ╔═╡ 366feb2b-b08e-4176-9623-149b6cbaf989
coeftable(fHMn2)

# ╔═╡ d0099447-0c7d-4573-b9b0-ff18398f36a4
let
	replace!(HMn.Gain, "Gain" => "获益框架", "Lose" => "损失框架")
	transform!(HMn,
		:Gain => (x -> categorical(x, ordered = true)) => :Gain
	)
	fig = Figure()
	ax = Axis(fig[1, 1],
		xlabel = "被试得分", ylabel = "被试情绪自评", xticks = -200:50:200)
	for (id, gg) in enumerate(unique(HMn.Gain))
		dt = HMn[HMn.Gain .== gg, :]
		xx = dt.Human .+ (iseven(id) ? 1 : -1) * 6.5
		yy = dt.Emotion
		boxplot!(ax, xx, yy, width = 15, label = gg, whiskerwidth = 0.5)
	end
	axislegend("框架", position = :rb)
	fig
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
FreqTables = "da1fdf0e-e0ff-5433-a45f-9bb5ff651cb1"
GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[compat]
CSV = "~0.10.4"
CairoMakie = "~0.7.5"
CategoricalArrays = "~0.10.5"
DataFrames = "~1.3.2"
FreqTables = "~0.4.5"
GLM = "~1.6.2"
PlutoUI = "~0.7.38"
"""

# ╔═╡ 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.AbstractFFTs]]
deps = ["ChainRulesCore", "LinearAlgebra"]
git-tree-sha1 = "6f1d9bc1c08f9f4a8fa92e3ea3cb50153a1b40d4"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.1.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

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

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

[[deps.Animations]]
deps = ["Colors"]
git-tree-sha1 = "e81c509d2c8e49592413bfb0bb3b08150056c79d"
uuid = "27a7e980-b3e6-11e9-2bcd-0b925532e340"
version = "0.4.1"

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

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

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

[[deps.Automa]]
deps = ["Printf", "ScanByte", "TranscodingStreams"]
git-tree-sha1 = "d50976f217489ce799e366d9561d56a98a30d7fe"
uuid = "67c07d97-cdcb-5c2c-af73-a7f9c32a568b"
version = "0.8.2"

[[deps.AxisAlgorithms]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"]
git-tree-sha1 = "66771c8d21c8ff5e3a93379480a2307ac36863f7"
uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950"
version = "1.0.1"

[[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.CEnum]]
git-tree-sha1 = "215a9aa4a1f23fbd05b92769fdd62559488d70e9"
uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
version = "0.4.1"

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

[[deps.Cairo]]
deps = ["Cairo_jll", "Colors", "Glib_jll", "Graphics", "Libdl", "Pango_jll"]
git-tree-sha1 = "d0b3f8b4ad16cb0a2988c6788646a5e6a17b6b1b"
uuid = "159f3aea-2a34-519c-b102-8c37f9878175"
version = "1.0.5"

[[deps.CairoMakie]]
deps = ["Base64", "Cairo", "Colors", "FFTW", "FileIO", "FreeType", "GeometryBasics", "LinearAlgebra", "Makie", "SHA", "StaticArrays"]
git-tree-sha1 = "4a0de4f5aa2d5d27a1efa293aeabb1a081e46b2b"
uuid = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
version = "0.7.5"

[[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.CategoricalArrays]]
deps = ["DataAPI", "Future", "Missings", "Printf", "Requires", "Statistics", "Unicode"]
git-tree-sha1 = "109664d3a6f2202b1225478335ea8fea3cd8706b"
uuid = "324d7699-5711-5eae-9e2f-1d82baa6b597"
version = "0.10.5"

[[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.ColorBrewer]]
deps = ["Colors", "JSON", "Test"]
git-tree-sha1 = "61c5334f33d91e570e1d0c3eb5465835242582c4"
uuid = "a2cac450-b92f-5266-8821-25eda20663c8"
version = "0.4.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.ColorVectorSpace]]
deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "TensorCore"]
git-tree-sha1 = "3f1f500312161f1ae067abe07d13b40f78f32e07"
uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4"
version = "0.9.8"

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

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

[[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.EllipsisNotation]]
deps = ["ArrayInterface"]
git-tree-sha1 = "d064b0340db45d48893e7604ec95e7a2dc9da904"
uuid = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
version = "1.5.0"

[[deps.Expat_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "bad72f730e9e91c08d9427d5e8db95478a3c323d"
uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
version = "2.4.8+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.FFTW]]
deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"]
git-tree-sha1 = "505876577b5481e50d089c1c68899dfb6faebc62"
uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
version = "1.4.6"

[[deps.FFTW_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c6033cc3892d0ef5bb9cd29b7f2f0331ea5184ea"
uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
version = "3.3.10+0"

[[deps.FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "80ced645013a5dbdc52cf70329399c35ce007fae"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.13.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.FreeType]]
deps = ["CEnum", "FreeType2_jll"]
git-tree-sha1 = "cabd77ab6a6fdff49bfd24af2ebe76e6e018a2b4"
uuid = "b38be410-82b0-50bf-ab77-7b57e271db43"
version = "4.0.0"

[[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.FreeTypeAbstraction]]
deps = ["ColorVectorSpace", "Colors", "FreeType", "GeometryBasics"]
git-tree-sha1 = "b5c7fe9cea653443736d264b85466bad8c574f4a"
uuid = "663a7486-cb36-511b-a19d-713bb74d65c9"
version = "0.9.9"

[[deps.FreqTables]]
deps = ["CategoricalArrays", "Missings", "NamedArrays", "Tables"]
git-tree-sha1 = "488ad2dab30fd2727ee65451f790c81ed454666d"
uuid = "da1fdf0e-e0ff-5433-a45f-9bb5ff651cb1"
version = "0.4.5"

[[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.GLM]]
deps = ["Distributions", "LinearAlgebra", "Printf", "Reexport", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "StatsModels"]
git-tree-sha1 = "609115155b0dc532fa5130de65ed086efd27bfbd"
uuid = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
version = "1.6.2"

[[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.Graphics]]
deps = ["Colors", "LinearAlgebra", "NaNMath"]
git-tree-sha1 = "1c5a84319923bea76fa145d49e93aa4394c73fc2"
uuid = "a2bd30eb-e257-5431-a919-1863eab51364"
version = "1.1.1"

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

[[deps.GridLayoutBase]]
deps = ["GeometryBasics", "InteractiveUtils", "Observables"]
git-tree-sha1 = "169c3dc5acae08835a573a8a3e25c62f689f8b5c"
uuid = "3955a311-db13-416c-9275-1d80ed98e5e9"
version = "0.6.5"

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

[[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.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.IfElse]]
git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.1"

[[deps.ImageCore]]
deps = ["AbstractFFTs", "ColorVectorSpace", "Colors", "FixedPointNumbers", "Graphics", "MappedArrays", "MosaicViews", "OffsetArrays", "PaddedViews", "Reexport"]
git-tree-sha1 = "9a5c62f231e5bba35695a20988fc7cd6de7eeb5a"
uuid = "a09fc81d-aa75-5fe9-8630-4744c3626534"
version = "0.9.3"

[[deps.ImageIO]]
deps = ["FileIO", "JpegTurbo", "Netpbm", "OpenEXR", "PNGFiles", "QOI", "Sixel", "TiffImages", "UUIDs"]
git-tree-sha1 = "464bdef044df52e6436f8c018bea2d48c40bb27b"
uuid = "82e4d734-157c-48bb-816b-45c225c6df19"
version = "0.6.1"

[[deps.Imath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "87f7662e03a649cffa2e05bf19c303e168732d3e"
uuid = "905a6f67-0a94-5f89-b386-d35d92009cd1"
version = "3.1.2+0"

[[deps.IndirectArrays]]
git-tree-sha1 = "012e604e1c7458645cb8b436f8fba789a51b257f"
uuid = "9b13fd28-a010-5f03-acff-a1bbcff69959"
version = "1.0.0"

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

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

[[deps.IntelOpenMP_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d979e54b71da82f3a65b62553da4fc3d18c9004c"
uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
version = "2018.0.3+2"

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

[[deps.Interpolations]]
deps = ["AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"]
git-tree-sha1 = "b15fc0a95c564ca2e0a7ae12c1f095ca848ceb31"
uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
version = "0.13.5"

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

[[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.Isoband]]
deps = ["isoband_jll"]
git-tree-sha1 = "f9b6d97355599074dc867318950adaa6f9946137"
uuid = "f1662d9f-8043-43de-a69a-05efc1cc6ff4"
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]]
deps = ["CEnum", "FileIO", "ImageCore", "JpegTurbo_jll", "TOML"]
git-tree-sha1 = "a77b273f1ddec645d1b7c4fd5fb98c8f90ad10a5"
uuid = "b835a17e-a41a-41e7-81f0-2f016b05efe0"
version = "0.1.1"

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

[[deps.KernelDensity]]
deps = ["Distributions", "DocStringExtensions", "FFTW", "Interpolations", "StatsBase"]
git-tree-sha1 = "591e8dc09ad18386189610acafb970032c519707"
uuid = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
version = "0.6.3"

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

[[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.LazyArtifacts]]
deps = ["Artifacts", "Pkg"]
uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"

[[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.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.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.MKL_jll]]
deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
git-tree-sha1 = "e595b205efd49508358f7dc670a940c790204629"
uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
version = "2022.0.0+0"

[[deps.Makie]]
deps = ["Animations", "Base64", "ColorBrewer", "ColorSchemes", "ColorTypes", "Colors", "Contour", "Distributions", "DocStringExtensions", "FFMPEG", "FileIO", "FixedPointNumbers", "Formatting", "FreeType", "FreeTypeAbstraction", "GeometryBasics", "GridLayoutBase", "ImageIO", "IntervalSets", "Isoband", "KernelDensity", "LaTeXStrings", "LinearAlgebra", "MakieCore", "Markdown", "Match", "MathTeXEngine", "Observables", "OffsetArrays", "Packing", "PlotUtils", "PolygonOps", "Printf", "Random", "RelocatableFolders", "Serialization", "Showoff", "SignedDistanceFields", "SparseArrays", "StaticArrays", "Statistics", "StatsBase", "StatsFuns", "StructArrays", "UnicodeFun"]
git-tree-sha1 = "63de3b8a5c1f764e4e3a036c7752a632b4f0b8d1"
uuid = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
version = "0.16.6"

[[deps.MakieCore]]
deps = ["Observables"]
git-tree-sha1 = "c5fb1bfac781db766f9e4aef96adc19a729bc9b2"
uuid = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
version = "0.2.1"

[[deps.MappedArrays]]
git-tree-sha1 = "e8b359ef06ec72e8c030463fe02efe5527ee5142"
uuid = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900"
version = "0.4.1"

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

[[deps.Match]]
git-tree-sha1 = "1d9bc5c1a6e7ee24effb93f175c9342f9154d97f"
uuid = "7eb4fadd-790c-5f42-8a69-bfa0b872bfbf"
version = "1.2.0"

[[deps.MathTeXEngine]]
deps = ["AbstractTrees", "Automa", "DataStructures", "FreeTypeAbstraction", "GeometryBasics", "LaTeXStrings", "REPL", "RelocatableFolders", "Test"]
git-tree-sha1 = "70e733037bbf02d691e78f95171a1fa08cdc6332"
uuid = "0a4f8689-d25c-4efe-a92b-7142dfc1aa53"
version = "0.2.1"

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

[[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.MosaicViews]]
deps = ["MappedArrays", "OffsetArrays", "PaddedViews", "StackViews"]
git-tree-sha1 = "b34e3bc3ca7c94914418637cb10cc4d1d80d877d"
uuid = "e94cdb99-869f-56ef-bcf0-1ae2bcbe0389"
version = "0.3.3"

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

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

[[deps.NamedArrays]]
deps = ["Combinatorics", "DataStructures", "DelimitedFiles", "InvertedIndices", "LinearAlgebra", "Random", "Requires", "SparseArrays", "Statistics"]
git-tree-sha1 = "2fd5787125d1a93fbe30961bd841707b8a80d75b"
uuid = "86f7a689-2022-50b4-a561-43c23ac3c673"
version = "0.9.6"

[[deps.Netpbm]]
deps = ["FileIO", "ImageCore"]
git-tree-sha1 = "18efc06f6ec36a8b801b23f076e3c6ac7c3bf153"
uuid = "f09324ee-3d7c-5217-9330-fc30815ba969"
version = "1.0.2"

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

[[deps.Observables]]
git-tree-sha1 = "fe29afdef3d0c4a8286128d4e45cc50621b1e43d"
uuid = "510215fc-4207-5dde-b226-833fc4488ee2"
version = "0.4.0"

[[deps.OffsetArrays]]
deps = ["Adapt"]
git-tree-sha1 = "043017e0bdeff61cfbb7afeb558ab29536bbb5ed"
uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
version = "1.10.8"

[[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.OpenEXR]]
deps = ["Colors", "FileIO", "OpenEXR_jll"]
git-tree-sha1 = "327f53360fdb54df7ecd01e96ef1983536d1e633"
uuid = "52e1d378-f018-4a11-a4be-720524705ac7"
version = "0.3.2"

[[deps.OpenEXR_jll]]
deps = ["Artifacts", "Imath_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "923319661e9a22712f24596ce81c54fc0366f304"
uuid = "18a262bb-aa17-5467-a713-aee519bc75cb"
version = "3.1.1+0"

[[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.PNGFiles]]
deps = ["Base64", "CEnum", "ImageCore", "IndirectArrays", "OffsetArrays", "libpng_jll"]
git-tree-sha1 = "eb4dbb8139f6125471aa3da98fb70f02dc58e49c"
uuid = "f57f5aa1-a3ce-4bc8-8ab9-96f992907883"
version = "0.3.14"

[[deps.Packing]]
deps = ["GeometryBasics"]
git-tree-sha1 = "1155f6f937fa2b94104162f01fa400e192e4272f"
uuid = "19eb6ba3-879d-56ad-ad62-d5c202156566"
version = "0.4.2"

[[deps.PaddedViews]]
deps = ["OffsetArrays"]
git-tree-sha1 = "03a7a85b76381a3d04c7a1656039197e70eda03d"
uuid = "5432bcbf-9aad-5242-b902-cca2824c8663"
version = "0.5.11"

[[deps.Pango_jll]]
deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "FriBidi_jll", "Glib_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "3a121dfbba67c94a5bec9dde613c3d0cbcf3a12b"
uuid = "36c8627f-9965-5494-a995-c6b170f724f3"
version = "1.50.3+0"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "621f4f3b4977325b9128d5fae7a8b4829a0c2222"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.4"

[[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.PkgVersion]]
deps = ["Pkg"]
git-tree-sha1 = "a7a7e1a88853564e551e4eba8650f8c38df79b37"
uuid = "eebad327-c553-4316-9ea0-9fa01ccd7688"
version = "0.1.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.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "670e559e5c8e191ded66fa9ea89c97f10376bb4c"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.38"

[[deps.PolygonOps]]
git-tree-sha1 = "77b3d3605fc1cd0b42d95eba87dfcd2bf67d5ff6"
uuid = "647866c9-e3ac-4575-94e7-e3d426903924"
version = "0.1.2"

[[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.ProgressMeter]]
deps = ["Distributed", "Printf"]
git-tree-sha1 = "d7a7aef8f8f2d537104f170139553b14dfe39fe9"
uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
version = "1.7.2"

[[deps.QOI]]
deps = ["ColorTypes", "FileIO", "FixedPointNumbers"]
git-tree-sha1 = "18e8f4d1426e965c7b532ddd260599e1510d26ce"
uuid = "4b34888f-f399-49d4-9bb3-47ed5cae4e65"
version = "1.0.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.Ratios]]
deps = ["Requires"]
git-tree-sha1 = "dc84268fe0e3335a62e315a3a7cf2afa7178a734"
uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439"
version = "0.4.3"

[[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.SIMD]]
git-tree-sha1 = "7dbc15af7ed5f751a82bf3ed37757adf76c32402"
uuid = "fdea26ae-647d-5447-a871-4b548cad5224"
version = "3.4.1"

[[deps.ScanByte]]
deps = ["Libdl", "SIMD"]
git-tree-sha1 = "9cc2955f2a254b18be655a4ee70bc4031b2b189e"
uuid = "7b38b023-a4d7-4c5e-8d43-3f3097f304eb"
version = "0.3.0"

[[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.ShiftedArrays]]
git-tree-sha1 = "22395afdcf37d6709a5a0766cc4a5ca52cb85ea0"
uuid = "1277b4bf-5013-50f5-be3d-901d8477a67a"
version = "1.0.0"

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

[[deps.SignedDistanceFields]]
deps = ["Random", "Statistics", "Test"]
git-tree-sha1 = "d263a08ec505853a5ff1c1ebde2070419e3f28e9"
uuid = "73760f76-fbc4-59ce-8f25-708e95d2df96"
version = "0.4.0"

[[deps.Sixel]]
deps = ["Dates", "FileIO", "ImageCore", "IndirectArrays", "OffsetArrays", "REPL", "libsixel_jll"]
git-tree-sha1 = "8fb59825be681d451c246a795117f317ecbcaa28"
uuid = "45858cf5-a6b0-47a3-bbea-62219f50df47"
version = "0.1.2"

[[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.StackViews]]
deps = ["OffsetArrays"]
git-tree-sha1 = "46e589465204cd0c08b4bd97385e4fa79a0c770c"
uuid = "cae243ae-269e-4f55-b966-ac2d0dc13c15"
version = "0.1.1"

[[deps.Static]]
deps = ["IfElse"]
git-tree-sha1 = "87e9954dfa33fd145694e42337bdd3d5b07021a6"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.6.0"

[[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 = "72e6abd6fc9ef0fa62a159713c83b7637a14b2b8"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.17"

[[deps.StatsModels]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Printf", "REPL", "ShiftedArrays", "SparseArrays", "StatsBase", "StatsFuns", "Tables"]
git-tree-sha1 = "03c99c7ef267c8526953cafe3c4239656693b8ab"
uuid = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
version = "0.6.29"

[[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.TensorCore]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "1feb45f88d133a655e001435632f019a9a1bcdb6"
uuid = "62fd8b95-f654-4bbd-a8a5-9c27f68ccd50"
version = "0.1.1"

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

[[deps.TiffImages]]
deps = ["ColorTypes", "DataStructures", "DocStringExtensions", "FileIO", "FixedPointNumbers", "IndirectArrays", "Inflate", "OffsetArrays", "PkgVersion", "ProgressMeter", "UUIDs"]
git-tree-sha1 = "aaa19086bc282630d82f818456bc40b4d314307d"
uuid = "731e570b-9d59-4bfa-96dc-6df516fadf69"
version = "0.5.4"

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

[[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.WeakRefStrings]]
deps = ["DataAPI", "InlineStrings", "Parsers"]
git-tree-sha1 = "b1be2855ed9ed8eac54e5caff2afcdb442d52c23"
uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
version = "1.4.2"

[[deps.WoodburyMatrices]]
deps = ["LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "de67fa59e33ad156a590055375a30b23c40299d3"
uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6"
version = "0.5.5"

[[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_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_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_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.isoband_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "51b5eeb3f98367157a7a12a1fb0aa5328946c03c"
uuid = "9a68df92-36a6-505f-a73e-abb412b6bfb4"
version = "0.2.3+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.libsixel_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "78736dab31ae7a53540a6b752efc61f77b304c5b"
uuid = "075b6546-f08a-558a-be8f-8157d0f608a5"
version = "1.8.6+1"

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

# ╔═╡ Cell order:
# ╟─7ac0c689-419a-428b-9006-69d3e6c1ea7e
# ╟─e8a9ae4d-438a-4496-b75c-ee9a5e91d740
# ╠═d0dfdb8d-2b7a-4956-802d-73922db32330
# ╠═b36a8b22-992c-4b50-a8ef-23852f73e26c
# ╟─b7a11e97-fba8-4656-b529-ff6f21a75a0f
# ╟─ab814eb8-b384-40d0-95b2-139031baaa35
# ╟─3e149cf5-4387-4a05-a13c-4b00167965b5
# ╠═2a9d311c-a766-4f32-ac8b-64a9b63cedcc
# ╠═7235b483-8484-4ca2-a700-2dbea497ac26
# ╠═91006353-4ea7-4fd0-98c2-f23e6a38a351
# ╟─06e153e3-f076-4e16-8ab6-141cc2b0d7bb
# ╠═30ad6818-714b-4362-b8b4-af193dcd70cf
# ╠═401c49d8-4096-480f-8271-5e37e9a307f5
# ╠═80ad622e-0783-4b7c-942d-73adb7eff2e8
# ╟─1ba1c644-a1d9-4233-9085-3ac871d95a9e
# ╠═9c8b0137-d0af-4a6e-b2ef-28d30bd80c28
# ╠═436f5420-3ef0-4a9c-8fc7-b20ea714a889
# ╠═e5d58b40-7c37-43c8-a83a-5127b9b0ea38
# ╠═a1e0d71c-01a4-42bd-af08-993369cd74b9
# ╠═83dfce03-65c2-4d0e-989d-5da22471a024
# ╠═b14a6b03-768d-4bac-8dab-96db011a2acd
# ╠═a61067b2-83b9-4dbb-aa5f-94aaa399ac43
# ╟─12d7c32f-4a62-4cae-a9b2-39ef597bd084
# ╟─76d4e45d-225e-4c4c-9f5e-3f30db7c53d9
# ╠═1abbbc4d-5213-4632-9f3f-f45ad638b015
# ╠═5bd3d6c1-2187-4d24-b354-244ffcb9d41b
# ╠═6b4a960b-205c-4e92-bc7e-ac4d36728316
# ╠═bda05cb8-6e31-4c5d-a0df-7b0eadefb1ce
# ╠═d0178ff1-4562-4be2-84d1-5f73ff561c91
# ╟─4b7d5041-7657-410b-abd6-2cbcb3ae7eb8
# ╟─396cbce1-4ab7-44f6-ae14-5ae041bf487a
# ╟─29c07409-5d62-42aa-a002-b14dfd7c1685
# ╟─3223d964-3424-4a91-87cd-690876844b4d
# ╠═af32e1d2-b575-11ec-34f6-25b42c355184
# ╠═e2d72e3d-b16e-4782-92fb-743fd332f202
# ╠═43f57dbc-3d7e-48e0-bce1-23e058938bd1
# ╠═218d3696-9ee4-4841-ad53-447849a238f9
# ╟─0f5bd14c-638a-4eb1-85e5-729767c5fc4f
# ╠═ad4a9897-84fa-4e42-a653-fb36bb2e7657
# ╠═a90937a1-1961-465a-bd7f-f8f35da9badc
# ╠═2eaf3273-0ac0-4cc8-8879-fd15903c75de
# ╠═699dd367-d4a3-431d-87a1-dd6ff276c84b
# ╠═a0691d8b-d80c-4b5c-b78c-163ee13faf85
# ╠═71e9a4b4-5dbb-4a67-b898-0cfffd90998d
# ╟─cc35db1a-715b-419a-b5a7-3f4a3f8faab4
# ╟─d5967d18-62d7-4870-b3b7-af7416acb2b4
# ╠═1ebe900d-2301-44ea-86b3-0fd7f03769b1
# ╠═a6d31357-bb6d-4fb2-a642-112f39d32832
# ╠═ac88d890-6514-4e82-b776-ff98467d8708
# ╠═241d663c-dc41-4039-b9a2-0f48d792b8a1
# ╠═2a111cc4-a80c-43d0-83b3-7c827761a4fa
# ╟─5d2801db-a987-48e5-b71d-6674c0f7af58
# ╟─ebc5cfbd-dd1d-4089-8798-e35357b38cac
# ╟─4e5c0d80-2871-4f57-a0a0-acbffb4ec199
# ╠═af3c116c-510d-4c5c-9b64-fb56a106e5eb
# ╟─f5d03c3a-3854-4c08-9cf7-074dd628393b
# ╠═fe698e99-939b-44d4-aaed-1bdfdc5c58f0
# ╠═b4f7a35f-0f74-465b-82a0-31c404b49a10
# ╠═7516d09f-7356-41ac-a565-2c4e2c21408d
# ╠═3e754434-5934-42df-861f-614bfa19643f
# ╠═00762ff1-c860-4d2c-b5ef-4a27b28b5d7f
# ╟─dccfe48e-4068-4f35-8446-8476033f90e6
# ╟─dec0718e-799d-48ce-b444-cf07daf803c6
# ╠═c36cb78e-09ce-4521-803f-3e94d5ebd7cd
# ╠═fa2bd5ca-1fc5-4706-86b1-c7baf5bd2a47
# ╠═4a90cefa-9793-481d-ae1b-863a379369ec
# ╠═dc676a84-430c-4316-9e42-e935d3dc9e74
# ╠═366feb2b-b08e-4176-9623-149b6cbaf989
# ╟─d0099447-0c7d-4573-b9b0-ff18398f36a4
# ╟─af57ef9e-0030-4b2d-8b1b-127f50d9754a
# ╟─ef982930-6388-466d-b0e1-8d715316a884
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using LinearAlgebra
using SparseArrays
using BandedMatrices

import Base: show

struct Basis
    t::AbstractKnotSet
    x::AbstractVector
    w::AbstractVector
    B::Vector{AbstractMatrix}
    ∂B::Vector{AbstractMatrix}
    bl
    br
end

# See http://pages.cs.wisc.edu/~deboor/pgs/bsplvb.f
function evaluate!(Bᵢ, t::AbstractKnotSet{T}, x::AbstractVector, bl::T, br::T) where T
    #=
    \[B_{i,1,t}(x)=\begin{cases}
    1, & t_i\leq x < t_{i+1}\\
    0, & \textrm{else}
    \end{cases}\]
    =#
    for i = 1:length(t)-1
        Bᵢ[1][:,i] = (t[i] .<= x .< t[i+1])
    end

    #=
    \[B_{i,k,t}(x)=\frac{x-t_i}{t_{i+k-1}-t_i}B_{i,k-1,t}(x)-
    \frac{t_{i+k}-x}{t_{i+k}-t_{i+1}}B_{i+1,k-1,t}(x)\]
    =#
    for kk = 2:order(t)
        for i = 1:length(t) - kk
            d₁ = t[i+kk-1]-t[i]
            d₂ = t[i+kk]-t[i+1]
            f₁ = d₁ == 0 ? 0 : (x .- t[i])/d₁
            f₂ = d₂ == 0 ? 0 : (t[i+kk] .- x)/d₂
            Bᵢ[kk][:,i] = f₁.*Bᵢ[kk-1][:,i] + f₂.*Bᵢ[kk-1][:,i+1]
        end
    end
    Bᵢ[end][:,1] *= bl
    Bᵢ[end][:,end] *= br
end

function evaluate!(∂Bᵢ, Bᵢ, t::BSplines.AbstractKnotSet{T}, x::AbstractVector,
                   bl::T, br::T) where T
    #=
    \[\partial B_{i,k,t}(x) = (k-1)
    \left[
    \frac{B_{i,k-1,t}(x)}{t_{i+k-1}-t_i}-
    \frac{B_{i+1,k-1,t}(x)}{t_{i+k}-t_{i+1}}
    \right]\]
    \[\begin{aligned}
    \partial^{(n)} B_{i,k,t}(x)
    &= (k-1)
    \left[
    \frac{\partial^{(n-1)} B_{i,k-1,t}(x)}{t_{i+k-1}-t_i}-
    \frac{\partial^{(n-1)} B_{i+1,k-1,t}(x)}{t_{i+k}-t_{i+1}}
    \right] \\
    &=
    (k-1)(k-2)
    \left[
    \frac{
    \frac{\partial^{(n-2)}B_{i,k-2,t}(x)}{t_{i+k-2}-t_i}-
    \frac{\partial^{(n-2)}B_{i+1,k-2,t}(x)}{t_{i+k-1}-t_{i+1}}
    }{t_{i+k-1}-t_i}-
    \frac{
    \frac{\partial^{(n-2)}B_{i+1,k-2,t}(x)}{t_{i+k-1}-t_{i+1}}-
    \frac{\partial^{(n-2)}B_{i+2,k-2,t}(x)}{t_{i+k}-t_{i+2}}
    }{t_{i+k}-t_{i+1}}
    \right]
    \end{aligned}\]

    de Boor X(15):
    \[ \partial^{(m)}\sum_j \alpha_j B_{j,k,t}
    =\sum_j \alpha_j^{(m+1)} B_{j,k-m,t}\]

    where [de Boor X(16)]

    \[\alpha_r^{(m+1)}=
    \begin{cases}
    \alpha_r, & m=0,\\
    \frac{\alpha_r^{(m)}-\alpha_{r-1}^{(m)}}{(t_{r+k-m}-t_r)/(k-m)}, & m>0.
    \end{cases}\]

    To calculate \(\partial^{(m)}B_{i,k,t}\), \(\forall m\in[1,k-1]\), we do

    - \(\forall i\in[1,N]\), set \(\alpha^{(1)}_j = \delta_{ij}\)
      - \(\forall m\in[1,k-1]\)
        - Compute \(\alpha^{(m+1)}\) using X(16)
        - Compute \(\partial^{(m)}B_{i,k,t}\) using X(15)
    =#

    evaluate!(Bᵢ, t, x, bl, br)

    k = order(t)
    N = length(t) - k
    α = Matrix{T}(undef, N, k)
    for i ∈ 1:N
        α[:,1] .= zero(T)
        α[i,1] = if i == 1
            bl
        elseif i == N
            br
        else
            one(T)
        end
        for m ∈ 1:k-1
            ∂Bᵢ[m][:,i] .= zero(T)
            for j ∈ 1:N
                δt = t[j+k-m] - t[j]
                if δt == 0
                    α[j,m+1] = 0
                    continue
                end
                α[j,m+1] = (k-m)*(α[j,m]-(j>1 ? α[j-1,m] : 0))/δt
                ∂Bᵢ[m][:,i] += α[j,m+1]*Bᵢ[k-m][:,j]
            end
        end
    end
end

function Basis(t::AbstractKnotSet{T}, bl::T=one(T), br::T=one(T)) where T
    x,w = lgwt(t)
    B = [spzeros(eltype(x), length(x), length(t)-kk)
         for kk = 1:order(t)]
    ∂B = [spzeros(eltype(x), length(x), length(t)-kk)
          for kk = 1:order(t)]
    evaluate!(∂B, B, t, x, bl, br)
    Basis(t, x, w, B, ∂B, bl, br)
end
Basis(t::AbstractKnotSet{T}, bl, br) where T = Basis(t, T(bl), T(br))

function (basis::Basis)(x::AbstractVector)
    Bᵢ = [spzeros(eltype(x), length(x), length(basis.t)-kk)
          for kk = 1:order(basis.t)]
    evaluate!(Bᵢ, basis.t, x, basis.bl, basis.br)
    Bᵢ
end

"""
    basis(f)

Generate the scalar operator corresponding to `f(x)` on the BSpline
`basis`.
"""
function (basis::Basis)(f::Function)
    m = size(basis.B[end],2)
    k = order(basis.t)
    B = BandedMatrix{eltype(basis.t)}(undef, m,m, 0, k-1)
    fwx = weights(basis) .* f.(locs(basis))
    for j = 1:m
        for i = max(1,j-k):j
            B[i,j] = dot(basis.B[end][:,i], (fwx .* basis.B[end][:,j]))
        end
    end
    Symmetric(B)
end

"""
    basis(I)

Return the overlap matrix of `basis`.
"""
(basis::Basis)(::UniformScaling) = basis(one)

function derop(basis::Basis, o::Integer)
    m = size(basis.B[end],2)
    k = order(basis.t)
    D = BandedMatrix{eltype(basis.t)}(undef, m,m, k-1, k-1)
    symmetric = o == 2 && basis.bl == 0 && basis.br == 0
    for j = 1:m
        for i = max(1,j-k):(symmetric ? j : min(j+k,m))
            D[i,j] = dot(basis.B[end][:,i], (weights(basis) .* basis.∂B[o][:,j]))
        end
    end
    symmetric ? Symmetric(D) : D
end

locs(basis::Basis) = basis.x
weights(basis::Basis) = basis.w

function show(io::IO, basis::Basis)
    write(io, "BSpline basis with $(basis.t)")
end

@recipe function plot(basis::Basis, n=100)
    x = range(first(basis.t), stop=last(basis.t),
              length=n*numintervals(basis.t))
    x, Matrix(basis(x)[end])
end

export locs, weights, derop

ext: jl


text:
function guides_position(fig::Figure, position)
    position = Symbol(position)
    supported = [:top, :bottom, :left, :right]
    position ∉ supported && throw(ArgumentError("Legend position $position ∉ $supported"))

    if position == :bottom
        legs_pos = fig[end+1,:]
    elseif position == :top
        legs_pos = fig[0,:]
    elseif position == :right
        legs_pos = fig[:,end+1]
    elseif position == :left
        legs_pos = fig[:,0]
    end

    legs_pos
end

default_orientation(position) = position in [:top, :bottom] ? :horizontal : :vertical

default_isvertical(position) = default_orientation(position) == :vertical


text:
using Setfield
using Setfield: PropertyLens, ComposedLens, IdentityLens, IndexLens, DynamicIndexLens

"""
    VarName{sym}(lens::Lens=IdentityLens())

A variable identifier for a symbol `sym` and lens `lens`.

The Julia variable in the model corresponding to `sym` can refer to a single value or to a
hierarchical array structure of univariate, multivariate or matrix variables. The field `lens`
stores the indices requires to access the random variable from the Julia variable indicated by `sym`
as a tuple of tuples. Each element of the tuple thereby contains the indices of one lens
operation.

`VarName`s can be manually constructed using the `VarName{sym}(lens)` constructor, or from an
lens expression through the [`@varname`](@ref) convenience macro.

# Examples

```jldoctest; setup=:(using Setfield)
julia> vn = VarName{:x}(Setfield.IndexLens((Colon(), 1)) ∘ Setfield.IndexLens((2, )))
x[:,1][2]

julia> getlens(vn)
(@lens _[Colon(), 1][2])

julia> @varname x[:, 1][1+1]
x[:,1][2]
```
"""
struct VarName{sym,T<:Lens}
    lens::T

    function VarName{sym}(lens=IdentityLens()) where {sym}
        # TODO: Should we completely disallow or just `@warn` of limited support?
        if !is_static_lens(lens)
            error("attempted to construct `VarName` with dynamic lens of type $(nameof(typeof(lens)))")
        end
        return new{sym,typeof(lens)}(lens)
    end
end

"""
    is_static_lens(l::Lens)

Return `true` if `l` does not require runtime information to be resolved.

In particular it returns `false` for `Setfield.DynamicLens` and `Setfield.FunctionLens`.
"""
is_static_lens(l::Lens) = is_static_lens(typeof(l))
is_static_lens(::Type{<:Lens}) = false
is_static_lens(::Type{<:Union{PropertyLens, IndexLens, IdentityLens}}) = true
function is_static_lens(::Type{ComposedLens{LO, LI}}) where {LO, LI}
    return is_static_lens(LO) && is_static_lens(LI)
end

# A bit of backwards compatibility.
VarName{sym}(indexing::Tuple) where {sym} = VarName{sym}(tupleindex2lens(indexing))

"""
    VarName(vn::VarName, lens::Lens)
    VarName(vn::VarName, indexing::Tuple)

Return a copy of `vn` with a new index `lens`/`indexing`.

```jldoctest; setup=:(using Setfield)
julia> VarName(@varname(x[1][2:3]), Setfield.IndexLens((2,)))
x[2]

julia> VarName(@varname(x[1][2:3]), ((2,),))
x[2]

julia> VarName(@varname(x[1][2:3]))
x
```
"""
VarName(vn::VarName, lens::Lens = IdentityLens()) = VarName{getsym(vn)}(lens)

function VarName(vn::VarName, indexing::Tuple)
    return VarName{getsym(vn)}(tupleindex2lens(indexing))
end

tupleindex2lens(indexing::Tuple{}) = IdentityLens()
tupleindex2lens(indexing::Tuple{<:Tuple}) = IndexLens(first(indexing))
function tupleindex2lens(indexing::Tuple)
    return IndexLens(first(indexing)) ∘ tupleindex2lens(indexing[2:end])
end

"""
    getsym(vn::VarName)

Return the symbol of the Julia variable used to generate `vn`.

## Examples

```jldoctest
julia> getsym(@varname(x[1][2:3]))
:x

julia> getsym(@varname(y))
:y
```
"""
getsym(vn::VarName{sym}) where {sym} = sym

"""
    getlens(vn::VarName)

Return the lens of the Julia variable used to generate `vn`.

## Examples

```jldoctest
julia> getlens(@varname(x[1][2:3]))
(@lens _[1][2:3])

julia> getlens(@varname(y))
(@lens _)
```
"""
getlens(vn::VarName) = vn.lens

@deprecate getindexing(vn::VarName) getlens(vn)

"""
    get(obj, vn::VarName{sym})

Alias for `get(obj, PropertyLens{sym}() ∘ getlens(vn))`.
"""
function Setfield.get(obj, vn::VarName{sym}) where {sym}
    return Setfield.get(obj, PropertyLens{sym}() ∘ getlens(vn))
end

"""
    set(obj, vn::VarName{sym}, value)

Alias for `set(obj, PropertyLens{sym}() ∘ getlens(vn), value)`.
"""
function Setfield.set(obj, vn::VarName{sym}, value) where {sym}
    return Setfield.set(obj, PropertyLens{sym}() ∘ getlens(vn), value)
end

Base.hash(vn::VarName, h::UInt) = hash((getsym(vn), getlens(vn)), h)
function Base.:(==)(x::VarName, y::VarName)
    return getsym(x) == getsym(y) && getlens(x) == getlens(y)
end

# Allow compositions with lenses.
function Base.:∘(vn::VarName{sym,<:Lens}, lens::Lens) where {sym}
    return VarName{sym}(getlens(vn) ∘ lens)
end

function Base.show(io::IO, vn::VarName{<:Any,<:Lens})
    # No need to check `Setfield.has_atlens_support` since
    # `VarName` does not allow dynamic lenses.
    print(io, getsym(vn))
    _print_application(io, getlens(vn))
end

# This is all just to allow to convert `Colon()` into `:`.
_print_application(io::IO, l::Lens) = Setfield.print_application(io, l)
function _print_application(io::IO, l::ComposedLens)
    _print_application(io, l.outer)
    _print_application(io, l.inner)
end
_print_application(io::IO, l::IndexLens) =
    print(io, "[", join(map(prettify_index, l.indices), ","), "]")
# This is a bit weird but whatever. We're almost always going to
# `concretize` anyways.
_print_application(io::IO, l::DynamicIndexLens) = print(io, l, "(_)")

prettify_index(x) = string(x)
prettify_index(::Colon) = ":"

"""
    Symbol(vn::VarName)

Return a `Symbol` represenation of the variable identifier `VarName`.

# Examples
```jldoctest
julia> Symbol(@varname(x[1][2:3]))
Symbol("x[1][2:3]")

julia> Symbol(@varname(x[1][:]))
Symbol("x[1][:]")
```
"""
Base.Symbol(vn::VarName) = Symbol(string(vn))  # simplified symbol

"""
    inspace(vn::Union{VarName, Symbol}, space::Tuple)

Check whether `vn`'s variable symbol is in `space`.  The empty tuple counts as the "universal space"
containing all variables. Subsumption (see [`subsume`](@ref)) is respected.

## Examples

```jldoctest
julia> inspace(@varname(x[1][2:3]), ())
true

julia> inspace(@varname(x[1][2:3]), (:x,))
true

julia> inspace(@varname(x[1][2:3]), (@varname(x),))
true

julia> inspace(@varname(x[1][2:3]), (@varname(x[1:10]), :y))
true

julia> inspace(@varname(x[1][2:3]), (@varname(x[:][2:4]), :y))
true

julia> inspace(@varname(x[1][2:3]), (@varname(x[1:10]),))
true
```
"""
inspace(vn, space::Tuple{}) = true # empty tuple is treated as universal space
inspace(vn, space::Tuple) = vn in space
inspace(vn::VarName, space::Tuple{}) = true
inspace(vn::VarName, space::Tuple) = any(_in(vn, s) for s in space)

_in(vn::VarName, s::Symbol) = getsym(vn) == s
_in(vn::VarName, s::VarName) = subsumes(s, vn)

"""
    subsumes(u::VarName, v::VarName)

Check whether the variable name `v` describes a sub-range of the variable `u`.  Supported
indexing:

  - Scalar:

  ```jldoctest
  julia> subsumes(@varname(x), @varname(x[1, 2]))
  true

  julia> subsumes(@varname(x[1, 2]), @varname(x[1, 2][3]))
  true
  ```

  - Array of scalar: basically everything that fulfills `issubset`.

  ```jldoctest
  julia> subsumes(@varname(x[[1, 2], 3]), @varname(x[1, 3]))
  true

  julia> subsumes(@varname(x[1:3]), @varname(x[2][1]))
  true
  ```

  - Slices:

  ```jldoctest
  julia> subsumes(@varname(x[2, :]), @varname(x[2, 10][1]))
  true
  ```

Currently _not_ supported are:

  - Boolean indexing, literal `CartesianIndex` (these could be added, though)
  - Linear indexing of multidimensional arrays: `x[4]` does not subsume `x[2, 2]` for a matrix `x`
  - Trailing ones: `x[2, 1]` does not subsume `x[2]` for a vector `x`
"""
function subsumes(u::VarName, v::VarName)
    return getsym(u) == getsym(v) && subsumes(u.lens, v.lens)
end

# Idea behind `subsumes` for `Lens` is that we traverse the two lenses in parallel,
# checking `subsumes` for every level. This for example means that if we are comparing
# `PropertyLens{:a}` and `PropertyLens{:b}` we immediately know that they do not subsume
# each other since at the same level/depth they access different properties.
# E.g. `x`, `x[1]`, i.e. `u` is always subsumed by `t`
subsumes(t::IdentityLens, u::Lens) = true
subsumes(t::Lens, u::IdentityLens) = false

subsumes(t::ComposedLens, u::ComposedLens) =
    subsumes(t.outer, u.outer) && subsumes(t.inner, u.inner)

# If `t` is still a composed lens, then there is no way it can subsume `u` since `u` is a
# leaf of the "lens-tree".
subsumes(t::ComposedLens, u::PropertyLens) = false
# Here we need to check if `u.outer` (i.e. the next lens to be applied from `u`) is
# subsumed by `t`, since this would mean that the rest of the composition is also subsumed
# by `t`.
subsumes(t::PropertyLens, u::ComposedLens) = subsumes(t, u.outer)

# For `PropertyLens` either they have the same `name` and thus they are indeed the same.
subsumes(t::PropertyLens{name}, u::PropertyLens{name}) where {name} = true
# Otherwise they represent different properties, and thus are not the same.
subsumes(t::PropertyLens, u::PropertyLens) = false

# Indices subsumes if they are subindices, i.e. we just call `_issubindex`.
# FIXME: Does not support `DynamicIndexLens`.
# FIXME: Does not correctly handle cases such as `subsumes(x, x[:])`
#        (but neither did old implementation).
subsumes(t::IndexLens, u::IndexLens) = _issubindex(t.indices, u.indices)
subsumes(t::ComposedLens{<:IndexLens}, u::ComposedLens{<:IndexLens}) = subsumes_index(t, u)
subsumes(t::IndexLens, u::ComposedLens{<:IndexLens}) = subsumes_index(t, u)
subsumes(t::ComposedLens{<:IndexLens}, u::IndexLens) = subsumes_index(t, u)

# Since expressions such as `x[:][:][:][1]` and `x[1]` are equal,
# the indexing behavior must be considered jointly.
# Therefore we must recurse until we reach something that is NOT
# indexing, and then consider the sequence of indices leading up to this.
"""
    subsumes_index(t::Lens, u::Lens)

Return `true` if the indexing represented by `t` subsumes `u`.

This is mostly useful for comparing compositions involving `IndexLens`
e.g. `_[1][2].a[2]` and `_[1][2].a`. In such a scenario we do the following:
1. Combine `[1][2]` into a `Tuple` of indices using [`combine_indices`](@ref).
2. Do the same for `[1][2]`.
3. Compare the two tuples from (1) and (2) using `subsumes_index`.
4. Since we're still undecided, we call `subsume(@lens(_.a[2]), @lens(_.a))`
   which then returns `false`.

# Example
```jldoctest; setup=:(using Setfield; using AbstractPPL: subsumes_index)
julia> t = @lens(_[1].a); u = @lens(_[1]);

julia> subsumes_index(t, u)
false

julia> subsumes_index(u, t)
true

julia> # `IdentityLens` subsumes all.
       subsumes_index(@lens(_), t)
true

julia> # None subsumes `IdentityLens`.
       subsumes_index(t, @lens(_))
false

julia> AbstractPPL.subsumes(@lens(_[1][2].a[2]), @lens(_[1][2].a))
false

julia> AbstractPPL.subsumes(@lens(_[1][2].a), @lens(_[1][2].a[2]))
true
```
"""
function subsumes_index(t::Lens, u::Lens)
    t_indices, t_next = combine_indices(t)
    u_indices, u_next = combine_indices(u)

    # If we already know that `u` is not subsumed by `t`, return early.
    if !subsumes_index(t_indices, u_indices)
        return false
    end

    if t_next === nothing
        # Means that there's nothing left for `t` and either nothing
        # or something left for `u`, i.e. `t` indeed `subsumes` `u`.
        return true
    elseif u_next === nothing
        # If `t_next` is not `nothing` but `u_ntext` is, then
        # `t` does not subsume `u`.
        return false
    end

    # If neither is `nothing` we continue.
    return subsumes(t_next, u_next)
end

"""
    combine_indices(lens)

Return sequential indexing into a single `Tuple` of indices,
e.g. `x[:][1][2]` becomes `((Colon(), ), (1, ), (2, ))`.

The result is compatible with [`subsumes_index`](@ref) for `Tuple` input.
"""
combine_indices(lens::Lens) = (), lens
combine_indices(lens::IndexLens) = (lens.indices,), nothing
function combine_indices(lens::ComposedLens{<:IndexLens})
    indices, next = combine_indices(lens.inner)
    return (lens.outer.indices, indices...), next
end

"""
    subsumes_index(left_index::Tuple, right_index::Tuple)

Return `true` if `right_index` is subsumed by `left_index`.

Currently _not_ supported are:
- Boolean indexing, literal `CartesianIndex` (these could be added, though)
- Linear indexing of multidimensional arrays: `x[4]` does not subsume `x[2, 2]` for a matrix `x`
- Trailing ones: `x[2, 1]` does not subsume `x[2]` for a vector `x`
- Dynamic indexing, e.g. `x[1]` does not subsume `x[begin]`.
"""
subsumes_index(::Tuple{}, ::Tuple{}) = true  # x subsumes x
subsumes_index(::Tuple{}, ::Tuple) = true    # x subsumes x[1]
subsumes_index(::Tuple, ::Tuple{}) = false   # x[1] does not subsume x
function subsumes_index(t::Tuple, u::Tuple)  # does x[i]... subsume x[j]...?
    return _issubindex(first(t), first(u)) && subsumes_index(Base.tail(t), Base.tail(u))
end

const AnyIndex = Union{Int,AbstractVector{Int},Colon}
_issubindex_(::Tuple{Vararg{AnyIndex}}, ::Tuple{Vararg{AnyIndex}}) = false
function _issubindex(t::NTuple{N,AnyIndex}, u::NTuple{N,AnyIndex}) where {N}
    return all(_issubrange(j, i) for (i, j) in zip(t, u))
end

const ConcreteIndex = Union{Int,AbstractVector{Int}} # this include all kinds of ranges

"""Determine whether indices `i` are contained in `j`, treating `:` as universal set."""
_issubrange(i::ConcreteIndex, j::ConcreteIndex) = issubset(i, j)
_issubrange(i::Colon, j::Colon) = true
_issubrange(i::ConcreteIndex, j::Colon) = true
# FIXME: [2021-07-31] This is wrong but we have tests in DPPL that tell
# us that it SHOULD be correct. I'll leave it as is for now to ensure that
# we preserve the status quo, but I'm confused.
_issubrange(i::Colon, j::ConcreteIndex) = true

"""
    concretize(l::Lens, x)

Return `l` instantiated on `x`, i.e. any runtime information evaluated using `x`.
"""
concretize(I::Lens, x) = I
concretize(I::DynamicIndexLens, x) = IndexLens(I.f(x))
function concretize(I::ComposedLens, x)
    x_inner = get(x, I.outer)
    return ComposedLens(concretize(I.outer, x), concretize(I.inner, x_inner))
end
"""
    concretize(vn::VarName, x)

Return `vn` instantiated on `x`, i.e. any runtime information evaluated using `x`.

# Examples
```jldoctest; setup=:(using Setfield)
julia> x = (a = [1.0 2.0;], );

julia> AbstractPPL.concretize(@varname(x.a[1, end][:]), x)
x.a[1,2][:]
```
"""
concretize(vn::VarName, x) = VarName(vn, concretize(getlens(vn), x))

"""
    @varname(expr)

A macro that returns an instance of [`VarName`](@ref) given a symbol or indexing expression `expr`.

If `concretize` is `true`, the resulting expression will be wrapped in a [`concretize`](@ref) call.

Note that expressions involving dynamic indexing, i.e. `begin` and/or `end`, will need to be
resolved as `VarName` only supports non-dynamic indexing as determined by
[`is_static_index`](@ref). See examples below.

## Examples
### Dynamic indexing
```jldoctest
julia> # Dynamic indexing is not allowed in `VarName`
       @varname(x[end])
ERROR: UndefVarError: x not defined
[...]

julia> # To be able to resolve `end` we need `x` to be available.
       x = randn(2); @varname(x[end])
x[2]

julia> # Note that "dynamic" here refers to usage of `begin` and/or `end`,
       # _not_ "information only available at runtime", i.e. the following works.
       [@varname(x[i]) for i = 1:length(x)][end]
x[2]
```

### General indexing

Under the hood Setfield.jl's `Lens` are used for the indexing:

```jldoctest
julia> getlens(@varname(x))
(@lens _)

julia> getlens(@varname(x[1]))
(@lens _[1])

julia> getlens(@varname(x[:, 1]))
(@lens _[Colon(), 1])

julia> getlens(@varname(x[:, 1][2]))
(@lens _[Colon(), 1][2])

julia> getlens(@varname(x[1,2][1+5][45][3]))
(@lens _[1, 2][6][45][3])
```

This also means that we support property access:

```jldoctest
julia> getlens(@varname(x.a))
(@lens _.a)

julia> getlens(@varname(x.a[1]))
(@lens _.a[1])

julia> x = (a = [(b = rand(2), )], ); getlens(@varname(x.a[1].b[end]))
(@lens _.a[1].b[2])
```

!!! compat "Julia 1.5"
    Using `begin` in an indexing expression to refer to the first index requires at least
    Julia 1.5.
"""
macro varname(expr::Union{Expr,Symbol})
    return varname(expr)
end

varname(sym::Symbol) = :($(AbstractPPL.VarName){$(QuoteNode(sym))}())
function varname(expr::Expr)
    if Meta.isexpr(expr, :ref) || Meta.isexpr(expr, :.)
        # Split into object/base symbol and lens.
        sym_escaped, lens = Setfield.parse_obj_lens(expr)
        # Setfield.jl escapes the return symbol, so we need to unescape
        # to call `QuoteNode` on it.
        sym = drop_escape(sym_escaped)

        return if Setfield.need_dynamic_lens(expr)
            :(
                $(AbstractPPL.VarName){$(QuoteNode(sym))}(
                    $(AbstractPPL.concretize)($lens, $sym_escaped)
                )
            )
        else
            :($(AbstractPPL.VarName){$(QuoteNode(sym))}($lens))
        end
    else
        error("Malformed variable name $(expr)!")
    end
end

drop_escape(x) = x
function drop_escape(expr::Expr)
    Meta.isexpr(expr, :escape) && return drop_escape(expr.args[1])
    return Expr(expr.head, map(x -> drop_escape(x), expr.args)...)
end

"""
    @vsym(expr)

A macro that returns the variable symbol given the input variable expression `expr`.
For example, `@vsym x[1]` returns `:x`.

## Examples

```jldoctest
julia> @vsym x
:x

julia> @vsym x[1,1][2,3]
:x

julia> @vsym x[end]
:x
```
"""
macro vsym(expr::Union{Expr,Symbol})
    return QuoteNode(vsym(expr))
end

"""
    vsym(expr)

Return name part of the [`@varname`](@ref)-compatible expression `expr` as a symbol for input of the
[`VarName`](@ref) constructor."""
function vsym end

vsym(expr::Symbol) = expr
function vsym(expr::Expr)
    if Meta.isexpr(expr, :ref) || Meta.isexpr(expr, :.)
        return vsym(expr.args[1])
    else
        error("Malformed variable name $(expr)!")
    end
end


text:
module SolovevAbstract

    import ..ElectromagneticFields
    import ..ElectromagneticFields: AnalyticEquilibrium, AnalyticPerturbation

    export AbstractSolovevEquilibrium

    abstract type AbstractSolovevEquilibrium <: AnalyticEquilibrium end

    R(x::AbstractVector, equ::AbstractSolovevEquilibrium) = x[1] * equ.R₀
    Z(x::AbstractVector, equ::AbstractSolovevEquilibrium) = x[2] * equ.R₀
    ϕ(x::AbstractVector, equ::AbstractSolovevEquilibrium) = x[3]

    X(x::AbstractVector, equ::AbstractSolovevEquilibrium) = R(x,equ) * cos(ϕ(x,equ))
    Y(x::AbstractVector, equ::AbstractSolovevEquilibrium) = R(x,equ) * sin(ϕ(x,equ))
    θ(x::AbstractVector, equ::AbstractSolovevEquilibrium) = atan(Z(x,equ), R(x,equ) - equ.R₀)

    r²(x::AbstractVector, equ::AbstractSolovevEquilibrium) = (R(x,equ) - equ.R₀)^2 + Z(x,equ)^2
    r(x::AbstractVector, equ::AbstractSolovevEquilibrium) = sqrt(r²(x, equ))

    ElectromagneticFields.J(x::AbstractVector, equ::AbstractSolovevEquilibrium) = R(x,equ) * equ.R₀^2

    ElectromagneticFields.A₁(x::AbstractVector, equ::AbstractSolovevEquilibrium) = + equ.B₀ * equ.R₀ * x[2] / x[1] / 2
    ElectromagneticFields.A₂(x::AbstractVector, equ::AbstractSolovevEquilibrium) = - equ.B₀ * equ.R₀ * log(x[1]) / 2

    ElectromagneticFields.x¹(ξ::AbstractVector, equ::AbstractSolovevEquilibrium) = X(ξ,equ)
    ElectromagneticFields.x²(ξ::AbstractVector, equ::AbstractSolovevEquilibrium) = Y(ξ,equ)
    ElectromagneticFields.x³(ξ::AbstractVector, equ::AbstractSolovevEquilibrium) = Z(ξ,equ)

    ElectromagneticFields.ξ¹(x::AbstractVector, equ::AbstractSolovevEquilibrium) = sqrt(x[1]^2 + x[2]^2) / equ.R₀
    ElectromagneticFields.ξ²(x::AbstractVector, equ::AbstractSolovevEquilibrium) = x[3] / equ.R₀
    ElectromagneticFields.ξ³(x::AbstractVector, equ::AbstractSolovevEquilibrium) = atan(x[2], x[1])

    ElectromagneticFields.g₁₁(x::AbstractVector, equ::AbstractSolovevEquilibrium) = equ.R₀^2
    ElectromagneticFields.g₂₂(x::AbstractVector, equ::AbstractSolovevEquilibrium) = equ.R₀^2
    ElectromagneticFields.g₃₃(x::AbstractVector, equ::AbstractSolovevEquilibrium) = R(x, equ)^2

    ElectromagneticFields.get_functions(::AbstractSolovevEquilibrium) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ, r²=r²)
    ElectromagneticFields.get_parameters(::AbstractSolovevEquilibrium) = (:R₀, :B₀, :ϵ, :κ, :δ, :α)

    function ElectromagneticFields.periodicity(x::AbstractArray{T,1}, ::AbstractSolovevEquilibrium) where {T <: Number}
        p = zero(x)
        p[3] = 2π
        return p
    end

end


text:
using BrainWave, Test

x = reshape(collect(1:27), 3,3,3)
@test isequal(deleteSlice3D(x, 2, 1), x[[1,3],:,:])
@test isequal(deleteSlice3D(x, [1,3], 2), x[:,[2], :])
@test isequal(deleteSlice3D(x, [3,1], 2), x[:,[2], :])
@test isequal(deleteSlice3D(x, [2,3], 3), x[:,:, [1]])


text:
# p17.jl - Helmholtz eq. u_xx + u_yy + (k^2)u = f
#         on [-1,1]x[-1,1]    (compare p16.jl)

# Set up spectral grid and tensor product Helmholtz operator:
N = 24; (D,x) = cheb(N); y = x;
xx = x[2:N]; yy = y[2:N];
f = @. exp(-10*((yy-1)^2+(xx'-.5)^2));
D2 = D^2; D2 = D2[2:N,2:N]; I = eye(N-1);
k = 9;
L = kron(I,D2) + kron(D2,I) + k^2*eye((N-1)^2);

# Solve for u, reshape to 2D grid, and plot:
u = L\f[:];
uu = zeros(N+1,N+1); uu[N:-1:2,N:-1:2] = reshape(u,N-1,N-1);
xxx = yyy = -1:.0333:1;
s = Spline2D(x[end:-1:1],y[end:-1:1],uu);
uuu = evalgrid(s,xxx,yyy);
figure(1); clf(); surf(xxx,yyy,uuu,rstride=1,cstride=1);
xlabel("x"); ylabel("y"); zlabel("u"); view(-37.5,30);
value = signif(uu[Int(N/2+1),Int(N/2+1)],10);
text3D(.2,1,.022,"u(0,0) = $value");
figure(2); clf(); contour(xxx,yyy,uuu,10); axis("square");


text:
divodd(a, b) = let u = a ÷ b
    if u < 3
        3
    elseif iseven(u)
        u + 1
    else
        u
    end
end

function mwmean(mean::Function, xs::Matrix, mw::Int)
    n, m = size(xs)
    result = Matrix{eltype(xs)}(undef, n-mw+1, m)
    for i in 1:size(xs, 1)-mw+1
        result[i,:] = mean(xs[(i-1).+(1:mw), :]; dims=1)
    end
    result
end

mwmean(xs::Matrix, mw::Int) = mwmean(mean, xs, mw)

function bin(frames::Array{Float64, 3}; mw=divodd(size(frames, 3), 10))
    n, m, t = size(frames)
    binned = bin(squish(frames); mw=mw)
    reshape(transpose(binned), n, m, t - mw + 1)
end

function bin(frames::Matrix{Float64}; mw=divodd(size(frames, 1), 10))
    bound = abs.(mwmean(mean, frames, mw)) .+ mwmean(std, frames, mw)
    binned = Int.(bound .< abs.(frames[(1+mw÷2):size(frames, 1)-mw÷2,:])) .+ 1
end

max_stars_repo_name: ELIFE-ASU/dglmoore

max_issues_repo_name: ELIFE-ASU/dglmoore

max_forks_repo_name: ELIFE-ASU/dglmoore


text:
using Test
using PETScX
using LinearAlgebra: norm
using MPI: mpiexec

@testset "vec tests" begin
    for petsclib in PETScX.petsclibs
        PETScX.Initialize(petsclib)
        PetscScalar = PETScX.scalartype(petsclib)
        PetscInt = PETScX.inttype(petsclib)

        julia_x = [PetscScalar(i) for i in 1:10]
        petsc_x = PETScX.Vec(petsclib, julia_x)

        # Check the basics
        @test length(petsc_x) == length(julia_x)
        @test norm(petsc_x) == norm(julia_x)
        #XXX: For some reason non-l2 is broken...
        # @test norm(petsc_x, PETScX.NORM_INFINITY) == norm(julia_x, Inf)
        display(petsc_x)

        # Check the view
        _stdout = stdout
        (rd, wr) = redirect_stdout()
        @show petsc_x
        @test readline(rd) == "petsc_x = Vec Object: 1 MPI processes"
        @test readline(rd) == "  type: seq"
        for x in julia_x
            @test readline(rd) == "$(Int(x))."
        end
        redirect_stdout(_stdout)

        # Test reading and writing values
        for (i, x) in enumerate(julia_x)
            @test petsc_x[i] === x
            petsc_x[i] ^= 2
            @test petsc_x[i] === x^2
        end

        # Test writing many values
        i = PetscInt.([0, 2, 4])
        v = PetscScalar.([-2, 4, 10])
        PETScX.setvalues!(petsc_x, i, v, PETScX.INSERT_VALUES)
        @test julia_x[i .+ 1] == v
        PETScX.setvalues!(petsc_x, i, v, PETScX.ADD_VALUES)
        @test julia_x[i .+ 1] == 2v

        # Test reading many values
        vals = fill!(similar(v), 0)
        PETScX.getvalues!(vals, petsc_x, i)
        @test vals == 2v

        PETScX.Finalize(petsclib)
    end
end


text:
#tests for "nuisance" inversions with HiQGA.transD_GP

## make options for a purely stationary properties GP
nmin, nmax = 2, 200
pnorm = 2.
λx, λy = 1, 1
x = 0:(0.005λx):λx
y = 0:(0.005λy):0.5λy
nmin, nmax = 2, 200
fbounds = [-2. 2]
δ = 0.1
sdev_prop = [0.1]
sdev_pos = [0.05, 0.05]
K = transD_GP.GP.Mat32()
λ = [0.05maximum(x), 0.05maximum(y)]
demean = false
sampledc = true
##
xall = zeros(2,length(x)*length(y))
for i in 1:size(xall,2)
    xid, yid = Tuple(CartesianIndices((length(x),length(y)))[i])
    xall[:,i] = [x[xid]; y[yid]]
end
xbounds = zeros(Float64,size(xall, 1), 2)
for dim in 1:size(xall, 1)
    xbounds[dim,:] = [minimum(xall[dim,:]), maximum(xall[dim,:])]
end

opt = transD_GP.OptionsStat(nmin = nmin,
                        nmax = nmax,
                        xbounds = xbounds,
                        fbounds = fbounds,
                        xall = xall,
                        λ = λ,
                        δ = δ,
                        demean = demean,
                        sampledc = sampledc,
                        sdev_prop = sdev_prop,
                        sdev_pos = sdev_pos,
                        pnorm = pnorm,
                        needλ²fromlog = false,
                        updatenonstat = false,
                        quasimultid = false,
                        K = K,
                        timesλ = 3,
                        )

@time m = transD_GP.init(opt)

@testset "Creation tests for OptionsNuisance" begin
    nu_sdev = [0.1, 0.2, 0.3]
    nu_bounds = [0. 4.; -1.5 1.0; -1. 0.]
    optn = transD_GP.OptionsNuisance(opt;
        sdev = nu_sdev,
        bounds = nu_bounds,
        updatenuisances = true,
    )
    @testset "Test parameters set properly" begin
        @test isapprox(optn.sdev, [0.4, 0.5, 0.3], rtol=1e-4) #sdev in the struct is in "absolute units"
        @test all(optn.bounds .== nu_bounds)
        @test isapprox(optn.rotatebounds, nu_bounds .- mean(nu_bounds, dims=2), rtol=1e-4) #no covariance matrix was supplied
        @test optn.nnu == 3
        @test optn.updatenuisances == true
    end
    @testset "Constructor does not modify arguments" begin
        @test all(nu_sdev .== [0.1, 0.2, 0.3])
        @test all(nu_bounds .== [0. 4.; -1.5 1.0; -1. 0.])
    end
end

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA


text:
import BioMetaheuristics.TestFunctions: Benchmark, evaluate

"""
    PSO <: PopulationBase

`PSO` is the type associated with the implementation for the
Particle Swarm Optimization with momentum. See [`Algorithms`](@ref implementations-docs)
for more information.
"""
struct PSO <: PopulationBase end

# To enable dispatch based on the type
function optimize(f, population, k_max, rng, ::PSO; kwargs...)
    return PSO(f, population, k_max, rng; kwargs...)
end

"""
    PSO(f::Function, population, k_max::Int, rng;
        w=0.9, c1=2.0, c2=2.0
    ) -> OptimizationResults
    PSO(f::Benchmark, population, k_max::Int, rng;
        w=0.9, c1=2.0, c2=2.0
    ) -> OptimizationResults

Method that implements `PSO` for a function `f` of type `Function`
or of type `Benchmark`.

Returns an `OptimizationResults` type with information relevant to the
run executed, see [`OptimizationResults`](@ref).

# Arguments

- `population`: can be any `AbstractArray` that contains [`Particle`](@ref)
instances, but it is expected to be generated by [`Population`](@ref).
- `k_max`: number of maximum iterations until "convergence" of the algorithm.
- `rng`: An object of type `AbstractRNG`.

# Keyword arguments

_It is recommended to use the default values provided._

- `w`: value that controls how much of the initial velocity is retained, i.e. an inertia term.
This value decays linearly over each iteration until it reaches
the default miminum value of 0.4.
- `c1`: balance between the influence of the individual's knowledge, i.e. the best inidividual solution so far.
- `c2`: balance between the influence of the population's knowledge, i.e. the best global solution so far.
- `seed`: an integer to be used as the seed for the pseudo random number generators.
If `nothing` is passed (the default), then a random seed will be taken from the system.

# Examples
```julia
using BioMetaheuristics
using Random

rng = MersenneTwister(1261234)
# Define the Sphere function
f_sphere(x) = sum(x.^2)

# Implement PSO for a 3-dimensional Sphere function, with
# 10000 iterations and 30 particles in the population.
val = PSO(f_sphere, Population(30, 3, -15.0, 15.0), 10_000, rng)
```
"""
function PSO(f::Function, population, k_max::Int, rng; w=0.9, c1=2.0, c2=2.0)

    val = _pso!(f, population, k_max, rng; w=w, c1=c1, c2=c2)

    optim_res = OptimizationResults(val, f(val), "PSO", k_max)

    return optim_res
end

function PSO(f::Benchmark, population, k_max::Int, rng; w=0.9, c1=2.0, c2=2.0)
    val = _pso!(
        x -> evaluate(f, x),
        population,
        k_max,
        rng;
        w=w,
        c1=c1,
        c2=c2
    )

    optim_res = OptimizationResults(val, evaluate(f, val), "PSO", k_max)

    return optim_res
end

function _pso!(f, population, k_max::Int, rng; w=0.9, c1=2.0, c2=2.0)
    # Obtain weight decay rate
    η = _weight_decay(w, k_max)

    # Evaluate initial costs
    dimension = length(population[1].x)

    x_best = similar(population[1].x_best) # Best candidate
    y_best = Inf # Best miminum
    @inbounds for particle in population
        y = f(particle.x)
        if y < y_best
            copyto!(x_best, particle.x)
            y_best = y
        end
    end

    # PSO main loop
    @inbounds for _ in 1:k_max
        _update!(f, population, w, c1, c2, dimension, x_best, y_best, rng)
        # Make the inertia weight decay over time
        w -= η
    end

    # Sort the solutions and return the first one
    return population[1].x_best
end

function _update!(f, population, w, c1, c2, n, x_best, y_best, rng)
    # Create container for random numbers
    rngs = zeros(2, n)

    for P in population
        rand!(rng, rngs)
        # Evaluate velocity
        P.v = (w * P.v) + (c1 * rngs[1, :] .* (P.x_best .- P.x))
        P.v .+= c2 * rngs[2, :] .* (x_best .- P.x)
        # Update position
        P.x .+= P.v
        # Apply boundary values to positions and velocities
        _clip_positions_velocities!(P)
        # Update values if they give lower cost
        y = f(P.x)
        if y < y_best
            copyto!(x_best, P.x)
            y_best = y
        end
        if y < f(P.x_best)
            copyto!(P.x_best, P.x)
        end
    end

    return nothing
end

"""
    Compute the corresponding weight decay depending the maximum
number of iterations and the initial value for it.
"""
@inline function _weight_decay(initial, itr_max)
    # Following the references, the minimum is 0.4
    stop = 0.4
    step_size = (initial - stop) / itr_max

    return step_size
end