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

using Markdown
using InteractiveUtils

# ╔═╡ 0150a3e2-f98e-4412-9e96-1a2db5a9421e
using Combinatorics, Plots

# ╔═╡ 936344da-213a-11ec-3e68-05f0f85cf2f3
md"""
# NP-complete problems

*STMO*

**Michiel Stock**
"""

# ╔═╡ c3b61d51-679c-4aba-957b-105c00037aee
md"In this document, we will go over some simple algorithms to solve the knapsack problem."

# ╔═╡ 62c4926d-0350-4088-8729-0ee5afb306be
md"""
## The knapsack problem

### Problem definition

> Given a set of items, each with a weight and a value, find the subset of items so that the total weight is less than or equal to a given limit and the total value is as large as possible.

Or, in symbols:

Given two sets of cardinality $n$ with the values $\{v_1,\ldots, v_n\}$ and the weights $\{w_1,\ldots, w_n\}$ and a capacity $W>0$ we want to determine $T\subseteq\{1,\ldots,n\}$ such that

$$\max_T \, \sum_{i\in T} v_i$$

$$\text{subject to } \sum_{i\in T} w_i \leq W\,.$$

We can represent this problem in Julia using a list of tuples, representing the items with their reprective values and weights, and a number containing the capacity $W$. Let us however use this opportunity to construct a simple structure to represent knapsack problems.
"""

# ╔═╡ 3c769076-fec9-4730-b993-2fefa09f0c8d
struct Knapsack{T,V<:Real,W<:Real}
  items::Array{Tuple{T,V,W},1}
  capacity::W
end

# ╔═╡ ada8c7c4-f6d8-41e7-a4e0-c137ce7ab151
C = 10

# ╔═╡ 64e7be40-cd0b-40c4-b476-51c6a66d40e1
md"""
Consider the Indiana Jones problem:

| i |  artifact     | $v_i$ |  $w_i$ |
|---|---------------|-------|--------|
| 1 |  statue 1     |   1   |  2     |
| 2 |  statue 2     |   1   |  2     |
| 3 |  statue 3     |   1   |  2     |
| 4 |  tablet 1     |  10   |  5     |
| 5 |  tablet 2     |  10   |  5     |
| 6 |  golden mask  |  13   |  7     |
| 7 |  golden plate |   7   |  3     |
"""

# ╔═╡ be803662-aa08-4ee3-a400-bbc18b7d060f
md"""The capacity is $C kg.

So we can define the instance as (feel free to modify)"""

# ╔═╡ 2ce221f2-5712-440e-9b34-168a3a5af65e
knapsack = Knapsack([("statue 1", 1, 2),
                     ("statue 2", 1, 2),
                     ("statue 2", 1, 2),
                     ("tablet 1", 10, 5),
                     ("tablet 2", 10, 5),
                     ("golden mask", 13, 7),
                     ("golden plate", 7, 3)],
                     C)

# ╔═╡ 2d9598de-dfd6-46c7-bf4f-c2f48134e558
md"We might as well define some useful function that work on `Knapsack` structures."

# ╔═╡ 265ac424-3e7b-4222-ab79-07907544a9fa
"""Getter for the items."""
items(knapsack::Knapsack) = knapsack.items

# ╔═╡ 876bd16d-89ac-4164-8055-e9ec4a4d875f
"""Getter for the capacity."""
capacity(knapsack::Knapsack) = knapsack.capacity

# ╔═╡ cca6890d-b0eb-4154-9b27-be79c0b8ab38
"""Number of items"""
Base.length(knapsack::Knapsack) = length(items(knapsack))

# ╔═╡ 94cb9a4c-aa3f-4ef0-a4e5-170b0d2ad69b
"""Weight of a set of items"""
weight(items::Array) = sum((w for (i, v, w) in items))

# ╔═╡ 1d07e4af-5419-464f-842d-3012920e8a1b
"""Weight of an item"""
weight(item::Tuple) = item[3]

# ╔═╡ f0734b83-dbbc-4d1d-a96e-9209e53b7382
"""Value of a set of items"""
value(items::Array) = sum((v for (i, v, w) in items))

# ╔═╡ 165d0c6a-ec73-46a9-a9f4-5085c65d184c
"""Value of an item"""
value(item::Tuple) = item[2]

# ╔═╡ e1fa7a85-5aab-4e17-b02f-c9a0c62f3f47
md"### Brute force"

# ╔═╡ 38e36108-85cf-4808-8e10-bc6f869522fc
function bruteforce(knapsack::Knapsack{T,V,W}) where {T,V<:Real,W<:Real}
  # number of items
  n = length(knapsack)
  cap = capacity(knapsack)
  knapsack_items = items(knapsack)
  best_value = zero(W)
  best_solution = eltype(knapsack_items)[]
  for solution in combinations(knapsack_items)
    if weight(solution) ≤ cap && value(solution) > best_value
      best_solution = solution
      best_value = value(solution)
    end
  end
  return best_solution, best_value
end

# ╔═╡ 9b555941-5ef2-48d5-8efd-0f4cf7e38af6
bruteforce(knapsack)

# ╔═╡ 96d84e8c-3b11-4b01-a036-e3434804f512
md"### Greedy"

# ╔═╡ f65c49f0-09bf-4b5f-99b9-21b2fe6f5760
function greedy(knapsack::Knapsack{T,V,W},
                    heuristic::Function) where {T,V<:Real,W<:Real}
  items_knapsack = items(knapsack)
  solution = eltype(items_knapsack)[]
  solution_weight = zero(W)
  solution_value = zero(V)
  cap = capacity(knapsack)
  for item in sort(items_knapsack, by=heuristic, rev=true)
    v, w = value(item), weight(item)
    if solution_weight + w ≤ cap
      push!(solution, item)
      solution_weight += w
      solution_value += v
    end
  end
  return solution, solution_value
end

# ╔═╡ 3e95cf39-5d8f-498a-bdb6-df99a51257f3
md"Here, `heuristic` is a function we can provide to guide the greedy search. For example, if we search by value."

# ╔═╡ 97465212-39d3-4b19-8dd0-cf654b80e3e2
greedy(knapsack, value)

# ╔═╡ 834e4c69-5e90-409a-a2ed-ac14609ae728
md"Or, search lightest items first."

# ╔═╡ 02e07bf0-741c-4df2-9a07-55093e2bb55d
greedy(knapsack, item -> -weight(item))

# ╔═╡ e53659a3-e791-4318-a738-678d8c267b86
md"Or, by value-density:"

# ╔═╡ fa75a27d-725f-49df-8527-125c90de1692
greedy(knapsack, item -> value(item) / weight(item))

# ╔═╡ 79be9317-762d-4ce5-85ed-983f4882fa0b
md"### Dynamics programming"

# ╔═╡ d500504b-5cfb-46b9-8b82-1fbfb843b479
function dynamicprogramming(knapsack::Knapsack{T,V,W}) where {T,V<:Real,W<:Int}
  items_knapsack = items(knapsack)
  n = length(knapsack)
  cap = capacity(knapsack)
  DP = zeros(V, n + 1, cap + 1)  # starts from zero
  # fill the DP matrix
  for w in 1:cap
    for (i, (itname, vᵢ, wᵢ)) in enumerate(items_knapsack)
      if wᵢ ≤ w  # this item can fit the bag
        DP[i+1,w+1] = max(DP[i,w+1],  # value without the item
                        DP[i,w-wᵢ+1] + vᵢ)  # value with the item, given room for the item
      else  # no room
        DP[i+1,w+1] = DP[i, w+1]  # take value without the item
      end
    end
  end
  # now, backtrack
  solution = eltype(items_knapsack)[]
  i, w = n, cap
  while i > 0
    if DP[i+1,w+1] > DP[i,w+1]  # item i was added
      item = items_knapsack[i]
      push!(solution, item)
      w -= weight(item)
    end
    i -= 1
    i == 0 || w == 0 && break
  end
  return solution, last(DP)
end

# ╔═╡ a7f61d0a-7706-47d1-8798-264881d032fa
md"Note the restriction on the type of the weigths!"

# ╔═╡ 67ccbd02-d63c-42e3-ae7b-3022d3ec60c0
dynamicprogramming(knapsack)

# ╔═╡ fe8bb7b2-6fd2-4b8f-9ef0-9d2fcc1b744a
md"""
## Exercise: brute-force TSP

Consider the capitals of the provice of Belgium 🇧🇪
"""

# ╔═╡ 66eb2e32-c0bd-42ab-bb67-05dbd01a1d6e
md"**Assignment** Loop over all possible tours to find the shortest one. A tour is just gives the order of the cities indices. For example"

# ╔═╡ 5afd91f6-9e0d-4774-a954-3f00b71dd2e7
md"## Appendix"

# ╔═╡ 8c685b23-0ee2-4027-ad5c-01c1dce6fe3a
stad_coor = Dict(
	"Antwerpen" => (51.219901727372466, 4.413766520738599),
	"Hasselt" => (50.929969726608796, 5.337669056378496),
	"Gent" => (51.046603059146115, 3.7227238267036165),
	"Leuven" => (50.87717306970613, 4.7013098809614045),
	"Brugge"=>(51.20994021593339, 3.2237458378594632),
	"Bergen" => (50.45178175960054, 3.95116527396219),
	"Luik" => (50.63563258400128, 5.58639186132504),
	"Aarlen" => (49.68489470742727, 5.812244206475732),
	"Namen" => (50.45840475199352, 4.858245535731698),
	"Waver" => (50.71402361189344, 4.596358571899753),
	"Brussel" => (50.86323347803335, 4.349845202806156)
);

# ╔═╡ 1a6ded1a-cfbe-4d7b-b70f-b06075747561
cities = keys(stad_coor) |> collect;

# ╔═╡ 8e5dcba4-21dc-4e7c-a3cb-7eec261a270c
cities

# ╔═╡ 15c2faec-655f-4d56-9e4c-5abd4c1d9026
nsteden = ncities = length(cities);

# ╔═╡ 45e7087e-9ad6-49b9-97ac-6a9038d5cb44
a_tour = 1:ncities  # order of the list

# ╔═╡ 827b5a0e-c838-42d2-be03-87a4dd34b24a
dist_sphere((lat1, long1), (lat2, long2); r=6_371) = 2r * asin(√(sin((lat2 -lat1)/2)^2 + cos(lat1)*cos(lat2)*sin((long2-long1)/2)^2));

# ╔═╡ 18e490a3-5823-4862-ae86-07546f556eab
D = [dist_sphere(stad_coor[h1] .* pi ./ 180, stad_coor[h2] .* pi ./ 180) for h1 in cities, h2 in cities];

# ╔═╡ 54db8c49-dee0-4ef1-bd26-78d10b03ae76
D  # distance matrix

# ╔═╡ c59969a6-eb3e-47c5-9bd8-f1af03fce595
function total_distance(route, incomplete=false)
	c = 0.0
	prev = incomplete ? 1 : route[end]
	for next in route
		c += D[prev, next]
		prev = next
	end
	incomplete && (c += D[route[end], 1])
	return c
end

# ╔═╡ b09bae94-4add-4d71-8bbc-5d5e4cdba2f6
total_distance(a_tour)

# ╔═╡ 137dc6ec-72c5-4858-a18a-9fb433899805
function plot_tour(route)
	routec = [route..., route[1]]
	p = plot()

	for h in cities
		lat, long = stad_coor[h]
		scatter!(p, [long], [lat], label="", color="blue")
		annotate!((long, lat, h))
	end

	for (i,j) in zip(routec[1:end-1], routec[2:end])
		h1, h2 = cities[i], cities[j]
		lat1, long1 = stad_coor[h1]
		lat2, long2 = stad_coor[h2]

		plot!(p, [long1, long2], [lat1, lat2], color="red", label="", lw=3, alpha=0.8)
	end
	return p
end;

# ╔═╡ ed475158-07c4-4116-b0eb-d23911bf6531
plot_tour([1 for i in 1:nsteden])

# ╔═╡ cedb1c31-859c-4c6a-a0fd-4c643959c4f7
plot_tour(a_tour)

# ╔═╡ fde5a843-5767-4fc7-a381-eb2ecc1fe801
begin
	myblue = "#304da5"
	mygreen = "#2a9d8f"
	myyellow = "#e9c46a"
	myorange = "#f4a261"
	myred = "#e76f51"
	myblack = "#50514F"

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

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

[compat]
Combinatorics = "~1.0.2"
Plots = "~1.22.3"
"""

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

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

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

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

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

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

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

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

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

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

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

[[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 = "31d0151f5716b655421d9d75b7fa74cc4e744df2"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.39.0"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

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

[[IterTools]]
git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.3.0"

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

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

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

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

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

[[LaTeXStrings]]
git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.2.1"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "9d8c00ef7a8d110787ff6f170579846f776133a9"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.4"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "1162ce4a6c4b7e31e0e6b14486a6986951c73be9"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.5.2"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# ╔═╡ Cell order:
# ╟─936344da-213a-11ec-3e68-05f0f85cf2f3
# ╟─c3b61d51-679c-4aba-957b-105c00037aee
# ╠═0150a3e2-f98e-4412-9e96-1a2db5a9421e
# ╟─62c4926d-0350-4088-8729-0ee5afb306be
# ╠═3c769076-fec9-4730-b993-2fefa09f0c8d
# ╠═ada8c7c4-f6d8-41e7-a4e0-c137ce7ab151
# ╟─64e7be40-cd0b-40c4-b476-51c6a66d40e1
# ╟─be803662-aa08-4ee3-a400-bbc18b7d060f
# ╠═2ce221f2-5712-440e-9b34-168a3a5af65e
# ╟─2d9598de-dfd6-46c7-bf4f-c2f48134e558
# ╠═265ac424-3e7b-4222-ab79-07907544a9fa
# ╠═876bd16d-89ac-4164-8055-e9ec4a4d875f
# ╠═cca6890d-b0eb-4154-9b27-be79c0b8ab38
# ╠═94cb9a4c-aa3f-4ef0-a4e5-170b0d2ad69b
# ╠═1d07e4af-5419-464f-842d-3012920e8a1b
# ╠═f0734b83-dbbc-4d1d-a96e-9209e53b7382
# ╠═165d0c6a-ec73-46a9-a9f4-5085c65d184c
# ╟─e1fa7a85-5aab-4e17-b02f-c9a0c62f3f47
# ╠═38e36108-85cf-4808-8e10-bc6f869522fc
# ╠═9b555941-5ef2-48d5-8efd-0f4cf7e38af6
# ╟─96d84e8c-3b11-4b01-a036-e3434804f512
# ╠═f65c49f0-09bf-4b5f-99b9-21b2fe6f5760
# ╟─3e95cf39-5d8f-498a-bdb6-df99a51257f3
# ╠═97465212-39d3-4b19-8dd0-cf654b80e3e2
# ╟─834e4c69-5e90-409a-a2ed-ac14609ae728
# ╠═02e07bf0-741c-4df2-9a07-55093e2bb55d
# ╟─e53659a3-e791-4318-a738-678d8c267b86
# ╠═fa75a27d-725f-49df-8527-125c90de1692
# ╟─79be9317-762d-4ce5-85ed-983f4882fa0b
# ╠═d500504b-5cfb-46b9-8b82-1fbfb843b479
# ╟─a7f61d0a-7706-47d1-8798-264881d032fa
# ╠═67ccbd02-d63c-42e3-ae7b-3022d3ec60c0
# ╟─fe8bb7b2-6fd2-4b8f-9ef0-9d2fcc1b744a
# ╠═ed475158-07c4-4116-b0eb-d23911bf6531
# ╠═8e5dcba4-21dc-4e7c-a3cb-7eec261a270c
# ╠═54db8c49-dee0-4ef1-bd26-78d10b03ae76
# ╟─66eb2e32-c0bd-42ab-bb67-05dbd01a1d6e
# ╠═45e7087e-9ad6-49b9-97ac-6a9038d5cb44
# ╠═cedb1c31-859c-4c6a-a0fd-4c643959c4f7
# ╠═b09bae94-4add-4d71-8bbc-5d5e4cdba2f6
# ╟─5afd91f6-9e0d-4774-a954-3f00b71dd2e7
# ╟─8c685b23-0ee2-4027-ad5c-01c1dce6fe3a
# ╟─1a6ded1a-cfbe-4d7b-b70f-b06075747561
# ╟─15c2faec-655f-4d56-9e4c-5abd4c1d9026
# ╟─827b5a0e-c838-42d2-be03-87a4dd34b24a
# ╟─c59969a6-eb3e-47c5-9bd8-f1af03fce595
# ╟─18e490a3-5823-4862-ae86-07546f556eab
# ╟─137dc6ec-72c5-4858-a18a-9fb433899805
# ╟─fde5a843-5767-4fc7-a381-eb2ecc1fe801
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: MichielStock/STMO

max_issues_repo_name: MichielStock/STMO

max_forks_repo_name: MichielStock/STMO


text:
abstract type AngularMomentum{label} <: Tensor{1,label} end

# * Orbital Angular Momentum

"""
    OrbitalAngularMomentum()

The angular momentum ``\\tensor{L}`` is associated with the spatial
coordinates ``\\theta`` and ``\\phi``.
"""
struct OrbitalAngularMomentum <: AngularMomentum{'L'} end
system(::Type{OrbitalAngularMomentum}) = OrbitalAngularMomentumSubSystem()

@tensor(OrbitalAngularMomentum) do
    ℓ′ == ℓ

    raw"""
    rme(ℓ′, 𝐋̂, ℓ)

Calculate the reduced matrix element of the orbital angular momentum:

```math
\begin{equation}
\tag{V13.2.67}
\redmatrixel{\ell'}{\tensor{L}}{\ell} \defd
\delta_{\ell'\ell}\sqrt{\ell(\ell+1)(2\ell+1)}
\end{equation}
```

# Examples

```jldoctest
julia> rme(1, OrbitalAngularMomentum(), 1)
2.449489742783178

julia> rme(1, OrbitalAngularMomentum(), 2)
0
```
"""
    √(ℓ*(ℓ+1)*(2ℓ+1))
end

# * Spin Angular Momentum

"""
    SpinAngularMomentum()

The spin angular momentum ``\\tensor{S}`` is the intrinsic angular
momentum associated with the coordinate ``s``.
"""
struct SpinAngularMomentum <: AngularMomentum{'S'} end
system(::Type{SpinAngularMomentum}) = SpinSubSystem()

@tensor(SpinAngularMomentum) do
    s′ == s

    raw"""
    rme(s′, ::SpinAngularMomentum, s)

Calculate the reduced matrix element of the spin angular momentum:

```math
\begin{equation}
\tag{V13.2.95}
\redmatrixel{s'}{\tensor{S}}{s} \defd
\delta_{ss'}
\sqrt{s(s+1)(2s+1)}
\end{equation}
```

# Examples

```jldoctest
julia> rme(half(1), SpinAngularMomentum(), half(1))
1.224744871391589

julia> rme(half(1), SpinAngularMomentum(), half(3))
0
```
"""
    √(s*(s+1)*(2s+1))
end

# * Total Angular Momentum

"""
    TotalAngularMomentum()

The total angular momentum ``\\tensor{J} = \\tensor{L} + \\tensor{S}``
results from the coupling of the orbital and spin angular momenta.
"""
struct TotalAngularMomentum <: AngularMomentum{'J'} end
system(::Type{TotalAngularMomentum}) = TotalAngularMomentumSubSystem()

@tensor(TotalAngularMomentum) do
    begin
        ℓ′ == ℓ
        s′ == s
        J′ == J
    end

    raw"""
    rme((ℓ′,s′,J′), ::TotalAngularMomentum, (ℓ,s,J))

Calculate the reduced matrix element of the total angular momentum:

```math
\begin{equation}
\tag{V13.2.38}
\redmatrixel{\ell's'J'}{\tensor{J}}{\ell s J} \defd
\delta_{\ell\ell'}\delta_{ss'}\delta_{JJ'}
\sqrt{J(J+1)(2J+1)}
\end{equation}
```

# Examples

```jldoctest
julia> rme((1,half(1),half(3)), TotalAngularMomentum(), (1,half(1),half(3)))
3.872983346207417

julia> rme((1,half(1),half(3)), TotalAngularMomentum(), (1,half(1),half(1)))
0
```

"""
    √(J*(J+1)*(2J+1))
end

export OrbitalAngularMomentum, SpinAngularMomentum, TotalAngularMomentum

ext: jl


text:
#!/usr/bin/env julia --startup-file=no

"""
Linear-Stability-Calculators: Quasi-Hydrostatic Model
=====================================================

Linear Stability Calculation of an inertially unstable jet using spectral methods

The basic state is that of a barotropic Bickley jet (to be generalized).
The perturbation is assumed to be periodic in both the zonal (x) and vertical (z) directions.

Calculations are done in parallel.  See inertial_instability.py for serial execution.
"""

# put this somewhere else?
"""
    split_count(N::Integer, n::Integer)
Return a vector of `n` integers which are approximately equally sized and sum to `N`.
"""
function split_count(N::Integer, n::Integer)
    q,r = divrem(N, n)
    return [i <= r ? q+1 : q for i = 1:n]
end

using MPI
using Printf

MPI.Init()

comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
comm_size = MPI.Comm_size(comm)

include("src/quasihydrostatic_setup.jl")
include("src/parameters.jl")
include("src/mpi_stuff.jl")
include("src/Compute_Eigenvalues.jl")
include("src/Plotting_Functions.jl")

files   = Files()
grid    = Grid(Ly = 1000e3, Lz = 3e3, Ny = 200, θ₀ = π/32, method=Cheb())
physics = Physics(N = 1e-4, ν = 0.26, θ₀ = grid.θ₀, NT = 1)
jet     = Jet(grid, physics)

rank == 0 ? Output_Parameters(grid, physics, jet, files.json) : nothing

Neigs  = 10
dk, Nk = 1e-6, 41
dm, Nm = 1e-4, 150

ms = collect(dm:dm:Nm*dm)

if rank == 0
    ks = Float64[(ik-1)*dk for ik in 1:Nk]
    Nk_counts = split_count(Nk, comm_size)
    ks_vbuf   = VBuffer(ks, Nk_counts)        # VBuffer for scatter

else
    Nk_counts = nothing
    ks_vbuf   = nothing
end

Nk_local, = MPI.Scatter!(Nk_counts, zeros(Int, 1), 0, comm)
ks_local = MPI.Scatterv!(ks_vbuf, zeros(Float64, Nk_local), 0, comm)

#@printf("rank = %s, Nk_local = %s, ks_local = %s\n", rank, Nk_local, ks_local)

ωs_local    = zeros(Complex{Float64}, (Neigs, Nm, length(ks_local)))
modes_local = zeros(Complex{Float64}, (3*grid.Ny+1, Neigs, Nm, length(ks_local)))

for (ik, k) in enumerate(ks_local)
    for (im, m) in enumerate(ms)
        A = build_matrix(k, m, grid, physics, jet)

        ωs_local[:, im, ik], modes_local[:, :, im, ik] = compute_spectrum(A, Neigs)

        #print("k = ", ks_local[ik], " m = ", ms[im], " growth = ", imag(ωs_local[1, im, ik]), "\n")
    end
end

MPI.Barrier(comm)

if rank == 0
    ωs    = zeros(Complex{Float64}, (Neigs, Nm, Nk))
    modes = zeros(Complex{Float64}, (3*grid.Ny+1, Neigs, Nm, Nk))

    Neigs_counts = [Neigs for k = 1:comm_size]
    Nm_counts    = [Nm    for j = 1:comm_size]

    ωs_sizes  = vcat(Neigs_counts', Nm_counts', Nk_counts')
    ωs_counts = vec(prod(ωs_sizes, dims=1))
    ωs_vbuf   = VBuffer(ωs, ωs_counts)

    Ny_counts = [3*grid.Ny+1 for k = 1:comm_size]

    modes_sizes  = vcat(Ny_counts', Neigs_counts', Nm_counts', Nk_counts')
    modes_counts = vec(prod(modes_sizes, dims=1))
    modes_vbuf   = VBuffer(modes, modes_counts)

else
    ωs_sizes = nothing
    ωs_vbuf  = nothing

    modes_sizes = nothing
    modes_vbuf  = nothing
end

MPI.Gatherv!(ωs_local,    ωs_vbuf,    0, comm)
MPI.Gatherv!(modes_local, modes_vbuf, 0, comm)

MPI.Barrier(comm)

if rank == 0
    save_spectrum(ωs, modes, ks, ms, Neigs, grid.y, grid.Ny, files.nc)

    plot_growth_slice(files.nc, files.json, files.plotslicem)
    plot_growth(files.nc, files.json, files.plotgrowth)

    plot_modes_1D(files.nc, files.json, files.plotmodes1D, Neigs)
    plot_modes_2D(files.nc, files.json, files.plotmodes2D, Neigs)
end

MPI.Finalize()


text:
function plot_coupler(c::Coupler, f₀ = 1.0)
    # Constants
    PTS = 1000

    # Analysis functions
    funs = analysis_functions(c)

    # Create x-axis
    f̄ = range(0.7, 1.3, length = PTS)
    f = f̄ .* f₀
    t = im .* tan.(π/4 .* f̄)

    # Labels
    titles = ["Return losses", "Directivity", "Insertion Losses", "Coupling"]
    labels = ["Γ", "D", "IL", "C"]

    # CREATE FIGURE
    fig = Figure()

    # CREATE TITLE
    supertitle = Label(fig[1,:], _title, textsize = 30)

    # CREATE AXES
    plot_grid = GridLayout()
    fig.layout[2,:] = plot_grid
    axs = [Axis(plot_grid[reverse(fldmod1(i,2))...], title = titles[i]) for i in 1:4]

    # CREATE LABELS
    ab_display = GridLayout()
    fig.layout[3,:] = ab_display
    a_labels = [ Label(ab_display[1,i],
                    @sprintf("%8.4f", a[i]),
                    textsize = 10,
                    tellwidth = false) for i in 1:10]
    b_labels = [ Label(ab_display[2,i],
                    @sprintf("%8.4f", b[i]),
                    textsize = 10,
                    tellwidth = false) for i in 1:10]

    # PLOT EVERYTHING
    _update_plot!(axs, a_labels, b_labels, c, labels, t, f)

    # RETURN WHOLE FIGURE
    return fig
end

function _update_plot!(axs, a_labels, b_labels, coupler, labels, t, f)
    funs = analysis_functions(coupler)

    for (ax_, fn_, lb_) in zip(axs, funs, labels)
        empty!(ax_)
        lines!(ax_, f, db.(fn_.(t)), label = lb_)
        axislegend(ax_)
        xlims!(ax_, f[1], f[end])
    end

    ylims!(axs[1], -50, -10)
    ylims!(axs[2], -50, -10)

    for i in eachindex(a_labels)
        checkindex(Bool, eachindex(coupler.a), i) ?
        a_labels[i].text = @sprintf("%8.4f", coupler.a[i]) :
        a_labels[i].text = ""
    end
    for i in eachindex(b_labels)
        checkindex(Bool, eachindex(coupler.b), i) ?
            (b_labels[i].text = @sprintf("%8.4f", coupler.b[i]);
                b_labels[i].color = (coupler.b[i] > 1) ? (:red) : (:blue)) :
            b_labels[i].text = ""

    end
end


text:
module ReconTPCF

# Write your package code here.

export loadim, SN_comp, S2_initialise, S2_finalise
export C2_initialise, naninf, blas_stat_st2
export pre_proc
export get_C2_S2, histrecon

include("./CorrFuncs.jl")
include("./histrecon.jl")

end


text:
@testset "inference" begin

    @testset "BinaryHDV" begin
        x = BinaryHDV([true, false, true, true])
        y = BinaryHDV([false, false, false, true])

        @test similarity(x, y) ≈ 1 / 3 ≈ jacc_sim(x.v, y.v)

        Π!(x, 2)
        @test similarity(x, y) ≈ jacc_sim(circshift(x.v, 2), y.v)

    end

    @testset "GradedHDV" begin
        x = GradedHDV([0.1, 0.4, 0.6, .8])
        y = GradedHDV([0.9, 0.8, 0.1, 0.3])

        @test similarity(x, y) ≈ jacc_sim(x.v, y.v)

        Π!(x, 2)
        @test similarity(x, y) ≈ jacc_sim(circshift(x.v, 2), y.v)
    end

    @testset "BipolarHDV" begin
        x = BipolarHDV([1, -1, 0, -1])
        y = BipolarHDV([-1, 0, 1, -1])

        @test similarity(x, y) ≈ cos_sim(x.v, y.v)

        Π!(x, 2)
        @test similarity(x, y) ≈ cos_sim(circshift(x.v, 2), y.v)
    end

    @testset "GradedBipolarHDV" begin
        x = GradedBipolarHDV([0.1, -0.4, 0.6, .8])
        y = GradedBipolarHDV([0.9, 0.8, -0.1, -0.3])

        @test similarity(x, y) ≈ cos_sim(x.v, y.v)

        Π!(x, 2)
        @test similarity(x, y) ≈ cos_sim(circshift(x.v, 2), y.v)
    end

    @testset "RealHDV" begin
        x = RealHDV([0.1, -0.4, 0.6, .8])
        y = RealHDV([0.9, 0.8, -0.1, -0.3])

        @test similarity(x, y) ≈ cos_sim(x.v, y.v)

        Π!(x, 2)
        @test similarity(x, y) ≈ cos_sim(circshift(x.v, 2), y.v)
    end
end


text:
using Bokeh

Bokeh.includejs(true)  # this is required as these files will be published online
Bokeh.noshow(true)
Bokeh.filewarnings(false)
exdir = dirname(@__FILE__)
cd(exdir)
for jl_src in filter(f -> endswith(f, ".jl") && f != "build_all.jl", readdir(exdir))
	glyphsize(6)
	hold(false)
	plotfile(replace(jl_src, ".jl", ".html"))
	evalfile(jl_src)
end


text:
#using TikzGraphs
#using TikzPictures

"""
    plot_dag(g, node_labels)

Plot DAG or the output of the PC algorithm.
"""
function plot_dag(g, node_labels::Array=[])
    plot_g = DiGraph(nv(g))

    if length(node_labels) != nv(g)
        node_labels = map(string, 1:nv(g))
    end

    node_style = "draw, rounded corners, fill=blue!10"
    options = "scale=2"

    styles_dict = Dict()

    for e in edges(g)
       if has_edge(plot_g, e.dst, e.src)
            styles_dict[(e.dst, e.src)] = "o-o"
       else
           add_edge!(plot_g, e.src, e.dst)
           push!(styles_dict, (e.src, e.dst)=>"->")
       end
    end

    #TikzGraphs.plot
    (plot_g, node_labels, edge_styles=styles_dict,
                    node_style=node_style, options=options)
end


text:
#!/usr/bin/env julia

function energy_makerecipe(result, recipes)
    for i in recipes[result]

    end
end

function solve(case, junctions, stones, recipes)
    energy = energy_makerecipe(1, recipes)
    println("Case #$case: $(energy > 10^12 ? -1 : energy)")
end

function main()
    T = parse(Int, readline(stdin))
    for case in 1:T
        N, M, S, R = parse.(Int, stdin |> readline |> split)

        # U -> (V, junction number) for each adjacent junction and vice versa
        junctions = Dict{Int,Vector{Int}}()
        for i in 1:M
            U, V = parse.(Int, stdin |> readline |> split)
            push!(get!(junctions, U, []), (V, i))
            push!(get!(junctions, V, []), (U, i))
        end

        # junction number -> stones
        stones = Dict{Int,Vector{Int}}()
        for i in 1:N
            C = parse.(Int, stdin |> readline |> split)
            map(x -> push!(get!(stones, C[x], []), i), 2:(C[1] + 1))
        end

        # resultant stone -> needed stones
        recipes = Dict{Int, Vector{Int}}()
        for _ in 1:R
            K = parse.(Int, stdin |> readline |> split)
            result = pop!(K)
            popfirst!(K)
            recipes[result] = K
        end

        solve(case, junctions, stones, recipes)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
using Meshing
using Test
using GeometryTypes
using ForwardDiff
using Profile
using Statistics: mean
using LinearAlgebra: dot, norm

@testset "meshing" begin
    @testset "surface nets" begin
        sdf_sphere = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v
            sqrt(sum(dot(v,v))) - 1 # sphere
        end
        sdf_torus = SignedDistanceField(HyperRectangle(Vec(-2,-2,-2.),Vec(4,4,4.)), 0.05) do v
            (sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25
        end

        snf_torus = HomogenousMesh(HyperRectangle(Vec(-2,-2,-2.),Vec(4,4,4.)), (81,81,81), NaiveSurfaceNets()) do v
            (sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25
        end

        snf_sphere = HomogenousMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), (21,21,21), NaiveSurfaceNets()) do v
            sqrt(sum(dot(v,v))) - 1 # sphere
        end

        # test convience constructors
        HomogenousMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), NaiveSurfaceNets()) do v
            sqrt(sum(dot(v,v))) - 1 # sphere
        end
        HomogenousMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), NaiveSurfaceNets(), size=(5,5,5)) do v
            sqrt(sum(dot(v,v))) - 1 # sphere
        end

        sphere = HomogenousMesh(sdf_sphere, NaiveSurfaceNets())
        torus = HomogenousMesh(sdf_torus, NaiveSurfaceNets())
        @test length(vertices(sphere)) == 1832
        @test length(vertices(torus)) == 5532
        @test length(faces(sphere)) == 1830
        @test length(faces(torus)) == 5532

        @test length(vertices(sphere)) == length(vertices(snf_sphere))
        @test length(vertices(torus)) == length(vertices(snf_torus))
        @test length(faces(sphere)) == length(faces(snf_sphere))
        @test length(faces(torus)) == length(faces(snf_torus))
    end

    @testset "noisy spheres" begin
        # Produce a level set function that is a noisy version of the distance from
        # the origin (such that level sets are noisy spheres).
        #
        # The noise should exercise marching tetrahedra's ability to produce a water-
        # tight surface in all cases (unlike standard marching cubes).
        #
        N = 10
        sigma = 1.0
        distance = Float32[ sqrt(Float32(i*i+j*j+k*k)) for i = -N:N, j = -N:N, k = -N:N ]
        distance = distance + sigma*rand(2*N+1,2*N+1,2*N+1)

        # Extract an isosurface.
        #
        lambda = N-2*sigma # isovalue

        msh = HomogenousMesh(distance, MarchingTetrahedra(lambda))

        s2 = SignedDistanceField(HyperRectangle(Vec(0,0,0.),Vec(1,1,1.))) do v
            sqrt(sum(dot(v,v))) - 1 # sphere
        end

        if "--profile" in ARGS
            HomogenousMesh(s2, MarchingTetrahedra())
            Profile.clear()
            @profile HomogenousMesh(s2, MarchingTetrahedra())
            #ProfileView.view()
        end

        msh = HomogenousMesh(s2, MarchingTetrahedra())
        @test length(vertices(msh)) == 973
        @test length(faces(msh)) == 1830
    end

    @testset "vertex interpolation" begin
        @test Meshing.vertex_interp(0, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,0.5,0)
        @test Meshing.vertex_interp(-1, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,0,0)
        @test Meshing.vertex_interp(1, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,1,0)
    end

    @testset "marching cubes" begin

        algo = MarchingCubes()

        sdf = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v
            sqrt(sum(dot(v,v))) - 1 # sphere
        end

        mf = SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)),algo, size=(21,21,21)) do v
            sqrt(sum(dot(v,v))) - 1 # sphere
        end

        mfrv = SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)),MarchingCubes(reduceverts=false), size=(21,21,21)) do v
            sqrt(sum(dot(v,v))) - 1 # sphere
        end

        # convience constructors
        SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), MarchingCubes()) do v
            sqrt(sum(dot(v,v))) - 1 # sphere
        end
        SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), MarchingCubes(), size=(5,6,7)) do v
            sqrt(sum(dot(v,v))) - 1 # sphere
        end

        @test length(vertices(mfrv)) == 10968
        m = SimpleMesh(sdf,algo)
        m2 = SimpleMesh(sdf,algo)
        @test length(vertices(m)) == 7320
        @test length(faces(m)) == 3656
        @test length(faces(mf)) == length(faces(mfrv))
        @test m == m2
        @test length(vertices(m)) == length(vertices(mf))
        @test length(faces(m)) == length(faces(mf))
    end

    @testset "respect origin" begin
        # verify that when we construct a mesh, that mesh:
        #   a) respects the origin of the SDF
        #   b) is correctly spaced so that symmetric functions create symmetric meshes
        f = x -> norm(x)
        bounds = HyperRectangle(Vec(-1, -1, -1), Vec(2, 2, 2))
        resolution = 0.1
        sdf = SignedDistanceField(f, bounds, resolution)

        for algorithm in (MarchingCubes(0.5), MarchingTetrahedra(0.5))
            mesh = @inferred GLNormalMesh(sdf, algorithm)
            # should be centered on the origin
            @test mean(vertices(mesh)) ≈ [0, 0, 0] atol=0.15*resolution
            # and should be symmetric about the origin
            @test maximum(vertices(mesh)) ≈ [0.5, 0.5, 0.5]
            @test minimum(vertices(mesh)) ≈ [-0.5, -0.5, -0.5]
        end

        # Test functional Meshing
        mesh = @inferred GLNormalMesh(f, HyperRectangle(Vec(-1, -1, -1), Vec(2, 2, 2)), (21,21,21), MarchingCubes(0.5))
        # should be centered on the origin
        @test mean(vertices(mesh)) ≈ [0, 0, 0] atol=0.15*resolution
        # and should be symmetric about the origin
        @test maximum(vertices(mesh)) ≈ [0.5, 0.5, 0.5]
        @test minimum(vertices(mesh)) ≈ [-0.5, -0.5, -0.5]

        # Naive Surface Nets has no accuracy guarantee, and is a weighted sum
        # so a larger tolerance is needed for this one. In addition,
        # quad -> triangle conversion is not functioning correctly
        # see: https://github.com/JuliaGeometry/GeometryTypes.jl/issues/169
        mesh = @inferred GLNormalMesh(sdf, NaiveSurfaceNets(0.5))
        # should be centered on the origin
        @test mean(vertices(mesh)) ≈ [0, 0, 0] atol=0.15*resolution
        # and should be symmetric about the origin
        @test maximum(vertices(mesh)) ≈ [0.5, 0.5, 0.5] atol=0.2
        @test minimum(vertices(mesh)) ≈ [-0.5, -0.5, -0.5] atol=0.2
    end

    @testset "AbstractMeshingAlgorithm interface" begin
        f = x -> norm(x) - 0.5
        bounds = HyperRectangle(Vec(-1, -1, -1), Vec(2, 2, 2))
        resolution = 0.1
        sdf = SignedDistanceField(f, bounds, resolution)

        @testset "marching cubes" begin
            @test_nowarn GLNormalMesh(sdf, MarchingCubes())
            @inferred GLNormalMesh(sdf, MarchingCubes())
        end

        @testset "marching tetrahedra" begin
            @test_nowarn GLNormalMesh(sdf, MarchingTetrahedra())
            @inferred GLNormalMesh(sdf, MarchingTetrahedra())
            @test_nowarn GLNormalMesh(sdf.data, MarchingTetrahedra(0.5))
            @inferred GLNormalMesh(sdf.data, MarchingTetrahedra(0.5))
        end
        @testset "naive surface nets" begin
            @test_nowarn GLNormalMesh(sdf, NaiveSurfaceNets())
            @inferred GLNormalMesh(sdf, NaiveSurfaceNets())
        end
    end

    @testset "mixed types" begin
        # https://github.com/JuliaGeometry/Meshing.jl/issues/9
        samples = randn(10, 10, 10)
        m1 = HomogenousMesh(SignedDistanceField(
            HyperRectangle(Vec(0,0,0), Vec(1, 1, 1)),
            samples), MarchingTetrahedra())
        m2 = HomogenousMesh(SignedDistanceField(
            HyperRectangle(Vec(0,0,0), Vec(1, 1, 1)),
            Float32.(samples)), MarchingTetrahedra())
        @test all(isapprox.(vertices(m1), vertices(m2), rtol=1e-6))
        @test all(faces(m1) == faces(m2))

        @testset "forward diff" begin
            # Demonstrate that we can even take derivatives through the meshing algorithm
            f = x -> norm(x)
            bounds = HyperRectangle(Vec(-1, -1, -1), Vec(2, 2, 2))
            resolution = 0.1
            sdf = SignedDistanceField(f, bounds, resolution)

            function surface_distance_from_origin(isovalue)
                # Compute the mean distance of each vertex in the isosurface
                # mesh from the origin. This should return a value equal to the
                # isovalue and should have a derivative of 1.0 w.r.t. the
                # isovalue. This function is just meant to serve as an example
                # of some quantity you might want to differentiate in a mesh,
                # and has the benefit for testing of having a trivial expected
                # solution.
                mesh = HomogenousMesh(sdf, MarchingTetrahedra(isovalue))
                mean(norm.(vertices(mesh)))
            end

            isoval = 0.85
            @test surface_distance_from_origin(isoval) ≈ isoval atol=1e-2
            @test ForwardDiff.derivative(surface_distance_from_origin, isoval) ≈ 1 atol=1e-2
        end

        @testset "type stability" begin
            # verify that we don't lose type stability just by mixing arguments
            # of different types
            data = randn(5, 5, 5)
            iso = 0.2
            eps = 1e-4
            @inferred(Meshing.marchingTetrahedra(data, iso, eps, Int))
            @inferred(Meshing.marchingTetrahedra(Float32.(data), Float64(iso), Float16(eps), Int32))
            @inferred(Meshing.marchingTetrahedra(Float64.(data), Float32(iso), Float64(eps), Int64))
        end
        @testset "Float16" begin
            sdf_torus = SignedDistanceField(HyperRectangle(Vec{3,Float16}(-2,-2,-2.),
                                                           Vec{3,Float16}(4,4,4.)),
                                            0.1, Float16) do v
                (sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25
            end
            @test typeof(HomogenousMesh(sdf_torus,NaiveSurfaceNets())) ==
                         PlainMesh{Float16,Face{4,Int}}
            m2 = HomogenousMesh(sdf_torus,MarchingTetrahedra())
            m3 = HomogenousMesh(sdf_torus,MarchingCubes())
        end
    end
end


text:
#=
 Train a simple network to predict formation energy per atom (downloaded from Materials Project).
=#
using CSV, DataFrames
using Random, Statistics
using Flux
using Flux: @epochs
using SimpleWeightedGraphs
using ChemistryFeaturization
using AtomicGraphNets

println("Setting things up...")

# data-related options
num_pts = 100 # how many points to use?
train_frac = 0.8 # what fraction for training?
num_epochs = 5 # how many epochs to train?
num_train = Int32(round(train_frac * num_pts))
num_test = num_pts - num_train
prop = "formation_energy_per_atom"
datadir = "$(@__DIR__)/data/"
id = "task_id" # field by which to label each input material

# atom featurization, pretty arbitrary choices for now
features = Symbol.(["Group", "Row", "Block", "Atomic mass", "Atomic radius", "X"])
num_bins = [18, 9, 4, 16, 10, 10]
num_features = sum(num_bins) # we'll use this later
logspaced = [false, false, false, true, true, false]
# returns actual vectors (in a dict with keys of elements) plus Vector of AtomFeat objects describing featurization metadata
atom_feature_vecs, featurization = make_feature_vectors(features, nbins=num_bins, logspaced=logspaced)

# model hyperparameters – keeping it pretty simple for now
num_conv = 3 # how many convolutional layers?
crys_fea_len = 32 # length of crystal feature vector after pooling (keep node dimension constant for now)
num_hidden_layers = 1 # how many fully-connected layers after convolution and pooling?
opt = ADAM(0.001) # optimizer

# dataset...first, read in outputs
info = CSV.read(string(datadir,prop,".csv"), DataFrame)
y = Array(Float32.(info[!, Symbol(prop)]))

# shuffle data and pick out subset
indices = shuffle(1:size(info,1))[1:num_pts]
info = info[indices,:]
output = y[indices]

# next, make graphs and build input features (matrices of dimension (# features, # nodes))
println("Building graphs and feature vectors from structures...")
inputs = AtomGraph[]

for r in eachrow(info)
    cifpath = string(datadir, prop, "_cifs/", r[Symbol(id)], ".cif")
    gr = build_graph(cifpath)
    feature_mat = hcat([atom_feature_vecs[e] for e in gr.elements]...)
    add_features!(gr, feature_mat, featurization)
    push!(inputs, gr)
end

# pick out train/test sets
println("Dividing into train/test sets...")
train_output = output[1:num_train]
test_output = output[num_train+1:end]
train_input = inputs[1:num_train]
test_input = inputs[num_train+1:end]
train_data = zip(train_input, train_output)

# build the model
println("Building the network...")
model = Xie_model(num_features, num_conv=num_conv, atom_conv_feature_length=crys_fea_len, pooled_feature_length=(Int(crys_fea_len/2)), num_hidden_layers=1)

# define loss function and a callback to monitor progress
loss(x,y) = Flux.Losses.mse(model(x), y)
evalcb() = @show(mean(loss.(test_input, test_output)))
evalcb()

# train
println("Training!")
@epochs num_epochs Flux.train!(loss, params(model), train_data, opt, cb = Flux.throttle(evalcb, 5))


text:
using CxxrtlCosim
using ModelingToolkit
using OrdinaryDiffEq

cd(@__DIR__)

if !isfile("blinky.so")
    println("compiling blinky.v")
    CxxrtlCosim.synth("blinky.cpp", "blinky.v")
    CxxrtlCosim.compile("blinky.so", "blinky.cpp")
end

top = CxxrtlCosim.cxxrtl_design_create("./blinky.so")

blink = CxxrtlCosim.cxxrtl_create(top)
CxxrtlCosim.cxxrtl_reset(blink)

@parameters t::Int
@variables clk(t)::Int clkout(t)::Int led(t)::Int

D = Difference(t; dt=1)

eqs = [
    D(clk) ~ mod(clk+1, 2),
    D(clkout) ~ CxxrtlCosim.rtlbind(t, clk, blink, "clk"),
    D(led) ~ CxxrtlCosim.rtlget(t, blink, "led"),
]

@named de = DiscreteSystem(eqs,t,[clk, clkout, led], [])

u0 = [0, 0, 0]
prob = DiscreteProblem(de, u0, (0, 10000))
sol = solve(prob, FunctionMap())

max_stars_repo_name: JuliaComputing/CxxrtlCosim

max_issues_repo_name: JuliaComputing/CxxrtlCosim

max_forks_repo_name: JuliaComputing/CxxrtlCosim


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

using Markdown
using InteractiveUtils

# ╔═╡ f36196a2-bd8a-11eb-35d3-77b0357b31b6
begin
    import Pkg
    Pkg.activate(mktempdir())

    Pkg.add([
        Pkg.PackageSpec(name="CairoMakie", version="0.5"),
        Pkg.PackageSpec(name="DataFrames", version="1"),
        Pkg.PackageSpec(name="PalmerPenguins", version="0.1"),
        Pkg.PackageSpec(name="AlgebraOfGraphics", version="0.4"),
    ])

    using CairoMakie, DataFrames, PalmerPenguins, AlgebraOfGraphics
	set_aog_theme!()
end

# ╔═╡ a75244d9-3db4-4f0e-ac6d-e21ce3d2e742
begin
	penguins = dropmissing(DataFrame(PalmerPenguins.load()))
	first(penguins, 6)
end

# ╔═╡ f901fd37-2ee0-49fc-8714-a487c35cff0e
begin
	axis = (width = 225, height = 225)
	penguin_frequency = data(penguins) * frequency() * mapping(:species)

	draw(penguin_frequency; axis)
end

# ╔═╡ 6efedce7-d3c6-42cc-a5a3-a8ac2c78ea8c
begin
	plt = penguin_frequency * mapping(color = :island, stack = :island)
	draw(plt; axis)
end

# ╔═╡ 134afff2-6ee9-40ff-8d88-213d67f5938d

# ╔═╡ Cell order:
# ╠═f36196a2-bd8a-11eb-35d3-77b0357b31b6
# ╠═a75244d9-3db4-4f0e-ac6d-e21ce3d2e742
# ╠═f901fd37-2ee0-49fc-8714-a487c35cff0e
# ╠═6efedce7-d3c6-42cc-a5a3-a8ac2c78ea8c
# ╠═134afff2-6ee9-40ff-8d88-213d67f5938d

max_stars_repo_name: StatisticalMice/franklin-pluto-demo

max_issues_repo_name: StatisticalMice/franklin-pluto-demo

max_forks_repo_name: StatisticalMice/franklin-pluto-demo


text:
function fib()
    a, b = 0, 1
    while true
        produce(b)
        a, b = b, a+b
    end
end

f = Task(fib)
consume(f) # 1
consume(f) # 1
consume(f) # 2
consume(f) # 3
consume(f) # 5
consume(f) # 8


text:
"""
    AbstractPartition{N,A}

Container to store a set of objects across different partitions. Needs a field
`x`.
"""
abstract type AbstractPartition{N,A} end

@inline Base.iterate(ff::AbstractPartition)         = iterate(ff.x)
@inline Base.iterate(ff::AbstractPartition, i::Int) = iterate(ff.x, i)

@inline Base.length(ff::AbstractPartition{N}) where{N} = N
@inline Base.firstindex(ff::AbstractPartition) = 1
@inline Base.lastindex(ff::AbstractPartition)  = length(ff)

@inline Base.getindex(ff::AbstractPartition, i::Int) = ff.x[i]

"""
    FlattenedVector{T,N,A} <: AbstractVector{T}

Container to store a set of `Arrays` as elements of an `NTuple`, or similar. The
idea is to treat them as a single column vector for all intents and purposes.
Inspired in `ArrayPartition` from `RecursiveArrayTools`. Needs a field `x`.
"""
abstract type FlattenedVector{T,N,A} <: AbstractVector{T} end

@inline Base.length(ff::FlattenedVector) = sum((length(x) for x in ff.x))
@inline Base.size(ff::FlattenedVector)   = (length(ff),)

# indexing. this is just a linear indexing through all the arrays. adapted from
# RecursiveArrayTools
@inline Base.firstindex(ff::FlattenedVector) = 1
@inline Base.lastindex(ff::FlattenedVector)  = length(ff)

@inline function Base.getindex(ff::FlattenedVector, i::Int)
    @inbounds for j in 1:length(ff.x)
        f  = ff.x[j]
        i -= length(f)
        if i <= 0
            return f[length(f)+i]
        end
    end
end
@inline function Base.setindex!(ff::FlattenedVector, v, i::Int)
    @inbounds for j in 1:length(ff.x)
        f  = ff.x[j]
        i -= length(f)
        if i <= 0
            f[length(f)+i] = v
            break
        end
    end
end


text:
struct Thing1
  other_thing::Thing2
end


text:
using Plots
using Statistics: mean
import MyterialColors: blue_grey_darker, blue_dark, green
import InfiniteOpt: termination_status
using CSV
import Colors: HSL
import DataFrames: DataFrame
using Term.progress
import Term.progress: with as withpbar
import Term: install_term_logger, install_stacktrace
import Printf: @sprintf
install_term_logger()
install_stacktrace()

using jcontrol
import jcontrol.comparisons: track_segments, TrackSegment
using jcontrol.visuals
import jcontrol: trimtrial, euclidean, FULLTRACK, Track, trim, toDict, ControlOptions
import jcontrol.Run: run_mtm
import jcontrol.forwardmodel: solution2state, solution2s
import jcontrol.io: PATHS

@Base.kwdef mutable struct SimTracker
    s0::Float64      = 0
    s::Float64      = 0
    iter::Int       = 0
    t::Float64      = 0.0
    Δt::Float64     = 0.01
    prevsol         = nothing
    hasfailed       = 0
end

@Base.kwdef mutable struct SolutionTracker
    t::Vector{Float64} = Vector{Float64}[]
    s::Vector{Float64} = Vector{Float64}[]
    x::Vector{Float64} = Vector{Float64}[]
    y::Vector{Float64} = Vector{Float64}[]
    θ::Vector{Float64} = Vector{Float64}[]
    δ::Vector{Float64} = Vector{Float64}[]
    u::Vector{Float64} = Vector{Float64}[]
    ω::Vector{Float64} = Vector{Float64}[]
    n::Vector{Float64} = Vector{Float64}[]
    ψ::Vector{Float64} = Vector{Float64}[]
    β::Vector{Float64} = Vector{Float64}[]
    v::Vector{Float64} = Vector{Float64}[]
    Fu::Vector{Float64} = Vector{Float64}[]
    δ̇::Vector{Float64} = Vector{Float64}[]
end

function add!(tracker::SolutionTracker, t::Float64, state::State, s::Float64)
    push!(tracker.t, t)

    for v in (:y, :x, :θ, :δ, :u, :ω, :n, :ψ, :β, :v, :Fu, :δ̇)
        push!(getfield(tracker, v), getfield(state, v)[1])
    end

    push!(tracker.s, s)
end

"""
take states at Δt interval while 0 < t < tf and store in tracker
"""
function extend_with_sol!(tracker::SolutionTracker, sol::Solution, Δt::Float64, tf::Float64)
    @info "Before extend with sol" length(tracker.s) tracker.s[end-4:end]
    Δt > tf && return

    for t in Δt:Δt:tf
        state = solution2state(t, sol; at=:time)
        s = solution2s(t, sol)
        add!(tracker, t, state, s)
    end
    @info "After extend with sol" length(tracker.s) tracker.s[end-4:end]
end

"""
Extends the tracker's solution with a whole other solution
"""
function extend_with_sol!(tracker::SolutionTracker, sol::Solution, Δt::Float64)
    N = (Int64 ∘ floor)(Δt / sol.δt)

    for i in 1:N:length(sol.s)
        t = sol.t[i]
        state = solution2state(t, sol; at=:time)
        s = solution2s(t, sol)
        add!(tracker, t, state, s)
    end
end

"""
Repeate MTM with windows of decreasing length
"""
function attempt_step(simtracker, control_model, s0, initial_state, planning_horizon)
    @warn "Could not solve $(simtracker.iter)" termination_status(control_model) s0

    # try again , either shortening or lengthening the planning window
    solution = nothing
    for i in 1:3
        _sign  = i < 3 ? +1 : -1
        len = sqrt(initial_state.v^2 + initial_state.u^2) * (planning_horizon + (_sign * .01 * i))
        if len < 4
            @warn "Length is $len<4"
            skip=true
            len = 4
        else
            skip = false
        end
        track = trim(FULLTRACK, s0, len)

        _, _, control_model, solution = run_mtm(
            :dynamics,
            1.75;
            track=track,
            icond=initial_state,
            fcond=:minimal,
            showplots=false,
            n_iter=5000,
            quiet=true,
        )

        converged(control_model) && break
        skip && break
    end
    return converged(control_model), solution
end

converged(control_model) = "LOCALLY_SOLVED" == string(termination_status(control_model))
fail() = nothing, nothing, true, nothing

"""
Perform a simulation step
"""
function step(simtracker, globalsolution, planning_horizon::Float64,)
    # get initial conditions
    if simtracker.iter == 1
        initial_state =  solution2state(simtracker.s0, globalsolution)
        simtracker.s = simtracker.s0
        simtracker.t += simtracker.Δt
    else
        initial_state =  solution2state(simtracker.Δt, simtracker.prevsol; at=:time)
        simtracker.s = solution2s(simtracker.Δt, simtracker.prevsol)

        idx = findfirst(simtracker.prevsol.t .>= simtracker.Δt)
        idx = isnothing(idx) ? length(simtracker.prevsol.t) : idx
        simtracker.t += simtracker.prevsol.t[idx]
    end

    # get planning window (trimmed track) & final conditions
    s0 =  max(0.001, simtracker.s)
    if s0 < 250
        len = max(sqrt(initial_state.v^2 + initial_state.u^2) * planning_horizon, 4)
        track = trim(FULLTRACK, s0, len)
        final_state = nothing
    else
        track = trim(FULLTRACK, s0, 200)
        final_state = State(; u=20, ψ=0)
    end

    # (attempt to) solve MTM problem over planning window
    control_model, solution = nothing, nothing
    try
        _, _, control_model, solution = run_mtm(
            :dynamics,
            1.75;
            track=track,
            icond=initial_state,
            fcond=final_state,
            showplots=false,
            n_iter=5000,
            quiet=true,
        )
    catch
        return initial_state, solution, true, track
    end

    # failed -> try alternative strategies
    if !converged(control_model)
        # try to recover a solution by shortening the planning window
        success, solution = attempt_step(simtracker, control_model, s0, initial_state, planning_horizon)
        success || return initial_state, solution, true, track
    end
    simtracker.prevsol = solution

    return initial_state, solution, false, track
end

"""
Run a simulation in which the model can only plan for `planning_horizon` seconds ahead.
"""
function run_simulation(; s0=0.0, sf=258, planning_horizon::Float64=.5, n_iter=1000, Δt=.01)
    # check if a solution was already saved and skip
    name = (@sprintf "s0_%.2f_horizon_length_%.2f" s0 planning_horizon)
    FOLDER = "D:\\Dropbox (UCL)\\Rotation_vte\\Locomotion\\analysis\\behavior\\horizons_mtm_sims"

    destination = joinpath(FOLDER, "$name.csv")
    isfile(destination) && return

    # run global solution
    track = Track(;start_waypoint=4, keep_n_waypoints=-1)

    _, bike, _, globalsolution = run_mtm(
        :dynamics,
        2.5;
        showtrials=nothing,
        track=track,
        n_iter=5000,
        fcond=State(; u=20, ψ=0),
        timed=false,
        showplots=false,
    )

    simtracker = SimTracker(s0=s0, Δt=Δt)
    solutiontracker = SolutionTracker()
    pbar = ProgressBar(; expand=true, columns=:detailed)
    job = addjob!(pbar; N=n_iter, description="Running horizon: $planning_horizon seconds")

    data = nothing
    withpbar(pbar) do
        anim = @animate for i in 1:n_iter
            simtracker.iter = i

            # run simulation and store results
            initial_state, solution, shouldstop, track = step(simtracker, globalsolution, planning_horizon)

            # extend tracked solution with planned solution
            if solution.s[end] >= (sf+5) || shouldstop
                @info "Stopping because step simulation reached the end of the curve or shouldstop"
                extend_with_sol!(solutiontracker, simtracker.prevsol, Δt)
                break
            end

            # extend tracked solution
            add!(solutiontracker, simtracker.t, initial_state, simtracker.s)

            # plot stuff
            draw(:arena)
            draw!(FULLTRACK; border_alpha=.0, alpha=0.0)

            color = "black"  # colors[i]
            draw!(track; color=color, border_lw=5, alpha=0.0)
            plot_bike_trajectory!(solution, bike; showbike=false, label=nothing, color=color, alpha=.6, lw=4)
            draw!(initial_state; color=color, alpha=1)
            plot!(; title="T: $(round(simtracker.t; digits=2))s | horizon: $planning_horizon s")

            simtracker.s > sf && begin
                @info "Stopping because we're beyond sf ($(sf))"
                extend_with_sol!(solutiontracker, simtracker.prevsol, Δt)
                break
            end
            update!(job)
            sleep(0.001)
        end

        if isnothing(simtracker.prevsol)
            @warn "Not saving because empty simulation"
        else
            @info "savin data and animations"
            # gifpath = joinpath(FOLDER, "$name.gif")
            # gif(anim, gifpath, fps=(Int ∘ round)(.25/Δt))

            # save video
            if simtracker.iter > 2
                vidpath = joinpath(FOLDER, "$name.mp4")
                mp4(anim, vidpath, fps=(Int ∘ round)(.25/Δt))
            end

            # save short horizon solution
            data = DataFrame(toDict(solutiontracker))
            CSV.write(destination, data)
        end
    end

    return data
end

horizons = vcat(collect(.08:.01:.38), collect(.38:.02:.54), collect(.5:.05:1.2))
starts = [0.0, 45.0, 98.0, 150.0]
ends = [40.0, 96.0, 144.0, 205.0]

for i in 1:length(starts)
    # if i != 4
    #     continue
    # end
    for horizon in horizons
        @info "Running horizon length $horizon seconds on curve $i"
        results = run_simulation(planning_horizon=horizon, s0=starts[i], sf=ends[i])
    end
end

# results = run_simulation(planning_horizon=.6, s0=0.0, sf=240, Δt=.025)

max_stars_repo_name: FedeClaudi/LocomotionControl

max_issues_repo_name: FedeClaudi/LocomotionControl

max_forks_repo_name: FedeClaudi/LocomotionControl


text:
# convert matrix to band storage
function band(A::AbstractMatrix,kl,ku)
    m, n = size(A)
    AB = zeros(eltype(A),kl+ku+1,n)
    for j = 1:n
        for i = max(1,j-ku):min(m,j+kl)
            AB[ku+1-j+i,j] = A[i,j]
        end
    end
    return AB
end

# convert band storage to general matrix
function unband(AB::AbstractMatrix,m,kl,ku)
    bm, n = size(AB)
    A = zeros(eltype(AB),m,n)
    for j = 1:n
        for i = max(1,j-ku):min(m,j+kl)
            A[i,j] = AB[ku+1-j+i,j]
        end
    end
    return A
end

# zero out elements not on matrix bands
function bandex(A::AbstractMatrix,kl,ku)
    m, n = size(A)
    AB = band(A,kl,ku)
    B = unband(AB,m,kl,ku)
    return B
end

const CublasFloat = Union{Float64,Float32,ComplexF64,ComplexF32}
const CublasReal = Union{Float64,Float32}
const CublasComplex = Union{ComplexF64,ComplexF32}

function Base.convert(::Type{cublasOperation_t}, trans::Char)
    if trans == 'N'
        return CUBLAS_OP_N
    elseif trans == 'T'
        return CUBLAS_OP_T
    elseif trans == 'C'
        return CUBLAS_OP_C
    else
        throw(ArgumentError("Unknown operation $trans"))
    end
end

function Base.convert(::Type{cublasFillMode_t}, uplo::Char)
    if uplo == 'U'
        return CUBLAS_FILL_MODE_UPPER
    elseif uplo == 'L'
        return CUBLAS_FILL_MODE_LOWER
    else
        throw(ArgumentError("Unknown fill mode $uplo"))
    end
end

function Base.convert(::Type{cublasDiagType_t}, diag::Char)
    if diag == 'U'
        return CUBLAS_DIAG_UNIT
    elseif diag == 'N'
        return CUBLAS_DIAG_NON_UNIT
    else
        throw(ArgumentError("Unknown diag mode $diag"))
    end
end

function Base.convert(::Type{cublasSideMode_t}, side::Char)
    if side == 'L'
        return CUBLAS_SIDE_LEFT
    elseif side == 'R'
        return CUBLAS_SIDE_RIGHT
    else
        throw(ArgumentError("Unknown side mode $side"))
    end
end


text:
function comb(n, x)
    s = 0
    for i in 1:n
        for j in i+1:n
            for k in j+1:n
                if i+j+k == x
                    s += 1
                end
            end
        end
    end
    return s
end

while true
    n, x = parse.(Int, split(readline()))
    n == x == 0 && break
    println(comb(n,x))
end


text:
function nn_law(ε)
    ae(ε, [20,20,20,3])
end

function nn_nonassociated_plasticity(epsilon, epsilon_n, sigma)
    coef = ae([epsilon epsilon_n sigma], [20,20,20,6])
    coef = [coef constant(ones(size(coef,1),1))]
    H = consistent_tangent_matrix(coef, Dc)
    return batch_matmul(H, epsilon-epsilon_n) + sigma
end

function nn_associated_plasticity(epsilon, epsilon_n, sigma)

    coef = ae([epsilon epsilon_n sigma], [20,20,20,3])

    coef = [coef coef constant(ones(size(coef,1),1))]
    H = consistent_tangent_matrix(coef, Dc)
    return batch_matmul(H, epsilon-epsilon_n) + sigma
end


text:
using Documenter, Weave
using Yao, YaoBlocks, YaoArrayRegister, YaoBase, BitBasis

# Preprocess weave

# const Examples = ["GHZ", "QFT", "Grover", "QCBM"]
const Examples = ["GHZ", "QFT", "QCBM"]

for each in Examples
    file_path = joinpath(@__DIR__, "src", "examples", join([each, ".jmd"]))
    @info "expanding $file_path to markdown"
    weave(file_path)
end

const PAGES = [
    "Home" => "index.md",
    "Examples" => map(x->joinpath("examples", x * ".md"), Examples),
    "Manual" => Any[
        "man/array_registers.md",
        "man/blocks.md",
        "man/base.md",
        "man/registers.md",
        "man/bitbasis.md",
        "man/extending_blocks.md",
        "man/benchmarks.md",
        ],
]

makedocs(
    modules = [Yao, YaoBase, YaoArrayRegister, YaoBlocks, BitBasis],
    format = Documenter.HTML(
        prettyurls = ("deploy" in ARGS),
        canonical = ("deploy" in ARGS) ? "https://quantumbfs.github.io/Yao.jl/latest/" : nothing,
        assets = ["assets/favicon.ico"],
        ),
    doctest = ("doctest=true" in ARGS),
    clean = false,
    sitename = "Yao.jl",
    linkcheck = !("skiplinks" in ARGS),
    pages = PAGES
)

deploydocs(
    repo = "github.com/QuantumBFS/Yao.jl.git",
    target = "build",
)


text:
export optimize_kahypar, optimize_kahypar_auto
export KaHyParBipartite

"""
    KaHyParBipartite{RT,IT,GM}
    KaHyParBipartite(; sc_target, imbalances=collect(0.0:0.005:0.8),
        max_group_size=40, greedy_config=GreedyMethod())

Optimize the einsum code contraction order using the KaHyPar + Greedy approach.
This program first recursively cuts the tensors into several groups using KaHyPar,
with maximum group size specifed by `max_group_size` and maximum space complexity specified by `sc_target`,
Then finds the contraction order inside each group with the greedy search algorithm. Other arguments are

* `sc_target` is the target space complexity, defined as `log2(number of elements in the largest tensor)`,
* `imbalances` is a KaHyPar parameter that controls the group sizes in hierarchical bipartition,
* `max_group_size` is the maximum size that allowed to used greedy search,
* `greedy_config` is a greedy optimizer.

### References
* [Hyper-optimized tensor network contraction](https://arxiv.org/abs/2002.01935)
* [Simulating the Sycamore quantum supremacy circuits](https://arxiv.org/abs/2103.03074)
"""
Base.@kwdef struct KaHyParBipartite{RT,IT,GM} <: CodeOptimizer
    sc_target::RT
    imbalances::IT = 0.0:0.005:0.8
    max_group_size::Int = 40
    greedy_config::GM = GreedyMethod()
end

function induced_subhypergraph(s::SparseMatrixCSC, group)
    s0 = s[group,:]
    nvs = vec(sum(s0, dims=1))
    remaining_edges = findall(!iszero, nvs)
    s0[:,remaining_edges], remaining_edges
end

function convert2int(sizes::AbstractVector)
    round.(Int, sizes .* 100)
end

function bipartite_sc(bipartiter::KaHyParBipartite, adj::SparseMatrixCSC, vertices, log2_sizes)
    n_v = length(vertices)
    subgraph, remaining_edges = induced_subhypergraph(adj, vertices)
    hypergraph = KaHyPar.HyperGraph(subgraph, ones(n_v), convert2int(log2_sizes[remaining_edges]))
    local parts
    min_sc = 999999
    for imbalance in bipartiter.imbalances
        parts = @suppress KaHyPar.partition(hypergraph, 2; imbalance=imbalance, configuration=:edge_cut)
        part0 = vertices[parts .== 0]
        part1 = vertices[parts .== 1]
        sc0, sc1 = group_sc(adj, part0, log2_sizes), group_sc(adj, part1, log2_sizes)
        sc = max(sc0, sc1)
        min_sc = min(sc, min_sc)
        @debug "imbalance $imbalance: sc = $sc, group = ($(length(part0)), $(length(part1)))"
        if sc <= bipartiter.sc_target
            return part0, part1
        end
    end
    error("fail to find a valid partition for `sc_target = $(bipartiter.sc_target)`, got minimum value `$min_sc` (imbalances = $(bipartiter.imbalances))")
end

# the space complexity (external degree of freedoms) if we contract this group
function group_sc(adj, group, log2_sizes)
    degree_in = sum(adj[group,:], dims=1)
    degree_all = sum(adj, dims=1)
    sum(i->(degree_in[i]!=0 && degree_in[i]!=degree_all[i] ? Float64(log2_sizes[i]) : 0.0), 1:size(adj,2))
end

function bipartition_recursive(bipartiter, adj::SparseMatrixCSC, vertices::AbstractVector{T}, log2_sizes) where T
    if length(vertices) > bipartiter.max_group_size
        parts = bipartite_sc(bipartiter, adj, vertices, log2_sizes)
        groups = Vector{T}[]
        for part in parts
            for component in _connected_components(adj, part)
                push!(groups, component)
            end
        end
        newparts = [bipartition_recursive(bipartiter, adj, groups[i], log2_sizes) for i=1:length(groups)]
        if length(groups) > 2
            tree = coarse_grained_optimize(adj, groups, log2_sizes, bipartiter.greedy_config)
            return map_tree_to_parts(tree, newparts)
        else
            return newparts
        end
    else
        return [vertices]
    end
end

function _connected_components(adj, part::AbstractVector{T}) where T
    A = adj[part,:]
    A = A * A'   # connectivility matrix
    n = length(part)
    visit_mask = zeros(Bool, n)
    groups = Vector{T}[]
    while !all(visit_mask)
        newset = Int[]
        push_connected!(newset, visit_mask, A, findfirst(==(false), visit_mask))
        push!(groups, getindex.(Ref(part), newset))
    end
    return groups
end

function push_connected!(set, visit_mask, adj, i)
    visit_mask[i] = true
    push!(set, i)
    for v = 1:size(adj, 2)
        if !visit_mask[v] && !iszero(adj[i,v])
            push_connected!(set, visit_mask, adj, v)
        end
    end
end

function coarse_grained_optimize(adj, parts, log2_sizes, greedy_config)
    incidence_list = get_coarse_grained_graph(adj, parts)
    log2_edge_sizes = Dict([i=>log2_sizes[i] for i=1:length(log2_sizes)])
    tree, _, _ = OMEinsum.tree_greedy(incidence_list, log2_edge_sizes; method=greedy_config.method, nrepeat=greedy_config.nrepeat)
    return tree
end

function get_coarse_grained_graph(adj, parts)
    ADJ = vcat([sum(adj[part,:], dims=1) for part in parts]...)
    degree_in = sum(ADJ, dims=1)
    degree_all = sum(adj, dims=1)
    openedges = filter(i->degree_in[i]!=0 && degree_in[i]!=degree_all[i], 1:size(adj, 2))
    v2e = Dict{Int,Vector{Int}}()
    for v=1:size(ADJ, 1)
        v2e[v] = findall(!iszero, view(ADJ,v,:))
    end
    incidence_list = OMEinsum.IncidenceList(v2e; openedges=openedges)
    return incidence_list
end

function map_tree_to_parts(tree, parts)
    if tree isa OMEinsum.ContractionTree
        [map_tree_to_parts(tree.left, parts), map_tree_to_parts(tree.right, parts)]
    else
        parts[tree]
    end
end

# KaHyPar
function adjacency_matrix(ixs::AbstractVector)
    rows = Int[]
    cols = Int[]
    edges = unique!([Iterators.flatten(ixs)...])
    for (i,ix) in enumerate(ixs)
        push!(rows, map(x->i, ix)...)
        push!(cols, map(x->findfirst(==(x), edges), ix)...)
    end
    return sparse(rows, cols, ones(Int, length(rows))), edges
end

# legacy interface
"""
    optimize_kahypar(code, size_dict; sc_target, max_group_size=40, imbalances=0.0:0.01:0.2, greedy_method=OMEinsum.MinSpaceOut(), greedy_nrepeat=10)

Optimize the einsum `code` contraction order using the KaHyPar + Greedy approach. `size_dict` is a dictionary that specifies leg dimensions.
Check the docstring of `KaHyParBipartite` for detailed explaination of other input arguments.
"""
function optimize_kahypar(code::EinCode, size_dict; sc_target, max_group_size=40, imbalances=0.0:0.01:0.2, greedy_method=OMEinsum.MinSpaceOut(), greedy_nrepeat=10)
    bipartiter = KaHyParBipartite(; sc_target=sc_target, max_group_size=max_group_size, imbalances=imbalances, greedy_config=GreedyMethod(method=greedy_method, nrepeat=greedy_nrepeat))
    recursive_bipartite_optimize(bipartiter, code, size_dict)
end

function recursive_bipartite_optimize(bipartiter, code::EinCode, size_dict)
    ixs, iy = getixsv(code), getiyv(code)
    ixv = [ixs..., iy]
    adj, edges = adjacency_matrix(ixv)
    vertices=collect(1:length(ixs))
    parts = bipartition_recursive(bipartiter, adj, vertices, [log2(size_dict[e]) for e in edges])
    recursive_construct_nestedeinsum(ixv, iy, parts, size_dict, 0, bipartiter.greedy_config)
end

function recursive_construct_nestedeinsum(ixs::AbstractVector{<:AbstractVector}, iy::AbstractVector{L}, parts::AbstractVector, size_dict, level, greedy_config) where L
    if length(parts) == 2
        # code is a nested einsum
        code1 = recursive_construct_nestedeinsum(ixs, iy, parts[1], size_dict, level+1, greedy_config)
        code2 = recursive_construct_nestedeinsum(ixs, iy, parts[2], size_dict, level+1, greedy_config)
        AB = recursive_flatten(parts[2]) ∪ recursive_flatten(parts[1])
        inset12, outset12 = ixs[AB], ixs[setdiff(1:length(ixs), AB)]
        iy12 = Iterators.flatten(inset12) ∩  (Iterators.flatten(outset12) ∪ iy)
        iy1, iy2 = getiyv(code1.eins), getiyv(code2.eins)
        return NestedEinsum([code1, code2], EinCode([iy1, iy2], L[(level==0 ? iy : iy12)...]))
    elseif length(parts) == 1
        return recursive_construct_nestedeinsum(ixs, iy, parts[1], size_dict, level, greedy_config)
    else
        error("not a bipartition, got size $(length(parts))")
    end
end

function recursive_construct_nestedeinsum(ixs::AbstractVector{<:AbstractVector}, iy::AbstractVector{L}, parts::AbstractVector{<:Integer}, size_dict, level, greedy_config) where L
    if isempty(parts)
        error("got empty group!")
    end
    inset, outset = ixs[parts], ixs[setdiff(1:length(ixs), parts)]
    iy1 = Iterators.flatten(inset) ∩  (Iterators.flatten(outset) ∪ iy)
    res = optimize_greedy(DynamicEinCode{L}, inset, iy1, size_dict; method=greedy_config.method, nrepeat=greedy_config.nrepeat)
    return maplocs(res, parts)
end

maplocs(ne::NestedEinsum{ET}, parts) where ET = isleaf(ne) ? NestedEinsum{ET}(parts[ne.tensorindex]) : NestedEinsum(maplocs.(ne.args, Ref(parts)), ne.eins)

function kahypar_recursive(ne::NestedEinsum; log2_size_dict, sc_target, min_size, imbalances=0.0:0.04:0.8)
    if length(ne.args >= min_size) && all(isleaf, ne.args)
        bipartite_eincode(adj, ne.args, ne.eins; log2_size_dict=log2_size_dict, sc_target=sc_target, min_size=min_size, imbalances=imbalances)
    end
    kahypar_recursive(ne.args; log2_size_dict, sc_target=sc_target, min_size=min_size, imbalances=imbalances)
end

recursive_flatten(obj::Tuple) = vcat(recursive_flatten.(obj)...)
recursive_flatten(obj::AbstractVector) = vcat(recursive_flatten.(obj)...)
recursive_flatten(obj) = obj

"""
    optimize_kahypar_auto(code, size_dict; max_group_size=40, greedy_method=OMEinsum.MinSpaceOut(), greedy_nrepeat=10)

Find the optimal contraction order automatically by determining the `sc_target` with bisection.
It can fail if the tree width of your graph is larger than `100`.
"""
function optimize_kahypar_auto(code::EinCode, size_dict; max_group_size=40, effort=500, greedy_method=OMEinsum.MinSpaceOut(), greedy_nrepeat=10)
    sc_high = 100
    sc_low = 1
    order_high = optimize_kahypar(code, size_dict; sc_target=sc_high, max_group_size=max_group_size, imbalances=0.0:0.6/effort*(sc_high-sc_low):0.6)
    _optimize_kahypar_auto(code, size_dict, sc_high, order_high, sc_low, max_group_size, effort, greedy_method, greedy_nrepeat)
end
function _optimize_kahypar_auto(code::EinCode, size_dict, sc_high, order_high, sc_low, max_group_size, effort, greedy_method, greedy_nrepeat)
    if sc_high <= sc_low + 1
        order_high
    else
        sc_mid = (sc_high + sc_low) ÷ 2
        try
            order_mid = optimize_kahypar(code, size_dict; sc_target=sc_mid, max_group_size=max_group_size, imbalances=0.0:0.6/effort*(sc_high-sc_low):0.6)
            order_high, sc_high = order_mid, sc_mid
            # `sc_target` too high
        catch
            # `sc_target` too low
            sc_low = sc_mid
        end
        _optimize_kahypar_auto(code, size_dict, sc_high, order_high, sc_low, max_group_size, effort, greedy_method, greedy_nrepeat)
    end
end


text:
using Documenter, SeisModels

makedocs(
    sitename = "SeisModels.jl documentation",
    pages = [
        "Home" => "index.md",
        "Inbuilt models" => "inbuilt_models.md",
        "Model types" => "types.md",
        "Examples" => "examples.md",
        "Model input-output" => "io.md",
        "Function index" => "function_index.md",
        "References" => "references.md"
        ]
    )

deploydocs(
    repo = "github.com/anowacki/SeisModels.jl.git",
)


text:
@system Ratio begin
    primordia ~ hold

    shoot_to_root_ratio => 0.7 ~ preserve

    root_to_shoot_ratio(shoot_to_root_ratio) => (1 - shoot_to_root_ratio) ~ preserve

    leaf_to_stem_ratio => 0.9 ~ preserve

    stem_to_leaf_ratio(leaf_to_stem_ratio) => (1 - leaf_to_stem_ratio) ~ preserve

    initial_leaf_ratio(shoot_to_root_ratio, leaf_to_stem_ratio, primordia) => begin
        #TODO how to handle primordia?
        shoot_to_root_ratio * leaf_to_stem_ratio / primordia
    end ~ preserve
end