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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 8b184833-59d9-4958-8460-369b1ea19b9e
using Printf

# ╔═╡ 0cac30f1-9101-4ba3-accc-52621bc1d16f
begin
	using PlutoUI, PlutoTest, PlutoTeachingTools
	using BenchmarkTools
	using Profile,  ProfileSVG, FlameGraphs
	using Plots # just needed for colorant when color FlameGraphs
	using LinearAlgebra, PDMats # used by periodogram.jl
	using Statistics: mean, std    # used by periodogram.jl
	using Random
	Random.seed!(123)
	eval(Meta.parse(code_for_check_type_funcs))
end;

# ╔═╡ f37ba7bb-0484-4a53-bc8b-94b06d685ac0
md"""
# Astro 528, Lab 5, Exercise 1
# Profiling & Code Inspection
"""

# ╔═╡ 9e04db5d-0662-4a39-af3b-85e68f948ebc
md"""
In this exercise, we'll profile some functions that we've used in [Lab 4](https://psuastro528.github.io/lab4-start/): `calc_ecc_anom`, `calc_rv_keplerian`, and  `calc_periodogram`.
You may want to look at the source code and/or the help information for each of these functions to remind yourself of what they do.  In order to make the profiler outpus easier to read, the functions that we'll profile have been placed into separate files inside this repository's `src` directory.
"""

# ╔═╡ f80f430b-45d0-4b2a-8be6-e39aabccd6ae
md"## Example with efficient code"

# ╔═╡ 7fcc3f2a-cf98-4a95-a52e-518ca880f31d
md"""
First, we'll demonstrate the syntax for how to load code from another file into a Pluto notebook.
"""

# ╔═╡ 0d25c024-ab14-44bb-85bf-15a86152e7a1
begin
	KeplerEqn = ingredients("./src/calc_rv.jl")
	import .KeplerEqn: calc_ecc_anom, calc_rv_keplerian
end

# ╔═╡ d1fae43a-b6a7-4ae2-aa5c-e495d59f5e77
md"""
Next, we'll demonstrate the syntax for using Julia's profiler.
"""

# ╔═╡ 1e4df7eb-d8de-4bae-912e-c2360f530b76
with_terminal() do
	calc_ecc_anom    # Tell Pluto to make sure our function is avaliable
	Profile.clear()  # Clear the data stored by the profiler
	@profile calc_ecc_anom.(π/4,0.5)
	Profile.print()
end

# ╔═╡ 3f55c815-65a2-4950-8b2f-bc6d19e606c6
md"""
The odds are that you got a warning message that there were no samples collected.  Julia uses a *statistical* profiler, meaning that it doesn't count every time a line is executed, but rather it checks what line of code is being executed every so often.  `calc_ecc_anom` returned so quickly that the profiler didn't get a chance to figure out where it was spending its time.  In this case, it's so fast that we can't make the delay between samples fast enough, so we'll combine both of the suggestions from the warning message to get usable profiling data.
"""

# ╔═╡ f5c19252-5c37-4688-b1f8-e29afa006684
begin
	calc_ecc_anom    # Tell Pluto to make sure our function is avaliable
	Profile.clear()  # Clear the data stored by the profiler
	Profile.init(delay=1/10^7)    # Check what code is running every 0.1μs
	for i in 1:10^5               # Accumulate data over multiple function calls
		mean_anom = 2π*rand()     # Generate random input values
		ecc = rand()
		# Actually run profiler
		@profile calc_ecc_anom(mean_anom,ecc)
	end
	# Explicilty store profiling results to avoid confusion later on
	retrieve_prof_calc_ecc_anom = Profile.retrieve()
end;

# ╔═╡ cce0518b-68fc-43a7-866a-d02febcf9901
with_terminal() do
	retrieve_prof_calc_ecc_anom
	Profile.print(retrieve_prof_calc_ecc_anom...)
end

# ╔═╡ fe3931f6-11cf-4853-8bc9-d621a3af29f9
md"""
Ok, now you should see lots of data.  Each line of output corresponds to one line of code.  On the right of each line is the file containing the relevant code, the line number and a  function.  The number immediately to the left of the text is the number of times that the profiles checked in and found that the computer was executing that line of code *including* any computations in functions that resulted from it.

This is a *tree view*, meaning that the lines are organized by what function calls what function.  The first several lines are because we're running inside a Pluto notebook.  Pluto wraps our cells in several functions to make it's reactive environment work.

Skip down to the line just below the line containing 'Profile.jl:28; macro expansion'.  That's where the work we're interested begins.  The number immediately to the left is the total number of profiler samples collected for our code (in this case `calc_ecc_anom`).

1a.  How many samples did you get for 'kepler\_eqn.jl:46; calc\_ecc\_anom'?
"""

# ╔═╡ 8894cd1c-5f94-4b07-94a0-757b2fdf75db
response_1a = missing # Replace with an integer

# ╔═╡ 640647b7-1388-41dd-9a9a-e440898be3ab
display_msg_if_fail(check_type_isa(:response_1a,response_1a,Integer))

# ╔═╡ 3005d8e2-f018-43eb-9759-8ee48fa9ea1a
md"""
1b.  Scroll down looking for other big numbers (e.g., at least a quarter of $response_1a).
Ignore any lines that refer to 'Profile.jl' or 'task.jl'.
Which line numer of code in kepler_eqn is taking the most time?
"""

# ╔═╡ 4fb20f7c-ba52-4381-b75d-1a726e62569c
response_1b = missing # Replace with an integer

# ╔═╡ 4a37d8c9-8766-4692-8c8d-4f14c67df4ee
begin
	if !@isdefined(response_1b)
		var_not_defined(:response_1b)
	elseif ismissing(response_1b)
		still_missing()
	elseif response_1b != 54
		keep_working(md"Please double check that.")
	else
		correct()
	end
end

# ╔═╡ 9d901dca-4660-43c8-9781-a87b02fd54e3
md"""
We can display the profiling results in other formats. For example, instead of the tree view, we can use a *flat view* and sort by the number of samples for that line of code.
"""

# ╔═╡ d1f872d2-0953-47d1-b4cc-e177496677b0
with_terminal() do
	Profile.print(retrieve_prof_calc_ecc_anom..., format=:flat, sortedby=:count)
end

# ╔═╡ dc0dce98-3824-44ef-bd58-d14c0a93f8ac
md"""
In this format, the first column contains the number of samples, the third colum contains the filename, the fourth column is the line number and the fifth column is the function within that line.

1c.  Scroll to the bottom of the above output, then go back up to find the first time (you see 'calc\_ecc\_anom' (starting from the bottom).  Look at the several lines above it (i.e., the lines of code account for a significant fraction of the total cost).  What line of code in 'update\_ecc\_anom\_laguerre' is taking the most time?
"""

# ╔═╡ ef76ab85-a2bd-4b5b-acc1-176eead528f5
response_1c = missing  # Replace with an integer

# ╔═╡ 23fbae8c-cf91-4946-b31d-86cd9c2acc65
begin
	if !@isdefined(response_1c)
		var_not_defined(:response_1c)
	elseif ismissing(response_1c)
		still_missing()
	elseif response_1c != 27
		keep_working(md"Please double check that.")
	else
		correct()
	end
end

# ╔═╡ 2869aaa6-9941-46b0-b240-cf9e39231525
md"1d. Look at the source code for 'update\_ecc\_anom\_laguerre'.  What is that line doing?"

# ╔═╡ 4ed00d92-eff3-4ff5-a97d-c844e7a48c67
response_1d = missing # Replace with md"Response"

# ╔═╡ 779d6fb0-4e6c-4e4b-9720-f26b0b77990d
display_msg_if_fail(check_type_isa(:response_1d,response_1d,Markdown.MD))

# ╔═╡ 182dc7ba-bb2d-4b95-b393-79b0860dd8d0
md"""
This is an example of a funciton with very little room for improvement.  Most of the time is being spent on computing sines, cosines, square roots and basic arithmetic.  Notably, we do *not* see significant time being devoted to memory allocations.  So it's unlikely that we'll be able to optimize this function further (unless we were to find a more efficient way to solve the problem by using a different algorithm).

Before we move on, we'll also try looking at a visual representation of the same data.
"""

# ╔═╡ 54fd03d0-9a90-429d-af80-824196cd40e6
begin  # dispaly configuration paramteres for ProfileSVG.view
	svg_fontsize = 14
	svg_width = 706
end;

# ╔═╡ 3a961a8b-7b61-48e8-b840-8ef6c5acafe8
ProfileSVG.view(data=retrieve_prof_calc_ecc_anom,fontsize=svg_fontsize,width=svg_width)

# ╔═╡ 049f048b-12ec-4ab3-891b-257c8b766a10
md"""
When interpretting such graphs, it can be helpful to specify how the cells will be color coded.  Here, we'll demonstrate how to color code cells by what module they're defined in.
"""

# ╔═╡ be9e0f8d-a424-4025-bb1e-baa876664633
md"""
Look in the tall part of the stack for the cell 'kepler\_eqn.jl:46; calc\_ecc\_anom' (just below the first lowest green cell).  Double click on that cell to zoom in on the time spent inside 'calc\_ecc\_anom'.  You can hover your mouse over other cells to see what line they refer to.
"""

# ╔═╡ 713265a1-2431-4e20-8e3d-965c924be1d1
md"""
1e. Which of the following operations is taking more of the computer's time?"""

# ╔═╡ c23f1cfb-de99-4119-ba00-14af073e33f1
@bind response_1e Radio(["subtraction","multiplication","division","trig functions","sqrt"])

# ╔═╡ 1c217bf9-5c67-42d7-8be7-b8b635a1ed93
begin
	if !@isdefined(response_1e)
		var_not_defined(:response_1e)
	elseif ismissing(response_1e)
		still_missing(md"Remember to respond to 1e.")
	elseif response_1e != "trig functions"
		keep_working(md"Please look again.")
	else
		correct()
	end
end

# ╔═╡ 6b5c76d9-da97-49a4-9d2c-5ae22e08b44a
md"## Example with an opportunity to improve efficiency"

# ╔═╡ e76ff922-f61f-4c80-836b-9cbbb0952d2d
md"""
Now, we'll consider two versions of functions to calculate a periodogram.  We'll use the profiler to help us find an opportunity to improve the efficiency of our code. First, we'll load two versions of the code into two modules named PeriodogramOrig (for the original code) and Periodogram (for the improved code).  This pattern of placing two versions of code to do the same thing in two different files can be very useful when optimizing code.  It allows you to keep the same function names, rather than changing them (and having to remember to update things consistently).  And it makes it easy to test that they give consistent results.
"""

# ╔═╡ b3d3f4c3-968d-4235-82b9-d6f9806109f8
begin
	PeriodogramOrig = ingredients("src/periodogram_orig.jl")
	calc_periodogram_orig = PeriodogramOrig.calc_periodogram
	Periodogram = ingredients("src/periodogram.jl")
	import .Periodogram: calc_periodogram
end

# ╔═╡ 3973c872-5817-4215-925a-b93e2186e33a
function color_by_module(mod::Module)
	pm = FlameGraphs.pm
	(@isdefined Periodogram) && (mod === Periodogram || pm(mod) === Periodogram) && return colorant"purple"
	(@isdefined PeriodogramOrig) && (mod === PeriodogramOrig || pm(mod) === PeriodogramOrig) && return colorant"purple"
	(@isdefined KeplerEqn) && (mod === KeplerEqn || pm(mod) === KeplerEqn) && return colorant"green"
	(@isdefined PDMats) && (mod === PDMats || pm(mod) === PDMats) && return colorant"yellow"
	(@isdefined LinearAlgebra) && (mod === LinearAlgebra || pm(mod) === LinearAlgebra) && return colorant"cyan"
	(mod === PlutoRunner || pm(mod) === PlutoRunner) && return colorant"grey40"
    (mod === Core.Compiler || pm(mod) === Core.Compiler) && return colorant"gray60"
    (mod === Core || pm(mod) === Core) && return colorant"gray20"
    (mod === Base || pm(mod) === Base) && return colorant"lightblue"
    return nothing
end

# ╔═╡ 78686ae5-4219-43b5-896b-3a86864ce049
ProfileSVG.view(data=retrieve_prof_calc_ecc_anom,StackFrameCategory(color_by_module))

# ╔═╡ 71205282-e4bc-479f-b7c0-db6783066d84
md"### Generate data to analyze"

# ╔═╡ 938e338e-9242-4343-accf-c9d5d1c0420a
md"""
Now, we'll generate a set of Keplerian orbital parameters and simulated data to use for profiling our code.
"""

# ╔═╡ 104828fd-61d1-48af-879c-719b95cc1074
begin
	P = 100.0
	K = 10.0
	ecc = rand()
	ω = 2π*rand()
	M0 = 2π*rand()
	param_true = [P,K,ecc,ω,M0]
end;

# ╔═╡ ffd9315b-d5ad-4636-ba2f-5380c22055d2
num_obs = 100;

# ╔═╡ b5f4992a-5233-4ade-8ca4-bf2e6da1a4ac
num_periods = 10000;

# ╔═╡ e36ebf1c-a750-40b4-8815-8bb3dee3910f
begin
	days_in_year = 365.2425
	times = 2*days_in_year*rand(num_obs)
	rvs_true = calc_rv_keplerian.(times,P,K,ecc,ω,M0)
	σ_rvs = 1*ones(num_obs)
	rvs_obs = rvs_true .+ σ_rvs .* randn(num_obs)
end;

# ╔═╡ 7cba065d-98b4-4e04-a3b4-76224878a0bd
md"### Profile original periodogram code"

# ╔═╡ c97e19e9-e818-432e-aaeb-618090372fad
if true
	calc_periodogram_orig(times,rvs_obs,σ_rvs,num_periods=num_periods)   # Make sure compiled before profiling
	Profile.clear()
	Profile.init(delay=1/10^5)
	for i in 1:10
		@profile calc_periodogram_orig(times,rvs_obs,σ_rvs,num_periods=num_periods)
	end
	# Explicitly store profiling results to avoid confusion later on
	retrieve_prof_periodogram_orig = Profile.retrieve()
end;

# ╔═╡ e54ae808-22f9-4ee5-9b75-d93aa56a8b3d
ProfileSVG.view(data=retrieve_prof_periodogram_orig,fontsize=svg_fontsize,width=svg_width,StackFrameCategory(color_by_module))

# ╔═╡ 74c3f2cb-6865-4cde-9c61-b75777854866

# ╔═╡ 488e83b4-4f5e-400f-85e2-f3b0716d7404
md"""
As before, there's about a dozen levels of functions at the bottom of the graph before we get to our function of interest, 'calc\_periodogram'.  Double click on the cell with 'calc\_periodogram' and line 31.  Look at how the code is spending its time.

1f.  Look for some cells near the top that take a non-trivial fraction of the time and aren't doing math.  What is taking time and could be avoided (or at least significantly reduced)?
"""

# ╔═╡ 5e48f2ae-bfc8-4130-a831-fae02b331bf7
response_1f = missing  # Replace with md"Response"

# ╔═╡ 31f2d123-8a9a-45c6-9f80-569b070c54e4
display_msg_if_fail(check_type_isa(:response_1f,response_1f,Markdown.MD))

# ╔═╡ 990f6bc0-f8d0-4374-b28c-82664e24984a
hint(md"The function `similar` creates an uninitialized array with the same size and element type as its argument.  The function `Array` allocates memory for an array.")

# ╔═╡ e5d41d8f-3b62-434b-9491-98a017e90658
md"### Updated perodogram code"

# ╔═╡ 40d92503-c7ab-4bf7-af62-02b08a284c17
md"""
Based on analyzing the above profiling information, I decided to pre-allocate memory for the [design matrix](https://en.wikipedia.org/wiki/Design_matrix) used for the generalized linear least squares regression (in the case of a periodogram the design matrix contains the values of sin(2π×t/P) and cos(2π×t/P) at each observation time for the putative orbital period being considered).  After implementing that, I saw that there was still significant time spend allocating memory.  Therefore, I preallocated a total of three matrices and four vectors.  I updated the periodogram code to use *in-place* functions like [`mul!`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.mul!), [`ldiv!`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.ldiv!) and [`lu!`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.lu!) so that calculations were written into an existing array/matrix, rather than allocating new memory for the result.
"""

# ╔═╡ 337d9aaf-e2f5-49cc-bd68-178d12110390
md"""
Before looking into any performance improvements, we should double check that we're still getting the same results.
"""

# ╔═╡ 2b93af3a-066b-443f-88c4-702b17e555b7
 begin
	result_orig = calc_periodogram_orig(times,rvs_obs,σ_rvs)
	result_new  = calc_periodogram(times,rvs_obs,σ_rvs)
	@test all(result_new.periodogram .== result_orig.periodogram) && all(result_new.predict .== result_orig.predict)
end

# ╔═╡ e16a4171-e82d-4c43-af26-1f183fd15b2f
md"""
Since we pre-allocated memory in order to reduce the memory being allocated, let's check to see if/how much of a difference our modifications made.
"""

# ╔═╡ e6ae2efe-d519-4c44-b706-a0ca95e6c050
calc_periodogram

# ╔═╡ 9de88c42-a73b-42f1-95be-e117b5beb77b
begin
	(calc_periodogram_orig, calc_periodogram )
	(times,rvs_obs,σ_rvs,num_periods)
	mem_pgram_orig = @allocated calc_periodogram_orig(times,rvs_obs,σ_rvs,num_periods=num_periods)
	mem_pgram_new  = @allocated calc_periodogram(times,rvs_obs,σ_rvs,num_periods=num_periods)
	mem_pgram_percent_str = @sprintf "%4.1f%c" mem_pgram_new/mem_pgram_orig*100 '%'
end;

# ╔═╡ 11bd6e51-6993-4b8e-945c-7e61ff970dd4
if @isdefined mem_pgram_percent_str
	md"""
	The new version only allocates $mem_pgram_percent_str as much memory as the original version!
	"""
end

# ╔═╡ cc4b05de-6d6a-4c61-bdca-a0f6992ecfca
md"Now, let's profile the updated periodogram code."

# ╔═╡ 2facd1bf-e9ac-46c1-9523-1244395a3e8b
if true
	calc_periodogram
	calc_periodogram(times,rvs_obs,σ_rvs,num_periods=num_periods)   # Make sure compiled before profiling
	Profile.clear()
	Profile.init(delay=1/10^5)
	for i in 1:10
		@profile calc_periodogram(times,rvs_obs,σ_rvs,num_periods=num_periods)
	end

	# Explicitly store profiling results to avoid confusion later on
	retrieve_prof_periodogram = Profile.retrieve()
end;

# ╔═╡ 63fbadee-2ad4-4f68-9ba1-5c9d0143eb1a
if @isdefined retrieve_prof_periodogram
	ProfileSVG.view(data=retrieve_prof_periodogram[1],
			fontsize=svg_fontsize,width=svg_width,StackFrameCategory(color_by_module))
end

# ╔═╡ 1d9033dd-e266-4192-930d-ef5be6780226
md"""
Double click on the cell with 'calc\_periodogram' and line 31.  Look at how the code is spending its time.

1g.  Do you notice any memory allocations inside calc_periodogram?  How does the time spent allocating memory compare to the original periodogram code?
"""

# ╔═╡ dba75149-3e4c-4046-8fc7-42e8bf7375a5
response_1g = missing  # Replace with md"Response"

# ╔═╡ a16b0dc3-56c4-4f1c-94f2-ab034bfc30f3
display_msg_if_fail(check_type_isa(:response_1g,response_1g,Markdown.MD))

# ╔═╡ e3fb296e-c338-460f-b73f-fc5a4ffe2368
md"Now, let's benchmark our original and updated periodogram functions."

# ╔═╡ c85762f2-6ad0-43aa-99a2-abd1e8249cdd
md"**Original Periodogram code**"

# ╔═╡ d30212db-5306-4d3c-a786-f568ca67725f
@benchmark calc_periodogram_orig($times,$rvs_obs,$σ_rvs,num_periods=$num_periods)

# ╔═╡ 9be239b8-6c91-47c7-97c6-c39cddd3bc5d
md"**Updated Periodogram code**"

# ╔═╡ 2c0adb00-71b6-483c-be10-f31642690184
(@isdefined calc_periodogram) &&
 @benchmark calc_periodogram($times,$rvs_obs,$σ_rvs,num_periods=$num_periods)

# ╔═╡ 7557b534-c348-44bd-8de2-557a290e1b59
md"1h.  How does the run time compare?"

# ╔═╡ 5d78bd1f-47b3-422c-ae34-11f5e2ae93c3
response_1h = missing # Replace with md"Response"

# ╔═╡ 5514f1a7-bc64-4042-8c5e-b3dc835ebd2d
display_msg_if_fail(check_type_isa(:response_1h,response_1h,Markdown.MD))

# ╔═╡ 0bf114a2-ec3a-4625-b60d-0287f447d2a4
md"For larger codes, sometimes the results are so big (or there is such a deep tree of functions), that you want the profiler output to be spread over more columns than the default.  Below, we'll see how to write the results to a file with a width of our choice."

# ╔═╡ 3d081d57-0bdb-448b-801f-b1a3bb8deb67
begin
	prof_periodogram_filename = "calc_periodogram.prof"
	open(prof_periodogram_filename, "w") do s
    	Profile.print(IOContext(s, :displaysize => (24, 500)),retrieve_prof_periodogram..., format=:tree )
	end
end

# ╔═╡ 5d715c51-147f-4a98-93cd-0c5fc2c51430
md"""
Once we've written the profiling results to a file, we can read each line and select the lines that refer to specific functions that we're specifically interested in.
"""

# ╔═╡ 3c9822b6-2643-4b0d-be1c-4716d701fbfa
prof_periodogram_lines = readlines(prof_periodogram_filename);

# ╔═╡ d253af17-92a7-4c62-a227-fb8cf9298939
prof_periodogram_lines[occursin.("calc_periodogram(",prof_periodogram_lines)]

# ╔═╡ 36de15f2-100d-4e21-8ca1-6a4b3a08a52d
prof_periodogram_lines[occursin.(r"; Array$",prof_periodogram_lines)]

# ╔═╡ ed673d82-55eb-4ba8-9e64-410ecc64e499
md"1i.  Based on the above results (be sure to click the triangles to expand the array of strings), estimate by the maximum percentage we could further speed up  `calc_periodogram`, if we were able to completely eliminate time spent allocating arrays."

# ╔═╡ 3d117fff-2447-4b3b-ab71-301b38e24d54
response_1i = missing  # Replace with md"Response"

# ╔═╡ 8a786af1-9902-4f58-ba58-4cf9dd3b9d3e
display_msg_if_fail(check_type_isa(:response_1i,response_1i,Markdown.MD))

# ╔═╡ 71ccf51e-2d8c-46d1-8b87-953bcbf8302f
md"Once you have a working serial implementation of your class project code, you'll wan to run a profiler on it to determine what are the most time consuming parts of the code.  That will help you identify what portions have the potential to give you a significant speed up."

# ╔═╡ 85ea82e9-35c6-41be-a575-ab823809eef7
md"Here's a good place to stop and move on to exercise 2.  If you have some extra time, then you could consider the following questions."

# ╔═╡ 10403496-5a93-4890-b399-93d6086ff3ee
md"## If you have time..."

# ╔═╡ 2ba5f81f-019b-4749-a0a2-35ead50dee5f
md"1j.  Do you think the difference in performance will become more or less significant if you increase 'num_periods' (the number of putative orbital periods searched)?  Try it.  Report what you find and try to explain your findings.
"

# ╔═╡ 1cd04b3f-3a4e-4c93-85f0-4a41773c4366
response_1j = missing  # Replace with md"Response"

# ╔═╡ 9bb894f5-6247-43b0-a632-0a034be4abe8
md"1k.  Do you think the difference in performance will become more or less significant if you increase 'num_obs' (the number of radial velocity observations)?  Try it.  Report what you find and try to explain your findings."

# ╔═╡ ca2acf91-28d1-4ddc-8150-f4eecef3dffa
response_1k = missing  # Replace with md"Response"

# ╔═╡ 2871b32a-180e-4d57-a1c8-44bc68ec1816
md"1l.  Can you think of a way to speed up the `calc_periodogram` by modifying the algorithm (to calculate the same values, perhaps with slightly different values due to round-off issues)?"

# ╔═╡ bc1e3f75-8d92-4311-877a-6a4d7f227830
response_1l = missing  # Replace with md"Response"

# ╔═╡ b0eccdc2-40bc-44d7-8a72-3e059a4214ab
md"# Helper Code"

# ╔═╡ 81d7a663-1079-4613-afd8-de95eb73f8ae
ChooseDisplayMode()

# ╔═╡ d93a8dca-9c63-4551-87e8-9939a765bef1
TableOfContents(aside=true)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
FlameGraphs = "08572546-2f56-4bcf-ba4e-bab62c3a3f89"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PlutoTeachingTools = "661c6b06-c737-4d37-b85c-46df65de6f69"
PlutoTest = "cb4044da-4d16-4ffa-a6a3-8cad7f73ebdc"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
Profile = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
ProfileSVG = "132c30aa-f267-4189-9183-c8a63c7e05e6"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[compat]
BenchmarkTools = "~1.1.4"
FlameGraphs = "~0.2.5"
PDMats = "~0.11.1"
Plots = "~1.22.1"
PlutoTeachingTools = "~0.1.4"
PlutoTest = "~0.1.0"
PlutoUI = "~0.7.9"
ProfileSVG = "~0.2.1"
"""

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

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

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

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

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

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

[[BenchmarkTools]]
deps = ["JSON", "Logging", "Printf", "Statistics", "UUIDs"]
git-tree-sha1 = "42ac5e523869a84eac9669eaceed9e4aa0e1587b"
uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
version = "1.1.4"

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

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

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

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

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

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "1a90210acd935f222ea19657f143004d2c2a1117"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.38.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"

[[FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "3c041d2ac0a52a12a27af2782b34900d9c3ee68c"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.11.1"

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

[[FlameGraphs]]
deps = ["AbstractTrees", "Colors", "FileIO", "FixedPointNumbers", "IndirectArrays", "LeftChildRightSiblingTrees", "Profile"]
git-tree-sha1 = "99c43a8765095efa6ef76233d44a89e68073bd10"
uuid = "08572546-2f56-4bcf-ba4e-bab62c3a3f89"
version = "0.2.5"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[HypertextLiteral]]
git-tree-sha1 = "1e3ccdc7a6f7b577623028e0095479f4727d8ec1"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.8.0"

[[IndirectArrays]]
git-tree-sha1 = "c2a145a145dc03a7620af1444e0264ef907bd44f"
uuid = "9b13fd28-a010-5f03-acff-a1bbcff69959"
version = "0.5.1"

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

[[LeftChildRightSiblingTrees]]
deps = ["AbstractTrees"]
git-tree-sha1 = "71be1eb5ad19cb4f61fa8c73395c0338fd092ae0"
uuid = "1d6d02ad-be62-4b6b-8a6d-2f90e265016e"
version = "0.1.2"

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

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

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "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 = "457b13497a3ea4deb33d273a6a5ea15c25c0ebd9"
uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
version = "1.22.2"

[[PlutoTeachingTools]]
deps = ["LaTeXStrings", "Markdown", "PlutoUI", "Random"]
git-tree-sha1 = "e2b63ee022e0b20f43fcd15cda3a9047f449e3b4"
uuid = "661c6b06-c737-4d37-b85c-46df65de6f69"
version = "0.1.4"

[[PlutoTest]]
deps = ["HypertextLiteral", "InteractiveUtils", "Markdown", "Test"]
git-tree-sha1 = "ada2eae88798ed6c93d9acb5e41e1671794bb8c8"
uuid = "cb4044da-4d16-4ffa-a6a3-8cad7f73ebdc"
version = "0.1.1"

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

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

[[Profile]]
deps = ["Printf"]
uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"

[[ProfileSVG]]
deps = ["Colors", "FlameGraphs", "Profile", "UUIDs"]
git-tree-sha1 = "e4df82a5dadc26736f106f8d7fc97c42cc6c91ae"
uuid = "132c30aa-f267-4189-9183-c8a63c7e05e6"
version = "0.2.1"

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

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

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

[[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:
# ╟─f37ba7bb-0484-4a53-bc8b-94b06d685ac0
# ╟─9e04db5d-0662-4a39-af3b-85e68f948ebc
# ╟─f80f430b-45d0-4b2a-8be6-e39aabccd6ae
# ╟─7fcc3f2a-cf98-4a95-a52e-518ca880f31d
# ╠═0d25c024-ab14-44bb-85bf-15a86152e7a1
# ╟─d1fae43a-b6a7-4ae2-aa5c-e495d59f5e77
# ╠═1e4df7eb-d8de-4bae-912e-c2360f530b76
# ╟─3f55c815-65a2-4950-8b2f-bc6d19e606c6
# ╠═f5c19252-5c37-4688-b1f8-e29afa006684
# ╠═cce0518b-68fc-43a7-866a-d02febcf9901
# ╟─fe3931f6-11cf-4853-8bc9-d621a3af29f9
# ╠═8894cd1c-5f94-4b07-94a0-757b2fdf75db
# ╟─640647b7-1388-41dd-9a9a-e440898be3ab
# ╟─3005d8e2-f018-43eb-9759-8ee48fa9ea1a
# ╠═4fb20f7c-ba52-4381-b75d-1a726e62569c
# ╟─4a37d8c9-8766-4692-8c8d-4f14c67df4ee
# ╟─9d901dca-4660-43c8-9781-a87b02fd54e3
# ╠═d1f872d2-0953-47d1-b4cc-e177496677b0
# ╟─dc0dce98-3824-44ef-bd58-d14c0a93f8ac
# ╠═ef76ab85-a2bd-4b5b-acc1-176eead528f5
# ╟─23fbae8c-cf91-4946-b31d-86cd9c2acc65
# ╟─2869aaa6-9941-46b0-b240-cf9e39231525
# ╠═4ed00d92-eff3-4ff5-a97d-c844e7a48c67
# ╟─779d6fb0-4e6c-4e4b-9720-f26b0b77990d
# ╟─182dc7ba-bb2d-4b95-b393-79b0860dd8d0
# ╟─54fd03d0-9a90-429d-af80-824196cd40e6
# ╠═3a961a8b-7b61-48e8-b840-8ef6c5acafe8
# ╟─049f048b-12ec-4ab3-891b-257c8b766a10
# ╠═3973c872-5817-4215-925a-b93e2186e33a
# ╠═78686ae5-4219-43b5-896b-3a86864ce049
# ╟─be9e0f8d-a424-4025-bb1e-baa876664633
# ╟─713265a1-2431-4e20-8e3d-965c924be1d1
# ╟─c23f1cfb-de99-4119-ba00-14af073e33f1
# ╟─1c217bf9-5c67-42d7-8be7-b8b635a1ed93
# ╟─6b5c76d9-da97-49a4-9d2c-5ae22e08b44a
# ╟─e76ff922-f61f-4c80-836b-9cbbb0952d2d
# ╠═b3d3f4c3-968d-4235-82b9-d6f9806109f8
# ╟─71205282-e4bc-479f-b7c0-db6783066d84
# ╟─938e338e-9242-4343-accf-c9d5d1c0420a
# ╠═104828fd-61d1-48af-879c-719b95cc1074
# ╠═ffd9315b-d5ad-4636-ba2f-5380c22055d2
# ╠═b5f4992a-5233-4ade-8ca4-bf2e6da1a4ac
# ╠═e36ebf1c-a750-40b4-8815-8bb3dee3910f
# ╟─7cba065d-98b4-4e04-a3b4-76224878a0bd
# ╠═c97e19e9-e818-432e-aaeb-618090372fad
# ╠═e54ae808-22f9-4ee5-9b75-d93aa56a8b3d
# ╠═74c3f2cb-6865-4cde-9c61-b75777854866
# ╟─488e83b4-4f5e-400f-85e2-f3b0716d7404
# ╠═5e48f2ae-bfc8-4130-a831-fae02b331bf7
# ╟─31f2d123-8a9a-45c6-9f80-569b070c54e4
# ╟─990f6bc0-f8d0-4374-b28c-82664e24984a
# ╟─e5d41d8f-3b62-434b-9491-98a017e90658
# ╟─40d92503-c7ab-4bf7-af62-02b08a284c17
# ╟─337d9aaf-e2f5-49cc-bd68-178d12110390
# ╠═2b93af3a-066b-443f-88c4-702b17e555b7
# ╟─e16a4171-e82d-4c43-af26-1f183fd15b2f
# ╠═e6ae2efe-d519-4c44-b706-a0ca95e6c050
# ╠═9de88c42-a73b-42f1-95be-e117b5beb77b
# ╟─11bd6e51-6993-4b8e-945c-7e61ff970dd4
# ╟─cc4b05de-6d6a-4c61-bdca-a0f6992ecfca
# ╠═2facd1bf-e9ac-46c1-9523-1244395a3e8b
# ╠═63fbadee-2ad4-4f68-9ba1-5c9d0143eb1a
# ╟─1d9033dd-e266-4192-930d-ef5be6780226
# ╠═dba75149-3e4c-4046-8fc7-42e8bf7375a5
# ╟─a16b0dc3-56c4-4f1c-94f2-ab034bfc30f3
# ╟─e3fb296e-c338-460f-b73f-fc5a4ffe2368
# ╟─c85762f2-6ad0-43aa-99a2-abd1e8249cdd
# ╠═d30212db-5306-4d3c-a786-f568ca67725f
# ╟─9be239b8-6c91-47c7-97c6-c39cddd3bc5d
# ╠═2c0adb00-71b6-483c-be10-f31642690184
# ╟─7557b534-c348-44bd-8de2-557a290e1b59
# ╠═5d78bd1f-47b3-422c-ae34-11f5e2ae93c3
# ╟─5514f1a7-bc64-4042-8c5e-b3dc835ebd2d
# ╟─0bf114a2-ec3a-4625-b60d-0287f447d2a4
# ╠═3d081d57-0bdb-448b-801f-b1a3bb8deb67
# ╟─5d715c51-147f-4a98-93cd-0c5fc2c51430
# ╠═3c9822b6-2643-4b0d-be1c-4716d701fbfa
# ╠═d253af17-92a7-4c62-a227-fb8cf9298939
# ╠═36de15f2-100d-4e21-8ca1-6a4b3a08a52d
# ╟─ed673d82-55eb-4ba8-9e64-410ecc64e499
# ╠═3d117fff-2447-4b3b-ab71-301b38e24d54
# ╟─8a786af1-9902-4f58-ba58-4cf9dd3b9d3e
# ╟─71ccf51e-2d8c-46d1-8b87-953bcbf8302f
# ╟─85ea82e9-35c6-41be-a575-ab823809eef7
# ╟─10403496-5a93-4890-b399-93d6086ff3ee
# ╟─2ba5f81f-019b-4749-a0a2-35ead50dee5f
# ╠═1cd04b3f-3a4e-4c93-85f0-4a41773c4366
# ╟─9bb894f5-6247-43b0-a632-0a034be4abe8
# ╠═ca2acf91-28d1-4ddc-8150-f4eecef3dffa
# ╟─2871b32a-180e-4d57-a1c8-44bc68ec1816
# ╠═bc1e3f75-8d92-4311-877a-6a4d7f227830
# ╟─b0eccdc2-40bc-44d7-8a72-3e059a4214ab
# ╠═81d7a663-1079-4613-afd8-de95eb73f8ae
# ╠═d93a8dca-9c63-4551-87e8-9939a765bef1
# ╠═8b184833-59d9-4958-8460-369b1ea19b9e
# ╠═0cac30f1-9101-4ba3-accc-52621bc1d16f
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using Test, Logging, Random, JuMP
using DecisionProgramming

function influence_diagram(rng::AbstractRNG, n_C::Int, n_D::Int, n_V::Int, m_C::Int, m_D::Int, states::Vector{Int}, n_inactive::Int)
    C, D, V = random_diagram(rng, n_C, n_D, n_V, m_C, m_D)
    S = States(rng, states, length(C) + length(D))
    X = [Probabilities(rng, c, S; n_inactive=n_inactive) for c in C]
    Y = [Consequences(rng, v, S; low=-1.0, high=1.0) for v in V]

    validate_influence_diagram(S, C, D, V)

    s_c = sortperm([c.j for c in C])
    s_d = sortperm([d.j for d in D])
    s_v = sortperm([v.j for v in V])
    C, D, V = C[s_c], D[s_d], V[s_v]
    X, Y = X[s_c], Y[s_v]
    P = DefaultPathProbability(C, X)
    U = DefaultPathUtility(V, Y)

    return D, S, P, U
end

function test_decision_model(D, S, P, U, n_inactive, probability_scale_factor, probability_cut)
    model = Model()

    @info "Testing DecisionVariables"
    z = DecisionVariables(model, S, D)

    @info "Testing PathCompatibilityVariables"
    x_s = PathCompatibilityVariables(model, z, S, P; probability_cut = probability_cut)

    @info "Testing PositivePathUtility"
    U′ = if probability_cut U else PositivePathUtility(S, U) end

    @info "Testing probability_cut"
    lazy_probability_cut(model, x_s, P)

    @info "Testing expected_value"
    if probability_scale_factor > 0
        EV = expected_value(model, x_s, U′, P; probability_scale_factor = probability_scale_factor)
    else
        @test_throws DomainError expected_value(model, x_s, U′, P; probability_scale_factor = probability_scale_factor)
    end

    @info "Testing conditional_value_at_risk"
    if probability_scale_factor > 0
        CVaR = conditional_value_at_risk(model, x_s, U′, P, 0.2; probability_scale_factor = probability_scale_factor)
    else
        @test_throws DomainError conditional_value_at_risk(model, x_s, U′, P, 0.2; probability_scale_factor = probability_scale_factor)
    end

    @test true
end

function test_analysis_and_printing(D, S, P, U)
    @info("Creating random decision strategy")
    Z_j = [LocalDecisionStrategy(rng, d, S) for d in D]
    Z = DecisionStrategy(D, Z_j)

    @info "Testing CompatiblePaths"
    @test all(true for s in CompatiblePaths(S, P.C, Z))
    @test_throws DomainError CompatiblePaths(S, P.C, Z, Dict(D[1].j => 1))
    node, state = (P.C[1].j, 1)
    @test all(s[node] == state for s in CompatiblePaths(S, P.C, Z, Dict(node => state)))

    @info "Testing UtilityDistribution"
    udist = UtilityDistribution(S, P, U, Z)

    @info "Testing StateProbabilities"
    sprobs = StateProbabilities(S, P, Z)

    @info "Testing conditional StateProbabilities"
    sprobs2 = StateProbabilities(S, P, Z, node, state, sprobs)

    @info "Testing "
    print_decision_strategy(S, Z)
    print_utility_distribution(udist)
    print_state_probabilities(sprobs, [c.j for c in P.C])
    print_state_probabilities(sprobs, [d.j for d in D])
    print_state_probabilities(sprobs2, [c.j for c in P.C])
    print_state_probabilities(sprobs2, [d.j for d in D])
    print_statistics(udist)
    print_risk_measures(udist, [0.0, 0.05, 0.1, 0.2, 1.0])

    @test true
end

@info "Testing model construction"
rng = MersenneTwister(4)
for (n_C, n_D, states, n_inactive, probability_scale_factor, probability_cut) in [
        (3, 2, [1, 2, 3], 0, 1.0, true),
        (3, 2, [1, 2, 3], 0, -1.0, true),
        (3, 2, [3], 1, 100.0, true),
        (3, 2, [1, 2, 3], 0, -1.0, false),
        (3, 2, [3], 1, 10.0, false)
    ]
    D, S, P, U = influence_diagram(rng, n_C, n_D, 2, 2, 2, states, n_inactive)
    test_decision_model(D, S, P, U, n_inactive, probability_scale_factor, probability_cut)
    test_analysis_and_printing(D, S, P, U)
end

ext: jl


text:
# Separable NMF

using NMF
using Printf

m, n, k = 7, 8, 3
W, H = NMF.separable_data(m, n, k)
println("Separable matrix X = ")
NMF.printf_mat(W*H)
println()

# Separable NMF
r = nnmf(W*H, k; alg=:spa, init=:spa)

# Results
println("Product of W and H = ")
NMF.printf_mat(r.W*r.H)
println()

println("Matrix W = ")
NMF.printf_mat(r.W)
println()

@printf("objvalue  = %.6e\n", r.objvalue)


text:
close("all")
transD_GP.getchi2forall(opt, fsize=6)
ax = gcf().axes;
ax[2].set_ylim(100, 200)
χ² = sum(2*csem.ndata)
ax[2].plot(xlim(), [χ²/2 , χ²/2], "--", color="gray")
savefig("csem_conv_s.png", dpi=300)
opt.xall[:] .= zall
transD_GP.plot_posterior(csem, opt, burninfrac=0.5, figsize=(4,4),
    cmappdf="inferno", qp1=0.1, qp2=0.9, nbins=100)
ax = gcf().axes
ax[1].step(log10.(ρ[2:end]), z[2:end], color="w", alpha=0.4)
ax[1].set_xlim(-1.0959789629712384, 2.5166552528237873)
savefig("csem_post_s.png", dpi=300)

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
@enum Severity::UInt8 Asymptomatic=1 Mild Severe Critical # 2bits

@enum HealthState::UInt8 Healthy Incubating Infectious MildSymptoms SevereSymptoms CriticalSymptoms Recovered Dead # 3 bits

@enum FreedomState::UInt8 Free HomeQuarantine HomeTreatment Hospitalized Released # 3 bits

@enum DetectionStatus::UInt8 Undetected UnderObservation TestPending Detected #2 bits

@enum ContactKind::UInt8 NoContact=0 HouseholdContact HospitalContact AgeCouplingContact ConstantKernelContact OutsideContact # 3 bits

@enum DetectionKind::UInt8 NoDetection=0 OutsideQuarantineDetection=1 FromQuarantineDetection FromTracingDetection

@enum TracingKind::UInt8 begin
  NotTraced=0
  MildCaseTraced #unused yet
  HospitalTraced #unused yet
  QuarantineTraced #unused yet
  ClassicalTraced
  PhoneTraced
end

@enum StrainKind::UInt8 NullStrain ChineseStrain BritishStrain DeltaStrain


text:
#==== Command line arguments ==================================================#

struct Constants
    K::Int
    R::Float64
    λ::Float64
    α::Float64
    β::Float64
    leps::Float64
    lleaps::Int
    heps::Float64
    hleaps::Int
    lλ::Float64  # used in prior
    αlβ::Float64 # calculations.
    Constants(K,R,λ,α,β,leps,heps,lleaps,hleaps) = new(
        K,R,λ,α,β,leps,heps,lleaps,hleaps,
        log(λ),
        α*log(β)-lgamma(α)
    )
end

#==== State ===================================================================#

mutable struct State
    k::Int
    locs::Vector{Float64}
    hts::Vector{Float64}
end

"""
Initialises a new state with evenly spaced step locations and mean heights.
"""
function State(Y::Vector{Float64}, k::Int, K::Int, R::Float64)
    S = State(k, fill(NaN,K+2), fill(NaN,K+1))
    S.locs[1:k+2] .= linspace(0.,R,k+2)
    s,y = 0.,1 # step location and data index.
    for i = 1:k+1
        sn = S.locs[i+1] # next step location.
        @views n = searchsortedfirst(Y[y:end], sn) - 1 # step data count.
        S.hts[i] = n/(sn-s)
        s,y = sn, y+n
    end
    return S
end

function header(K)
    io = IOBuffer()
    write(io, "k,")
    for k = 0:K+1
        write(io, "s$k,")
    end
    for k = 0:K
        write(io, "h$k")
        if (k < K) write(io, ",") end
    end
    return String(take!(io))
end

function string(S::State)
    k = S.k
    io = IOBuffer()
    write(io, string(k,","))
    for s in S.locs
        if !isnan(s) write(io, string(s)) end
        write(io, ",")
    end
    for (i,h) in enumerate(S.hts)
        if !isnan(h) write(io, string(h)) end
        if (i != length(S.hts)) write(io, ",") end
    end
    return String(take!(io))
end

#==== Extended state ==========================================================#

struct Auxiliary
    locs::Vector{Float64}
    hts::Vector{Float64}
    widths::Vector{Float64}
    counts::Vector{Int}
    tlocs::Vector{Float64} # t = logit(s/R).
    lhts::Vector{Float64}  # j = ln(h).
    ∇tlocs::Vector{Float64}
    ∇lhts::Vector{Float64}
    ptlocs::Vector{Float64}
    plhts::Vector{Float64}
    ix::Vector{Int} # index vector used for sorting.
    Auxiliary(K) = new(
        Vector{Float64}(K+2),
        Vector{Float64}(K+1),
        Vector{Float64}(K+1),
        Vector{Int}(K+1),
        Vector{Float64}(K),
        Vector{Float64}(K+1),
        Vector{Float64}(K),
        Vector{Float64}(K+1),
        Vector{Float64}(K),
        Vector{Float64}(K+1),
        collect(1:K+1)
    )
end

#==== Birth/death probabilities ===============================================#

struct Jumps
    bprob::Function
    dprob::Function
end

"""
Initialises a new pair of jump proposals such that the sum of their returned
probabilities never exceeds the given value.
"""
function Jumps(b::Function, d::Function, kmin::Int, kmax::Int, s::Float64=.9)
    c = max(b(kmin),d(kmax))
    if kmin < kmax-1
        c = max(c, maximum(b(k)+d(k) for k = kmin+1:kmax-1))
    end
    c = s/c
    Jumps(k -> (k < kmax) ? c*b(k) : 0,
          k -> (k > kmin) ? c*d(k) : 0)
end

#==== Move statistics =========================================================#

mutable struct Counts
    bprops::Int
    dprops::Int
    lprops::Int
    hprops::Int
    baccepts::Int
    daccepts::Int
    laccepts::Int
    haccepts::Int
    lsorts::Int
    ks::Vector{Int}
    Counts(K) = new(0,0,0,0,0,0,0,0,0,zeros(Int,K+1))
end

function show(io::IO, C::Counts)

    K = length(C.ks) - 1
    println(io, "No. of change points:")
    for k = 0:K
        m = C.ks[k+1]
        (m > 0) && println(io, "  $k: $m")
    end
    println(io, "Mean: $(dot(C.ks,0:K) / sum(C.ks))")

    println(io, "Acceptance rates:")
    @printf(io, "  birth: %d/%d (%.2f%%)\n",
        C.baccepts, C.bprops, (C.bprops == 0) ? 0 : 100*C.baccepts/C.bprops)
    @printf(io, "  death: %d/%d (%.2f%%)\n",
        C.daccepts, C.dprops, (C.dprops == 0) ? 0 : 100*C.daccepts/C.dprops)
    @printf(io, "  location: %d/%d (%.2f%%)\n",
        C.laccepts, C.lprops, (C.lprops == 0) ? 0 : 100*C.laccepts/C.lprops)
    @printf(io, "  height: %d/%d (%.2f%%)",
        C.haccepts, C.hprops, (C.hprops == 0) ? 0 : 100*C.haccepts/C.hprops)
    if C.lsorts > 0
        print(io, "\nLocation sorts: $(C.lsorts)")
    end

end

max_stars_repo_name: kevdur/rjcp

max_issues_repo_name: kevdur/rjcp

max_forks_repo_name: kevdur/rjcp


text:
"""
Creates Ohms constraints (yt post fix indicates that Y and T values are in rectangular form)
```
p[f_idx] ==  (g+g_fr)/tm*v[f_bus]^2 + (-g*tr+b*ti)/tm^2*(v[f_bus]*v[t_bus]*cos(t[f_bus]-t[t_bus])) + (-b*tr-g*ti)/tm^2*(v[f_bus]*v[t_bus]*sin(t[f_bus]-t[t_bus]))
q[f_idx] == -(b+b_fr)/tm*v[f_bus]^2 - (-b*tr-g*ti)/tm^2*(v[f_bus]*v[t_bus]*cos(t[f_bus]-t[t_bus])) + (-g*tr+b*ti)/tm^2*(v[f_bus]*v[t_bus]*sin(t[f_bus]-t[t_bus]))
```
"""
function constraint_ohms_from_variable_impedance(pm::_PM.AbstractACPModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tr, ti, tm)
    (l,i,j) = f_idx
    p_fr  = _PM.var(pm, n,  :p, f_idx)
    q_fr  = _PM.var(pm, n,  :q, f_idx)
    vm_fr = _PM.var(pm, n, :vm, f_bus)
    vm_to = _PM.var(pm, n, :vm, t_bus)
    va_fr = _PM.var(pm, n, :va, f_bus)
    va_to = _PM.var(pm, n, :va, t_bus)

    b = _PM.var(pm, n, :b, l)
    g = _PM.var(pm, n, :g, l)    # g = y*r
    b_fr = _PM.var(pm, n, :b_shunt, l)
    g_fr = _PM.var(pm, n, :g_shunt, l)

    JuMP.@NLconstraint(pm.model, p_fr ==  (g+g_fr)/tm^2*vm_fr^2
    + (-g*tr+b*ti)/tm^2*(vm_fr*vm_to*cos(va_fr-va_to))
    + (-b*tr-g*ti)/tm^2*(vm_fr*vm_to*sin(va_fr-va_to)) )
    JuMP.@NLconstraint(pm.model, q_fr == -(b+b_fr)/tm^2*vm_fr^2
    - (-b*tr-g*ti)/tm^2*(vm_fr*vm_to*cos(va_fr-va_to))
    + (-g*tr+b*ti)/tm^2*(vm_fr*vm_to*sin(va_fr-va_to)) )
end

"""
Creates Ohms constraints (yt post fix indicates that Y and T values are in rectangular form)
```
p[t_idx] ==  (g+g_to)*v[t_bus]^2 + (-g*tr-b*ti)/tm^2*(v[t_bus]*v[f_bus]*cos(t[t_bus]-t[f_bus])) + (-b*tr+g*ti)/tm^2*(v[t_bus]*v[f_bus]*sin(t[t_bus]-t[f_bus]))
q[t_idx] == -(b+b_to)*v[t_bus]^2 - (-b*tr+g*ti)/tm^2*(v[t_bus]*v[f_bus]*cos(t[f_bus]-t[t_bus])) + (-g*tr-b*ti)/tm^2*(v[t_bus]*v[f_bus]*sin(t[t_bus]-t[f_bus]))
```
"""
function constraint_ohms_to_variable_impedance(pm::_PM.AbstractACPModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_to, b_to, tr, ti, tm)
    (l,i,j) = f_idx
    p_to  = _PM.var(pm, n,  :p, t_idx)
    q_to  = _PM.var(pm, n,  :q, t_idx)
    vm_fr = _PM.var(pm, n, :vm, f_bus)
    vm_to = _PM.var(pm, n, :vm, t_bus)
    va_fr = _PM.var(pm, n, :va, f_bus)
    va_to = _PM.var(pm, n, :va, t_bus)

    # r = g/b
    b = _PM.var(pm, n, :b, l)
    g = _PM.var(pm, n, :g, l)
    b_to = _PM.var(pm, n, :b_shunt, l)
    g_to = _PM.var(pm, n, :g_shunt, l)
    # g = y*r

    JuMP.@NLconstraint(pm.model, p_to ==  (g+g_to)*vm_to^2
    + (-g*tr-b*ti)/tm^2*(vm_to*vm_fr*cos(va_to-va_fr))
    + (-b*tr+g*ti)/tm^2*(vm_to*vm_fr*sin(va_to-va_fr)) )
    JuMP.@NLconstraint(pm.model, q_to == -(b+b_to)*vm_to^2
    - (-b*tr+g*ti)/tm^2*(vm_to*vm_fr*cos(va_to-va_fr))
    + (-g*tr-b*ti)/tm^2*(vm_to*vm_fr*sin(va_to-va_fr)) )
end

"""
Defines the faithfullness in terms of grid losses
|actual_loss - reference_loss|/(reference_loss) <= beta
we split this in two constraints to deal with the absolute value
|x| <= y iff x <=y and -x<=y
"""
function constraint_loss_faithfulness(pm::_PM.AbstractPowerModel, n::Int, ref_loss, beta)
    arcs_from = _PM.ref(pm, :arcs_from)
    p = _PM.var(pm, :p)

    loss = sum(p[(l,i,j)] + p[(l,j,i)] for (l,i,j) in arcs_from)

    JuMP.@constraint(pm.model, (loss - ref_loss)/(ref_loss) <= beta)
    JuMP.@constraint(pm.model, (ref_loss - loss)/(ref_loss) <= beta)
end

function constraint_cost_faithfulness(pm::_PM.AbstractPowerModel, n::Int, ref_cost, beta)
    cost = _PM.var(pm, :cost)

    JuMP.@constraint(pm.model, (cost - ref_cost)/(ref_cost) <= beta)
    JuMP.@constraint(pm.model, (ref_cost - cost)/(ref_cost) <= beta)
end

"Supporting only quadratic cost functions for generators breaks compatibility with Matpower to some extent"
function constraint_fuel_cost_quadratic(pm::_PM.AbstractPowerModel)
    gen_cost = Dict()
    for (n, nw_ref) in _PM.nws(pm)
        for (i,gen) in nw_ref[:gen]
            pg = sum( _PM.var(pm, n, :pg, i)[c] for c in _PM.conductor_ids(pm, n) )

            if length(gen["cost"]) == 1
                gen_cost[(n,i)] = gen["cost"][1]
            elseif length(gen["cost"]) == 2
                gen_cost[(n,i)] = gen["cost"][1]*pg + gen["cost"][2]
            elseif length(gen["cost"]) == 3
                gen_cost[(n,i)] = gen["cost"][1]*pg^2 + gen["cost"][2]*pg + gen["cost"][3]
            else
                gen_cost[(n,i)] = 0.0
            end
        end
    end

    cost = _PM.var(pm, :cost)
    JuMP.@constraint(pm.model, cost ==
        sum(
            sum( gen_cost[(n,i)] for (i,gen) in nw_ref[:gen] )
        for (n, nw_ref) in _PM.nws(pm))
            )
end

function constraint_gen_bounds_cc(pm::_PM.AbstractPowerModel, n::Int, i, pmin, pmax, qmin, qmax, eta)
    pg = _PM.var(pm, n, :pg, i)
    qg = _PM.var(pm, n, :qg, i)
    d = Distributions.Normal()
    z = Distributions.quantile(d,eta)

    # JuMP.@constraint(pm.model, (cost - ref_cost)/(ref_cost) <= beta)
end


text:
import NiLang: i_dot, i_norm2, i_axpy!
abstract type Point{T,N} end
export P2, P3

struct P2{T} <: Point{T,2}
    x::T
    y::T
end

struct P3{T} <: Point{T,3}
    x::T
    y::T
    z::T
end

Base.zero(::Type{P2{T}}) where T = P2(zero(T), zero(T))
Base.zero(::Type{P3{T}}) where T = P3(zero(T), zero(T), zero(T))
Base.zero(::P2{T}) where T = P2(zero(T), zero(T))
Base.zero(::P3{T}) where T = P3(zero(T), zero(T), zero(T))
NiLang.AD.GVar(x::P2) = P2(GVar(x.x), GVar(x.y))
(_::Type{Inv{GVar}})(x::P2) = P2((~GVar)(x.x), (~GVar)(x.y))
NiLang.AD.GVar(x::P2, g::P2) = P2(GVar(x.x, g.x), GVar(x.y, g.y))

NiLang.AD.GVar(x::P3) = P3(GVar(x.x), GVar(x.y), GVar(x.z))
(_::Type{Inv{GVar}})(x::P3) = P3((~GVar)(x.x), (~GVar)(x.y), (~GVar)(x.z))
NiLang.AD.GVar(x::P3, g::P3) = P3(GVar(x.x, g.x), GVar(x.y, g.y), GVar(x.z, g.z))

StaticArrays.SVector(p::P3) = SVector(p.x, p.y, p.z)
StaticArrays.SVector(p::P2) = SVector(p.x, p.y)

function vec2cam(v::AbstractVector)
    Camera(P3(v[1:3]...), P3(v[4:6]...), v[7], P2(v[8:9]...), P2(v[10:11]...))
end

rand_camera() = vec2cam(randn(11))
#mzeros(n) = MVector{n}(zeros(n))
NiLang.AD.GVar(cam::Camera) = Camera(GVar(cam.r), GVar(cam.c), GVar(cam.f), GVar(cam.x0), GVar(cam.κ))

##################### objective #############################

@i @inline function i_dot(out!, x::P2, y::P2)
    out! += x.x * y.x
    out! += x.y * y.y
end

@i @inline function i_dot(out!, x::P3, y::P3)
    out! += x.x * y.x
    out! += x.y * y.y
    out! += x.z * y.z
end

@i @inline function i_norm2(out!, x::P2)
    out! += abs2(x.x)
    out! += abs2(x.y)
end

@i @inline function i_norm2(out!, x::P3)
    out! += abs2(x.x)
    out! += abs2(x.y)
    out! += abs2(x.z)
end

@i @inline function i_cross!(out!, x, y)
    @inbounds begin
        out![1] += x[2] * y[3]
        out![1] -= x[3] * y[2]
        out![2] += x[3] * y[1]
        out![2] -= x[1] * y[3]
        out![3] += x[1] * y[2]
        out![3] -= x[2] * y[1]
    end
end

@i @inline function i_cross!(out!, x::P3, y::P3)
    @inbounds begin
        out!.x += x.y * y.z
        out!.x -= x.z * y.y
        out!.y += x.z * y.x
        out!.y -= x.x * y.z
        out!.z += x.x * y.y
        out!.z -= x.y * y.x
    end
end

@i @inline function i_axpy!(a, X::P2, Y::P2)
    Y.x += a * X.x
    Y.y += a * X.y
end

@i @inline function i_axpy!(a, X::P3, Y::P3)
    Y.x += a * X.x
    Y.y += a * X.y
    Y.z += a * X.z
end

@i @inline function rodrigues_rotate_point(out!::P3{T}, rot::P3{T}, X::P3{T}) where T
    @invcheckoff sqtheta ← zero(T)
    i_norm2(sqtheta, rot)
    @invcheckoff if (sqtheta > 1e-10, ~)
        @routine begin
            @zeros T tmp1 tmp2 wX s c theta invtheta
            @zeros P3{T} w w_cross_X
            theta += sqtheta ^ 0.5
            (s, c) += sincos(theta)
            tmp1 += 1 - c
            invtheta += 1 / theta
            i_axpy!(invtheta, rot, w)
            i_cross!(w_cross_X, w, X)
            i_dot(wX, w, X)
            tmp2 += wX * tmp1
        end
        i_axpy!(c, X, out!)
        i_axpy!(s, w_cross_X, out!)
        i_axpy!(tmp2, w, out!)
        ~@routine
    else
        @routine @zeros P3{T} w_cross_X
        i_cross!(w_cross_X, rot, X)
        out! += X + w_cross_X
        ~@routine
    end
    (~i_norm2)(sqtheta, rot)
    @invcheckoff sqtheta → zero(T)
end

@i @inline function radial_distort(out!::P2{T}, rad_params::P2{T}, proj::P2{T}) where T
    @routine @invcheckoff begin
        @zeros T rsq L rsqsq
        i_norm2(rsq, proj)
        rsqsq += rsq ^ 2
        @inbounds L += rad_params.y * rsqsq
        @inbounds L += rad_params.x * rsq
        L += 1
    end
    i_axpy!(L, proj, out!)
    ~@routine
end

@i @inline function :(+=)(identity)(x::P2, y::P2)
    x.x += y.x
    x.y += y.y
end

@i @inline function :(+=)(+)(x::P2, y::P2, z::P2)
    x.x += y.x + z.x
    x.y += y.y + z.y
end

@i @inline function :(+=)(-)(x::P2, y::P2, z::P2)
    x.x += y.x - z.x
    x.y += y.y - z.y
end

@i @inline function :(+=)(identity)(x::P3, y::P3)
    x.x += y.x
    x.y += y.y
    x.z += y.z
end

@i @inline function :(+=)(+)(x::P3, y::P3, z::P3)
    x.x += y.x + z.x
    x.y += y.y + z.y
    x.z += y.z + z.z
end

@i @inline function :(+=)(-)(x::P3, y::P3, z::P3)
    x.x += y.x - z.x
    x.y += y.y - z.y
    x.z += y.z - z.z
end

@i function project(out!, cam::Camera{T}, X) where T
    @routine @invcheckoff begin
        Xcam ← zero(cam.r)
        Xcam2 ← zero(cam.x0)
        distorted ← zero(cam.x0)
        X -= cam.c
        rodrigues_rotate_point(Xcam, cam.r, X)
        Xcam2.x += Xcam.x / Xcam.z
        Xcam2.y += Xcam.y / Xcam.z
        radial_distort(distorted, cam.κ, Xcam2)
    end
    i_axpy!(cam.f, distorted, out!)
    out! += cam.x0
    ~@routine
end

@i function compute_reproj_err(out!, diff, cam, X, w, feat)
    project(diff, cam, X)
    diff -= feat
    i_axpy!(w, diff, out!)
end

@i @inline function compute_w_err(out!, w)
    out! += 1.0
    out! -= w ^ 2
end

function compute_ba_J(::Val{:NiLang}, cams::AbstractArray{<:Camera{T}}, X::AbstractArray{<:P3}, w, obs, feats::AbstractArray{<:P2}) where T
    p = size(obs,2)
    reproj_err_d = zeros(2*p, 15)
    for i in 1:p
        idx = (2*(i-1))+1
        @inbounds j = obs[1,i]
        @inbounds l = obs[2,i]
        cam = cams[j]
        out! = zero(cam.x0)
        diff = zero(cam.x0)
        res = compute_reproj_err(out!, diff, cam, X[l], w[i], feats[i])
        _, _, gcam, gX, gw, _ = (~compute_reproj_err)(P2(GVar(res[1].x, 1.0), GVar(res[1].y, 0.0)),
                        GVar(res[2]), GVar(res[3]), GVar(res[4]), GVar(res[5]), GVar(res[6]))
        packgrad!(reproj_err_d, idx, gcam, gX, gw)
        _, _, gcam, gX, gw, _ = (~compute_reproj_err)(P2(GVar(res[1].x, 0.0), GVar(res[1].y, 1.0)),
                        GVar(res[2]), GVar(res[3]), GVar(res[4]), GVar(res[5]), GVar(res[6]))
        packgrad!(reproj_err_d, idx+1, gcam, gX, gw)
    end
    w_err_d = zeros(1,p)
    for i in 1:p
        @inbounds w_err_d[i] = grad(Grad(compute_w_err)(Val(1), 0.0, w[i])[3])
    end
    (reproj_err_d, w_err_d)
end

@inline function packgrad!(out, idx, gcam, gX, gw)
    @inbounds begin
        out[idx,1] = grad(gcam.r.x)
        out[idx,2] = grad(gcam.r.y)
        out[idx,3] = grad(gcam.r.z)
        out[idx,4] = grad(gcam.c.x)
        out[idx,5] = grad(gcam.c.y)
        out[idx,6] = grad(gcam.c.z)
        out[idx,7] = grad(gcam.f)
        out[idx,8] = grad(gcam.x0.x)
        out[idx,9] = grad(gcam.x0.y)
        out[idx,10] = grad(gcam.κ.x)
        out[idx,11] = grad(gcam.κ.y)
        out[idx,12] = grad(gX.x)
        out[idx,13] = grad(gX.y)
        out[idx,14] = grad(gX.z)
        out[idx,15] = grad(gw)
    end
end

@i function ba_objective!(reproj_err!, w_err!, reproj_err_cache!, cams, X, w, obs, feats)
    for i in 1:length(feats)
        compute_reproj_err(reproj_err![i], reproj_err_cache![i], cams[obs[1,i]], X[obs[2,i]], w[i], feats[i])
    end
    @inbounds for i=1:length(w_err!)
        w_err![i] += 1.0
        w_err![i] -= w[i] ^ 2
    end
end


text:
include("testsystem1.jl")
srand(12)
Y, X = sample(20, 100, M)

function testmean(X, μ)
    @test norm(mean(X) - μ) < 1.96*std(X)/sqrt(length(X))
end
# Test the test: [try testmean(randn(1000), 0.0).value catch false end for i in 1:1000]

testmean(X[1, :, 1], 1.0)

function tostatic(M::LinearHomogSystem)
    d2, d = Kalman.dims(M)
    LinearHomogSystem(SVector{d}(M.x0),
    SMatrix{d,d}(M.P0),
    SMatrix{d,d}(M.Phi),
    SVector{d}(M.b),
    SMatrix{d,d}(M.Q),
    SVector{d2}(M.y),
    SMatrix{d2,d}(M.H),
    SMatrix{d2,d2}(M.R))
end
M2 = tostatic(M)

include("testsystem.jl")
Y, X = sample(20, 100, M)

# Reproducibility
@test sample(MersenneTwister(1), 5, 5, M) == sample(MersenneTwister(1), 5, 5, M)


text:
module ToroidalHarmonics

include("util.jl")
include("TORH1.jl")
export TORH1

end # module


text:
#   This file is part of Fatou.jl. It is licensed under the MIT license
#   Copyright (C) 2017 Michael Reed

export orbit

"""
    orbit(K::Fatou.Define)
    orbit(u::Function, ∂, orbit, depth, n)

Plot funciton compositions of the primary `Fatou.Define` function up to any `depth` including cobweb plot with an `orbit` depth from the `x0` start point.

# Examples
```Julia
julia> juliafill("z^2-0.67",∂=[-1.25,1.5],x0=1.25,orbit=17,depth=3) |> orbit
```
"""
function orbit(K::Define{FT,QT,CT}) where {FT,QT,CT}
    !isdefined(Fatou,:UnicodePlots) && !isdefined(Fatou,:PyPlot) && !isdefined(Fatou,:Makie) && throw(error("Requires `using PyPlot` or `using UnicodePlots` or `using Makie`"))
    K.x0 == nothing ? (bi = K.Ω.∂[1:2]') : (bi = [K.Ω.∂[1:2]...,K.x0]')
    orbit(K.E,z->K.F(z,0),convert(Array{Float64},bi),K.orbit,K.depth,Int(K.Ω.n))
end

function real_orb(E,f::Function,bi::Matrix{Float64},orb::Int=0,depth::Int=1,incr::Int=384)
    # initalize array to depth
    N = zeros(incr,depth+1)
    bis = zeros(3)
    bis[1:length(bi)] = bi[:]
    # set x-axis coordinate set
    x = range(bi[1],stop=bi[2],length=incr)
    N[:,1] = x[:]
    # loop over all discrete x-axis points in set
    # loop function composition at x val
    for t ∈ 1:depth
        N[:,t+1] = broadcast(f,N[:,t])
    end
    # initialize for orbit cobweb
    N2 = zeros(orb+1)
    N2[1] = bis[3];
    # loop over orbit cobweb
    for t ∈ 1:orb
        # evaluate function composition
        N2[t+1] = f(N2[t])
    end
    # interleave cobweb orbit data
    siz = 3*(length(N2)-1)
    orbit = zeros(siz,2)
    orbit[1:3:siz,1] = N2[1:end-1]
    orbit[2:3:siz,1] = N2[1:end-1]
    orbit[3:3:siz,1] = N2[2:end]
    orbit[1:3:siz,2] = N2[1:end-1]
    orbit[2:3:siz,2] = N2[2:end]
    orbit[3:3:siz,2] = N2[2:end]
    return x,N,N2,orbit,bis
end


text:
xvals = 0.01:0.01:8.0
yvals = map(xvals) do x
    comp1 = pdf(LogNormal(), x)
    comp2 = pdf(Uniform(2.0, 3.0), x)
    0.5 * comp1 + 0.5 * comp2
end

p = Plots.plot()
Plots.plot!(p, xvals, yvals, labels = "Real distribution")

kde = KernelDensity.kde(xs, bandwidth = 0.1)
Plots.plot!(p, kde.x, kde.density, labels = "KDE")

Plots.plot!(p, xvals, yvals, labels = "Real distribution")
Plots.xlims!(p, 0.0, 8.0)
Plots.savefig(joinpath(base_img, "compare_kde.pdf"))


text:
export PositiveDefR

struct PositiveDefR{VT,MT,C}
    Wr::MT
    Wi::MT
    b::VT
    σ::C
end
functor(d::PositiveDefR) = (Wr=d.Wr, Wi=d.Wi, b=d.b), t -> PositiveDefR(t.Wr, t.Wi, t.b, d.σ)

# NOTE This is based upon the fact that σr/c are real so that conj(W*σc) == conj(W)*σc
(l::PositiveDefR)((σr, σc)) = l.σ.(       l.Wr * (σr .+ σc) .+
                                    im .* l.Wi * (σr .- σc) .+ l.b)

PositiveDefR(in::Integer, args...;kwargs...) =
    PositiveDefR(STD_REAL_PREC, in, args...;kwargs...)
function PositiveDefR(T::Type{<:Real}, in::Integer, out::Integer, σ = identity;
               initW = glorot_uniform, initb = glorot_uniform)
  return PositiveDefR(initW(T, out, in), initW(T, out, in),
                      initb(T, out), σ)
end

function Base.show(io::IO, l::PositiveDefR)
  print(io, "PositiveDefR{$(eltype(l.Wr))}(", size(l.Wr, 2), ", ", size(l.Wr, 1))
  l.σ == identity || print(io, ", ", l.σ)
  print(io, ")")
end

struct PositiveDefRCache{Ta,Tb,Tc,Td,Te}
    σr::Tc
    σc::Tc
    ∑σ::Tc
    Δσ::Tc

    out::Tb

    δℒℒ::Td
    δℒr::Te
    δℒc::Te

    θ::Tb
    θr::Ta
    θr2::Ta
    valid::Bool
end

cache(l::PositiveDefR{Ta,Tb}, arr_T, in_T, in_sz) where {Ta,Tb} =
    PositiveDefRCache(
               similar(l.Wr, size(l.Wr,2)),
               similar(l.Wr, size(l.Wr,2)),
               similar(l.Wr, size(l.Wr,2)),
               similar(l.Wr, size(l.Wr,2)),

               similar(l.b, Complex{eltype(l.b)}),

               similar(l.Wr, Complex{eltype(l.b)}, size(l.Wr,1)),
               similar(l.Wr, Complex{eltype(l.b)}, 1, size(l.Wr, 2)),
               similar(l.Wr, Complex{eltype(l.b)}, 1, size(l.Wr, 2)),

               similar(l.b, Complex{eltype(l.b)}),
               similar(l.b),
               similar(l.b),
               false)

function layer_out_type_size(l::PositiveDefR, in_T ,in_sz)
    T1     = promote_type(in_T, Complex{eltype(l.Wr)})
    out_T  = promote_type(T1, Complex{eltype(l.b)})
    out_sz = size(l.b)
    return out_T, out_sz
end

function (l::PositiveDefR)(c::PositiveDefRCache, (xr, xc))
    # The preallocated caches
    out  = c.out
    θ    = c.θ
    θᵣ₁  = c.θr
    θᵣ₂  = c.θr2

    # Store the input to this layer for the backpropagation
    σr = copyto!(c.σr, xr)
    σc = copyto!(c.σc, xc)

    c.∑σ .= σr .+ σc
    c.Δσ .= σr .- σc

    #θ .= net.b .+ net.W * x
    mul!(θᵣ₁, l.Wr, c.∑σ)
    mul!(θᵣ₂, l.Wi, c.Δσ)
    # NOTE This is based upon the fact that σr/c are real so that conj(W*σc) == conj(W)*σc
    θ .= θᵣ₁ .+ im .* θᵣ₂ .+ l.b

    # Apply the nonlinear function
    out  .= l.σ.(θ)
    return out
end

function backprop(∇, l::PositiveDefR, c::PositiveDefRCache, δℒ)
    # The preallocated caches
    θ = c.θ
    δℒℒ = c.δℒℒ

    # Compute the actual sensitivity
    copyto!(δℒℒ, δℒ)
    δℒℒ .*= fwd_der.(l.σ, θ)

    ∇.Wr .= δℒℒ .*       transpose(c.∑σ)
    ∇.Wi .= δℒℒ .* im .* transpose(c.Δσ)
    ∇.b .= δℒℒ

    # TODO alloc of transpose
    δℒr = 0 #mul!(c.δℒr, transpose(δℒℒ), l.Wr)
    δℒc = 0 #mul!(c.δℒc, transpose(δℒℒ), l.Wc)
    #c.δℒc.*= im

    return (δℒr, δℒc)
end


text:
# Inception modules
struct InceptionA
    c1_alone::ConvBN
    c1_before_3::ConvBN
    c1_before_d3::ConvBN
    c1_after_pool::ConvBN

    c3::ConvBN
    cd3_1::ConvBN
    cd3_2::ConvBN
    pool_mode::Int
end

function InceptionA(cx, num_1, num_1_before_3, num_3, num_1_before_d3, num_d3, num_1_after_pool, pool_mode)
    c1_alone = ConvBN(1, 1, cx, num_1)
    c1_before_3 = ConvBN(1, 1, cx, num_1_before_3)
    c1_before_d3 = ConvBN(1, 1, cx, num_1_before_d3)
    c1_after_pool = ConvBN(1, 1, cx, num_1_after_pool)

    c3 = ConvBN(3, 3, num_1_before_3, num_3, padding=1)
    cd3_1 = ConvBN(3, 3, num_1_before_d3, num_d3, padding=1)
    cd3_2 = ConvBN(3, 3, num_d3, num_d3, padding=1)
    return InceptionA(c1_alone, c1_before_3, c1_before_d3, c1_after_pool, c3, cd3_1, cd3_2, pool_mode)
end

function (i::InceptionA)(x)
    y1 = i.c1_alone(x)
    y2 = i.c3(i.c1_before_3(x))
    y3 = i.cd3_2(i.cd3_1(i.c1_before_d3(x)))
    y4 = i.c1_after_pool(pool(x, window=3, stride=1, padding=1, mode=i.pool_mode))
    old_size = size(y1)

    y1 = reshape(y1, :, size(y1)[end])
    y2 = reshape(y2, :, size(y2)[end])
    y3 = reshape(y3, :, size(y3)[end])
    y4 = reshape(y4, :, size(y4)[end])

    y_2d = vcat(y1, y2, y3, y4)
    return reshape(y_2d, old_size[1], old_size[2], :, old_size[4])
end

struct InceptionB
    c1_before_3::ConvBN
    c1_before_d3::ConvBN

    c3::ConvBN
    cd3_1::ConvBN
    cd3_2::ConvBN
end
function InceptionB(cx, num_1_before_3, num_3, num_1_before_d3, num_d3)
    c1_before_3 = ConvBN(1, 1, cx, num_1_before_3)
    c1_before_d3 = ConvBN(1, 1, cx, num_1_before_d3)

    c3 = ConvBN(3, 3, num_1_before_3, num_3, padding=1, stride=2)
    cd3_1 = ConvBN(3, 3, num_1_before_d3, num_d3, padding=1)
    cd3_2 = ConvBN(3, 3, num_d3, num_d3, padding=1, stride=2)
    return InceptionB(c1_before_3, c1_before_d3, c3, cd3_1, cd3_2)
end

function (i::InceptionB)(x)
    y1 = i.c3(i.c1_before_3(x))
    y2 = i.cd3_2(i.cd3_1(i.c1_before_d3(x)))
    y3 = pool(x, window=3, stride=2, padding=1)
    old_size = size(y1)

    y1 = reshape(y1, :, size(y1)[end])
    y2 = reshape(y2, :, size(y2)[end])
    y3 = reshape(y3, :, size(y3)[end])

    y_2d = vcat(y1, y2, y3)
    return reshape(y_2d, old_size[1], old_size[2], :, old_size[4])
end

#------------------------ Net2DeeperNet versions -------------------------------

# Deeper version of InceptionA
struct InceptionADeeper
    c1_before_3::ConvBN
    c1_before_d3::ConvBN
    c1_after_pool::ConvBN
    c1_alone::ConvBN

    c3s::Array{ConvBN}
    cd3s::Array{ConvBN}

    pool_mode::Int
end

# This constructor is used during random initialization (baseline for Net2DeeperNet)
function InceptionADeeper(cx, num_1, num_1_before_3, num_3, num_1_before_d3, num_d3, num_1_after_pool, pool_mode, deepening_factor::Int=2)
    c1_before_3 = ConvBN(1, 1, cx, num_1_before_3)
    c1_before_d3 = ConvBN(1, 1, cx, num_1_before_d3)
    c1_after_pool = ConvBN(1, 1, cx, num_1_after_pool)
    c1_alone = ConvBN(1, 1, cx, num_1)

    c3s = [ConvBN(3, 3, num_1_before_3, num_3, padding=1)]
    cd3s = [ConvBN(3, 3, num_1_before_d3, num_d3, padding=1), ConvBN(3, 3, num_d3, num_d3, padding=1)]

    for i in 1:deepening_factor
        push!(c3s, ConvBN(3, 3, num_3, num_3, padding=1))
        push!(cd3s, ConvBN(3, 3, num_d3, num_d3, padding=1), ConvBN(3, 3, num_d3, num_d3, padding=1))
    end

    return InceptionADeeper(c1_before_3, c1_before_d3, c1_after_pool, c1_alone, c3s, cd3s, pool_mode)
end

function (i::InceptionADeeper)(x)
    y1 = i.c1_alone(x)
    y2 = i.c1_before_3(x); for l in i.c3s; y2 = l(y2); end;
    y3 = i.c1_before_d3(x); for l in i.cd3s; y3 = l(y3); end;
    y4 = i.c1_after_pool(pool(x, window=3, stride=1, padding=1, mode=i.pool_mode))
    old_size = size(y1)

    y1 = reshape(y1, :, size(y1)[end])
    y2 = reshape(y2, :, size(y2)[end])
    y3 = reshape(y3, :, size(y3)[end])
    y4 = reshape(y4, :, size(y4)[end])

    y_2d = vcat(y1, y2, y3, y4)
    return reshape(y_2d, old_size[1], old_size[2], :, old_size[4])
end

# Deeper version of InceptionB
struct InceptionBDeeper
    c1_before_3::ConvBN
    c1_before_d3::ConvBN

    c3s::Array{ConvBN}
    cd3s::Array{ConvBN}
end

# This constructor is used during random initialization (baseline for Net2DeeperNet)
function InceptionBDeeper(cx, num_1_before_3, num_3, num_1_before_d3, num_d3, deepening_factor::Int=2)
    c1_before_3 = ConvBN(1, 1, cx, num_1_before_3)
    c1_before_d3 = ConvBN(1, 1, cx, num_1_before_d3)

    c3s = [ConvBN(3, 3, num_1_before_3, num_3, padding=1, stride=2)]
    cd3s = [ConvBN(3, 3, num_1_before_d3, num_d3, padding=1), ConvBN(3, 3, num_d3, num_d3, padding=1, stride=2)]

    for i in 1:deepening_factor
        push!(c3s, ConvBN(3, 3, num_3, num_3, padding=1))
        push!(cd3s, ConvBN(3, 3, num_d3, num_d3, padding=1), ConvBN(3, 3, num_d3, num_d3, padding=1))
    end
    return InceptionBDeeper(c1_before_3, c1_before_d3, c3s, cd3s)
end

function (i::InceptionBDeeper)(x)
    y1 = i.c1_before_3(x); for l in i.c3s; y1 = l(y1); end;
    y2 = i.c1_before_d3(x); for l in i.cd3s; y2 = l(y2); end;
    y3 = pool(x, window=3, stride=2, padding=1)
    old_size = size(y1)

    y1 = reshape(y1, :, size(y1)[end])
    y2 = reshape(y2, :, size(y2)[end])
    y3 = reshape(y3, :, size(y3)[end])

    y_2d = vcat(y1, y2, y3)
    return reshape(y_2d, old_size[1], old_size[2], :, old_size[4])
end

#------------------------------ Network creator functions ----------------------

"Builds an original Inception-BN network model"
function create_inception_bn_model(num_channels::Int, num_classes::Int)
    Chain(
        ConvBN(7, 7, num_channels, 64, padding=3, stride=2),
        Pool(3, 2, 0, 0),

        ConvBN(1, 1, 64, 64),
        ConvBN(3, 3, 64, 192, padding=1),
        Pool(3, 2, 0, 0),

        InceptionA(192, 64, 64, 64, 64, 96, 32, 2),
        InceptionA(256, 64, 64, 96, 64, 96, 64, 2),
        InceptionB(320, 128, 160, 64, 96),

        InceptionA(576, 224, 64, 96, 96, 128, 128, 2),
        InceptionA(576, 192, 96, 128, 96, 128, 128, 2),
        InceptionA(576, 160, 128, 160, 128, 160, 128, 2),
        InceptionA(608, 96, 128, 192, 160, 192, 128, 2),
        InceptionB(608, 128, 192, 192, 256),

        InceptionA(1056, 352, 192, 320, 160, 224, 128, 2),
        InceptionA(1024, 352, 192, 320, 192, 224, 128, 0),
        Pool(7, 1, 0, 2),
        Dense(1024, num_classes)
    )
end

"""
Builds a smaller version of Inception-BN network model modified to take 32x32 images.
This is the version used in the experiments.
"""
function create_inception_bn_sm_model(num_channels::Int, num_classes::Int, deeper=false)
    IA = deeper ? InceptionADeeper : InceptionA
    IB = deeper ? InceptionBDeeper : InceptionB

    Chain(
        ConvBN(3, 3, num_channels, 48),
        ConvBN(3, 3, 48, 96),

        IA(96, 26, 26, 26, 26, 35, 13, 2),
        IA(100, 26, 26, 35, 26, 35, 26, 2),
        IB(122, 53, 70, 26, 35),

        Pool(5, 3, 0, 2),
        ConvBN(1, 1, 227, 64),
        Flatten(1024),
        Linear(1024, num_classes)
    )
end


text:
"""
This module provides `tick()`, `tock()`, and `tok()` functions.

- `tick()`        start a new timer and start counting
- `tock()`        stop counting, show total elapsed time in canonical form
- `tok()`         stop counting, show and return total elapsed time in seconds
- `peektimer()    continue counting, show and return elapsed seconds so far
- `peektimers()   continue counting, show and return elapsed seconds so far
- `laptimer()     continue counting, show elapsed time so far in canonical form
- `alarm(h, m, s) set an alarm in h/m/s from now
- `alarm(dt)      set an alarm for `dt`
- `alarmlist()    list alarms

Don't use these for timing code execution! Julia provides
much better facilities for measuring performance, ranging
from the `@time` and `@elapsed` macros to packages such as
[BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl).
"""
module TickTock

export tick, tock, tok, peektimer, peektimers, laptimer, alarm, alarmlist

using Dates

const alarm_list = Tuple[]

function __init__()
    if ! haskey(ENV, "TICKTOCK_MESSAGES")
         ENV["TICKTOCK_MESSAGES"] = true
    end
end

"""
    tick()

Start a timer.

Other functions:

- `tock()` (stop counting and show canonical)

- `tok()` (stop counting and return seconds)

- `peektimer()` (continue counting, return elapsed seconds)

- `laptimer()` (continue counting, show canonical)

If `ENV["TICKTOCK_MESSAGES"]` = `false`, `tick()` will not display messages.
"""
function tick()
    t0 = time_ns()
    task_local_storage(:TIMERS, (t0, get(task_local_storage(), :TIMERS, ())))
    if ENV["TICKTOCK_MESSAGES"] == "true"
        @info " started timer at: $(now())"
    end
end

function printcanonical(tnow)
    canondc = Dates.canonicalize(
        Dates.CompoundPeriod(Dates.Second(floor(tnow)),
        Dates.Millisecond(floor((tnow - floor(tnow)) * 1000))))
    @info "$(lpad(tnow, 20))s: $canondc"
end

function gettimers()
    result = UInt64[]
    timers = get(task_local_storage(), :TIMERS, ())
    if timers === ()
        @warn "Use `tick()` to start a timer."
    else
        while timers[2] != ()
            push!(result, timers[1]::UInt64)
            timers = timers[2]
        end
        push!(result, timers[1]::UInt64)
        return result
    end
end

"""
    peektimer()

Return the elapsed seconds counted by the most recent timer, without stopping it.
"""
function peektimer()
    t1 = time_ns()
    timers = get(task_local_storage(), :TIMERS, ())
    if timers === ()
        @warn "Use `tick()` to start a timer."
    else
        t0 = timers[1]::UInt64
        return (t1 - t0)/1e9
    end
end

"""
    peektimers()

Return the elapsed seconds counted by all timers, without
stopping them, as an array.

## Example
```
julia> laptimer()
3         [ Info:         16.380142899s: 16 seconds, 380 milliseconds
2         [ Info:         21.480662472s: 21 seconds, 480 milliseconds
1         [ Info:         23.888862411s: 23 seconds, 888 milliseconds

julia> peektimers()
3-element Vector{Float64}:
 18.56985403
 23.670373603
 26.078573542
```
"""
function peektimers()
    t1 = time_ns()
    timers = gettimers()
    result = Float64[]
    if isnothing(timers)
        @warn "Use `tick()` to start a timer."
    else
        timernumber = length(timers)
        for t0 in timers
            tnow = (t1 - t0)/1e9
            push!(result, (t1 - t0)/1e9)
        end
    end
    if length(TickTock.alarm_list) > 1
        @info "show current alarms with `alarmlist()`"
    end
    if !isempty(result)
        return result
    end
end

"""
    tok()

Return the elapsed seconds counted by the most recent timer,
then stop counting.
"""
function tok()
    timers = get(task_local_storage(), :TIMERS, ())
    if timers === ()
        error("Use `tick()` to start a timer.")
    end
    t = peektimer()
    task_local_storage(:TIMERS, timers[2])
    return t
end

"""
    tock()

Print the elapsed time, in canonical form, of the most
recent timer, then stop counting.
"""
function tock()
    t = tok()
    printcanonical(t)
end

function showtimes(;canonical=true)
    t1 = time_ns()
    timers = gettimers()
    timernumber = length(timers)
    for t0 in timers
        tnow = (t1 - t0)/1e9
        print(rpad(timernumber, 10))
        canonical ? printcanonical(tnow) : @info(tnow)
        timernumber -= 1
    end
end

"""
    laptimer()

Print the elapsed time, in canonical form, of all active
timers, and continue counting.
"""
laptimer() = showtimes(canonical=true)

"""
    alarm(hours, minutes, seconds;
        action = () -> @info("TickTock: time's up"))

Set an alarm, with the option of providing a anonymous
function that executes when alarm fires.

```
using Dates

@async alarm(0, 5, 0, action = ()-> println("TickTock.jl: 5 minutes is up!"))

@async alarm(0, 0, 5, action = () -> run(`say "your alarm is ringing, sir"`), alarmname="tea's up") # uses macOS speech
```
"""
function alarm(hours, minutes, seconds;
        action=() -> @info("TickTock: alarm"),
        alarmname="TickTock alarm")
    tick();
    t = now()
    push!(TickTock.alarm_list, (
        "$(lpad(hour(t), 2, '0')):$(lpad(minute(t), 2, '0')):$(lpad(second(t), 2, '0'))",
        "$(lpad(hours, 2, '0')):$(lpad(minutes, 2, '0')):$(lpad(seconds, 2, '0'))",
        alarmname))
    while true
        sleep(1)
        if peektimer() > hours * 60 * 60 + minutes * 60 + seconds
            @info alarmname
            action()
            tock()
            break
        end
    end
end

"""
    alarm(dt::DateTime;
        action = () -> @info("TickTock: time's up"))

Run an alarm that fires at time `dt`, with the option of
providing a function that executes when alarm fires.

# Examples

```
using Dates

dt = now() + Dates.Minute(1)

@async alarm(dt, action = () -> println("TickTock.jl: Ready!"))
```

```
@async alarm(now() + Dates.Minute(10) + Dates.Second(0), action = () -> run(`say "wake up, 10 minutes is up"`)) # macOS speech command
```

"""
function alarm(dt::DateTime;
        action = () -> @info("TickTock: alarm"),
        alarmname = "TickTock alarm")
    p = Dates.Period(dt - now())
    secs = round(p, Dates.Second).value
    m, s = divrem(secs, 60)
    h, m = divrem(m, 60)
    alarm(h, m, s, action=action, alarmname=alarmname)
    @info "TickTock: \"$(alarmname)\" alarm for $h hours, $m minutes, $s seconds"
end

"""
@async alarm(now() + Dates.Second(5), action = () ->  TickTock.alarmnotify("time's up"))

MacOS only at the moment.
"""
function alarmnotify(subtitle="TickTock alarm")
    !Sys.isapple() && exit()
    command = """
    display notification with title "TickTock.jl" subtitle \"$(subtitle)\" sound name "frog"
    """
    chomp(read(`osascript -e $command`, String))
end

function alarmlist()
    if length(TickTock.alarm_list) > 0
        println("\nstart    | duration | finish     | name")
        for alarm in TickTock.alarm_list
            sh, sm, ss = parse.(Int, split(alarm[1], ":"))
            eh, em, es = parse.(Int, split(alarm[2], ":"))
            st = Time(sh, sm, ss)
            et = st + Hour(eh) + Minute(em) + Second(es)
            (Time(now()) > et) ? f = "✓" : f = " "
            println("$(alarm[1]) | $(alarm[2]) | $(Dates.format(et, "HH:MM:SS")) $f | $(alarm[3])" )
        end
    end
end

end # module

max_stars_repo_name: cormullion/TickTock

max_issues_repo_name: cormullion/TickTock

max_forks_repo_name: cormullion/TickTock


text:
struct Selector
    gid :: UInt64
    tag :: Symbol # :default, :invalid, :Gibbs, :HMC, etc.
    rerun :: Bool
end
function Selector(tag::Symbol = :default, rerun = tag != :default)
    return Selector(time_ns(), tag, rerun)
end
function Selector(gid::Integer, tag::Symbol = :default)
    return Selector(gid, tag, tag != :default)
end
hash(s::Selector) = hash(s.gid)
==(s1::Selector, s2::Selector) = s1.gid == s2.gid