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

using Markdown
using InteractiveUtils

# ╔═╡ 7e7dc954-866d-11eb-14df-c37b4fe266bf
using Calculus

# ╔═╡ 80e8c8fc-85c7-11eb-191f-23c76544a1ae
using Plots; import PlotlyJS; plotlyjs()

# ╔═╡ 8c0015c6-868d-11eb-1633-eb919d27d3ac
using PlutoUI

# ╔═╡ 0af7ed02-d4ab-4aab-8ed5-f03eb8cbadce
md"[Reference Chapter](https://www.deeplearningbook.org/contents/numerical.html)"

# ╔═╡ ddeed08a-8885-11eb-3681-fb6118238332
md"""
!!! example "Example 1"

	$y = x^4 - 3x^3 - 1$

	its derivative with respect to $x$ :

	$\frac{dy}{dx} = 4x^3 - 9x^2$

	equating them to zero to find their minimum:

	$x^2(4x - 9) = 0$

	We see that $y$ is minimum at $x=2.25$
"""

# ╔═╡ f123d468-8630-11eb-3daf-d35082dc5458
md"# Gradient Descent"

# ╔═╡ 946e6300-85c7-11eb-0e3d-abc305941caf
f(x) = x^4 - 3x^3 - 1

# ╔═╡ a4135446-85c7-11eb-1ca5-835e8460acd8
let
	x = -3:0.001:4.5
	plot(x, f)
end

# ╔═╡ c51e07a0-85c8-11eb-3247-c345ab59be62
function gradient_descent(α, x_0, f)
	x = copy(x_0)
	t = copy(x_0)
	df(x) = Calculus.gradient(f, x)
	for i = 1:500
		x = x - α * df(x)
		t = [t; x]
		if(round.(df(x), digits=5) == zeros(length(x_0)))
			break;
		end
	end
	return reshape(t, length(x_0), :)'
end

# ╔═╡ ed02fa0c-acf3-4a1a-b594-19747a69e757
md"""
#### What happened here?

I just took the gradient function from the book, and made it into a function in julia. It returns all the values of x, so we can can visualise the descent later.
"""

# ╔═╡ 2ba433f4-85cc-11eb-2043-0b75815bfec6
points = gradient_descent(0.01, 3.0, f)[1:end]

# ╔═╡ 4c51b95e-85cd-11eb-3432-47fc75f137eb
let
	x = -2:0.01:3
	plot(x, f)
	plot!(points, f, line = (:dash, :green), marker=(:circle,2))
end

# ╔═╡ 11908454-8630-11eb-1430-d37daf6d38d8
md"# Newton's Method"

# ╔═╡ 4e5269ea-8633-11eb-1a9e-db09365517f1
g(x) = x[1]^3 + x[2]^3 - 6x[1]*x[2]

# ╔═╡ 7408d648-8894-11eb-35f2-09ab2d5fbe52
md"""
!!! example "Example 2"

	$z = x^3 + y^3 - 6xy$

	its derivatives with respect to $x$ and $y$:

	$\frac{dz}{dx} = 3x^2 - 6y, \frac{dz}{dy} = 3y^2 - 6x$

	equating them to zero to find their minimum:

	$x^2 = 2y \;and\;y^2  = 2x$

	$\frac{y^4}{4} = 2y \;and\;\frac{x^4}{4} = 2x$

	$y^4 = 8y \;and\;x^4 = 8x$

	$y = 2 \;and\;x = 2$

	We see that $z$ is minimum at $x=y=2$
"""

# ╔═╡ 28a59930-8631-11eb-301d-6147e94972da
function newton(x_0, g)
	x = vec(Float64.(x_0))
	t = x_0
	for i in 1:500
		x = x - (inv(Calculus.hessian(g, x)) * Calculus.derivative(g, x))
		t = [t; x']
		if(round.(Calculus.derivative(g, x), digits=9) == zeros(length(x_0)))
			break
		end
	end
	return t
end

# ╔═╡ f8338402-866b-11eb-0d5c-9d66645ec5a2
newpoints = newton([4 4], g)

# ╔═╡ 5ed42dac-9016-11eb-3ba1-935c00ec9c99
f1(x,y) = x^3 + y^3 - 6x*y

# ╔═╡ c51c1cc8-8689-11eb-3bea-ad122713fa65
let
	x = y = -3:0.5:3
	surface(x, y, f1, fillalpha=0.5)
	plot!(newpoints[:,1], newpoints[:,2], f1.(newpoints[:,1],newpoints[:,2]), line=(:dot, 4), marker=(:circle, 1))
end

# ╔═╡ 6845fc55-81aa-423a-a08a-732ba0583a86
md"""
**Note: _You can find many such pre-defined methods in_** [Optim.jl](https://julianlsolvers.github.io/Optim.jl/stable/)
"""

# ╔═╡ ff78ac97-d8ba-45dd-abe5-e2518272de6c
md"In both the above methods, try setting the inital values to $-1.0$ for the Gradient Descent method, and $[-2, -2]$ for Newton's method. Lets see that below."

# ╔═╡ 2109fed6-b891-4be7-a078-cca774d9eed1
md"##### For Gradient Descent"

# ╔═╡ afbec578-7aa8-465d-8ab5-1cdcec6b4750
saddlepoints1 = gradient_descent(0.01, -1.0, f)

# ╔═╡ e1653b17-92e7-4b43-bc41-065d00d8c8f6
let
	x = -2:0.01:3
	plot(x, f, linecolor = :gray)
	plot!(saddlepoints1, f, linestyle = :dash, marker=(:circle,1), linecolor=:red)
end

# ╔═╡ ad1b745d-551d-4e9b-a943-577a1ab34dfb
md"##### For Newton's Method"

# ╔═╡ 842c8484-816d-42ac-be11-573cd39373a0
saddlepoints2 = newton([-2 -2], g)

# ╔═╡ cbc8ecb8-03c5-4e0f-aea1-0ad3807917d1
let
	x = y = -3:0.5:3
	surface(x, y, f1, fillalpha=0.5, c = :winter)
	plot!(saddlepoints2[:,1], saddlepoints2[:,2], f1.(saddlepoints2[:,1],saddlepoints2[:,2]), line=(:dash, :red, 2), marker=(:circle,1))
end

# ╔═╡ a60233dd-ef16-4101-a97e-f4db59ca67d8
md"""
Thats not the global minimum. Those points are called **Saddle Points**.
"""

# ╔═╡ 10a3c5bc-8cfa-11eb-0ff7-195329f7b738
md"## Try Yourself"

# ╔═╡ 87052830-9016-11eb-1082-111250655ad8
md"Try the linear least squares example from the book. Make the function and try plugging it into Newton's step method. As given, you should get your output within $1$ step. Also try to use Optim.jl to optimize this funciton again."

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
Calculus = "~0.5.1"
PlotlyJS = "~0.18.7"
Plots = "~1.21.3"
PlutoUI = "~0.7.9"
"""

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

julia_version = "1.7.0-beta4"
manifest_format = "2.0"

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

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

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

[[deps.AssetRegistry]]
deps = ["Distributed", "JSON", "Pidfile", "SHA", "Test"]
git-tree-sha1 = "b25e88db7944f98789130d7b503276bc34bc098e"
uuid = "bf4720bc-e11a-5d0c-854e-bdca1663c893"
version = "0.1.0"

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

[[deps.BinDeps]]
deps = ["Libdl", "Pkg", "SHA", "URIParser", "Unicode"]
git-tree-sha1 = "1289b57e8cf019aede076edab0587eb9644175bd"
uuid = "9e28174c-4ba2-5203-b857-d8d62c4213ee"
version = "1.0.2"

[[deps.Blink]]
deps = ["Base64", "BinDeps", "Distributed", "JSExpr", "JSON", "Lazy", "Logging", "MacroTools", "Mustache", "Mux", "Reexport", "Sockets", "WebIO", "WebSockets"]
git-tree-sha1 = "08d0b679fd7caa49e2bca9214b131289e19808c0"
uuid = "ad839575-38b3-5650-b840-f874b8c74a25"
version = "0.12.5"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

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

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

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

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

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

[[deps.FunctionalCollections]]
deps = ["Test"]
git-tree-sha1 = "04cb9cfaa6ba5311973994fe3496ddec19b6292a"
uuid = "de31a74c-ac4f-5751-b3fd-e18cd04993ca"
version = "0.5.0"

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

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

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

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

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

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

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

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

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

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

[[deps.Hiccup]]
deps = ["MacroTools", "Test"]
git-tree-sha1 = "6187bb2d5fcbb2007c39e7ac53308b0d371124bd"
uuid = "9fb69e20-1954-56bb-a84f-559cc56a8ff7"
version = "0.2.2"

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

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

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

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

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

[[deps.JSExpr]]
deps = ["JSON", "MacroTools", "Observables", "WebIO"]
git-tree-sha1 = "bd6c034156b1e7295450a219c4340e32e50b08b1"
uuid = "97c1335a-c9c5-57fe-bc5d-ec35cebe8660"
version = "0.5.3"

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

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

[[deps.Kaleido_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "2ef87eeaa28713cb010f9fb0be288b6c1a4ecd53"
uuid = "f7e6163d-2fa5-5f23-b69c-1db539e41963"
version = "0.1.0+0"

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

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

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

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

[[deps.Lazy]]
deps = ["MacroTools"]
git-tree-sha1 = "1370f8202dac30758f3c345f9909b97f53d87d3f"
uuid = "50d2b5c4-7a5e-59d5-8109-a42b560f39c0"
version = "0.15.1"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.Mustache]]
deps = ["Printf", "Tables"]
git-tree-sha1 = "36995ef0d532fe08119d70b2365b7b03d4e00f48"
uuid = "ffc61752-8dc7-55ee-8c37-f3e9cdd09e70"
version = "1.0.10"

[[deps.Mux]]
deps = ["AssetRegistry", "Base64", "HTTP", "Hiccup", "Pkg", "Sockets", "WebSockets"]
git-tree-sha1 = "82dfb2cead9895e10ee1b0ca37a01088456c4364"
uuid = "a975b10e-0019-58db-a62f-e48ff68538c9"
version = "0.7.6"

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

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

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

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

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

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

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

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

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

[[deps.Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "2276ac65f1e236e0a6ea70baff3f62ad4c625345"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.2"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "438d35d2d95ae2c5e8780b330592b6de8494e779"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.3"

[[deps.Pidfile]]
deps = ["FileWatching", "Test"]
git-tree-sha1 = "1be8660b2064893cd2dae4bd004b589278e4440d"
uuid = "fa939f87-e72e-5be4-a000-7fc836dbe307"
version = "1.2.0"

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

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

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

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

[[deps.PlotlyBase]]
deps = ["ColorSchemes", "Dates", "DelimitedFiles", "DocStringExtensions", "JSON", "LaTeXStrings", "Logging", "Parameters", "Pkg", "REPL", "Requires", "Statistics", "UUIDs"]
git-tree-sha1 = "7eb4ec38e1c4e00fea999256e9eb11ee7ede0c69"
uuid = "a03496cd-edff-5a9b-9e67-9cda94a718b5"
version = "0.8.16"

[[deps.PlotlyJS]]
deps = ["Base64", "Blink", "DelimitedFiles", "JSExpr", "JSON", "Kaleido_jll", "Markdown", "Pkg", "PlotlyBase", "REPL", "Reexport", "Requires", "WebIO"]
git-tree-sha1 = "ec6bc7270269be2d565b272116ca74ca2f8bd9ab"
uuid = "f0f68f2c-4968-5e81-91da-67840de0976a"
version = "0.18.7"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.URIParser]]
deps = ["Unicode"]
git-tree-sha1 = "53a9f49546b8d2dd2e688d216421d050c9a31d0d"
uuid = "30578b45-9adc-5946-b283-645ec420af67"
version = "0.4.1"

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

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

[[deps.UnPack]]
git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
version = "1.0.2"

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

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

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

[[deps.WebIO]]
deps = ["AssetRegistry", "Base64", "Distributed", "FunctionalCollections", "JSON", "Logging", "Observables", "Pkg", "Random", "Requires", "Sockets", "UUIDs", "WebSockets", "Widgets"]
git-tree-sha1 = "adc25e225bc334c7df6eec3b39496edfc451cc38"
uuid = "0f1e0344-ec1d-5b48-a673-e5cf874b6c29"
version = "0.8.15"

[[deps.WebSockets]]
deps = ["Base64", "Dates", "HTTP", "Logging", "Sockets"]
git-tree-sha1 = "f91a602e25fe6b89afc93cf02a4ae18ee9384ce3"
uuid = "104b5d7c-a370-577a-8038-80a2059c5097"
version = "1.5.9"

[[deps.Widgets]]
deps = ["Colors", "Dates", "Observables", "OrderedCollections"]
git-tree-sha1 = "eae2fbbc34a79ffd57fb4c972b08ce50b8f6a00d"
uuid = "cc8bc4a8-27d6-5769-a93b-9d913e69aa62"
version = "0.6.3"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# ╔═╡ Cell order:
# ╟─0af7ed02-d4ab-4aab-8ed5-f03eb8cbadce
# ╠═7e7dc954-866d-11eb-14df-c37b4fe266bf
# ╠═80e8c8fc-85c7-11eb-191f-23c76544a1ae
# ╟─ddeed08a-8885-11eb-3681-fb6118238332
# ╟─f123d468-8630-11eb-3daf-d35082dc5458
# ╠═946e6300-85c7-11eb-0e3d-abc305941caf
# ╠═a4135446-85c7-11eb-1ca5-835e8460acd8
# ╠═c51e07a0-85c8-11eb-3247-c345ab59be62
# ╟─ed02fa0c-acf3-4a1a-b594-19747a69e757
# ╠═2ba433f4-85cc-11eb-2043-0b75815bfec6
# ╠═4c51b95e-85cd-11eb-3432-47fc75f137eb
# ╟─11908454-8630-11eb-1430-d37daf6d38d8
# ╠═4e5269ea-8633-11eb-1a9e-db09365517f1
# ╟─7408d648-8894-11eb-35f2-09ab2d5fbe52
# ╠═28a59930-8631-11eb-301d-6147e94972da
# ╠═f8338402-866b-11eb-0d5c-9d66645ec5a2
# ╠═8c0015c6-868d-11eb-1633-eb919d27d3ac
# ╠═5ed42dac-9016-11eb-3ba1-935c00ec9c99
# ╠═c51c1cc8-8689-11eb-3bea-ad122713fa65
# ╟─6845fc55-81aa-423a-a08a-732ba0583a86
# ╟─ff78ac97-d8ba-45dd-abe5-e2518272de6c
# ╟─2109fed6-b891-4be7-a078-cca774d9eed1
# ╠═afbec578-7aa8-465d-8ab5-1cdcec6b4750
# ╠═e1653b17-92e7-4b43-bc41-065d00d8c8f6
# ╟─ad1b745d-551d-4e9b-a943-577a1ab34dfb
# ╠═842c8484-816d-42ac-be11-573cd39373a0
# ╠═cbc8ecb8-03c5-4e0f-aea1-0ad3807917d1
# ╟─a60233dd-ef16-4101-a97e-f4db59ca67d8
# ╟─10a3c5bc-8cfa-11eb-0ff7-195329f7b738
# ╟─87052830-9016-11eb-1082-111250655ad8
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: VarLad/MLTutorials

max_issues_repo_name: VarLad/MLTutorials

max_forks_repo_name: VarLad/MLTutorials


text:
module CFor

export @cfor

macro cfor(init, cond, inc, body)
    entry_label = gensym("entry")
    let_block = Expr(:block)
    while_block = Expr(:block)
    if !Meta.isexpr(init, :null) || (Meta.isexpr(inc, :block) && !isempty(inc.args))
        push!(while_block.args, inc)
        push!(while_block.args, Expr(:symboliclabel, entry_label))
    end
    if !Meta.isexpr(cond, :null) || (Meta.isexpr(cond, :block) && !isempty(cond.args))
        push!(while_block.args, Expr(:||, cond, Expr(:break)))
    end
    push!(while_block.args, body)

    if !Meta.isexpr(init, :null) || (Meta.isexpr(init, :block) && !isempty(init.args))
        push!(let_block.args, init)
    end
    push!(let_block.args, Expr(:symbolicgoto, entry_label))
    push!(let_block.args, Expr(:while, true, while_block))
    let_expr = Expr(:let, Expr(:block), let_block)
    return esc(let_expr)
end

end # module

ext: jl


text:
using Oceananigans
using Oceananigans.Units
using OceanTurbulenceParameterEstimation
using OceanTurbulenceParameterEstimation: Transformation
using LinearAlgebra, CairoMakie, DataDeps, Distributions

using ElectronDisplay

using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities:
    CATKEVerticalDiffusivity,
    SurfaceTKEFlux,
    MixingLength

#####
##### Compile LESbrary
#####

case_path(case) = @datadep_str("four_day_suite_1m/$(case)_instantaneous_statistics.jld2")

Δz = 2
#times = [96hours - 20minutes, 96hours]
times = [12hours, 96hours]
field_names = (:b, :e, :u, :v)
regrid_size = (1, 1, Int(256/Δz))

k₁ = floor(Int, 120 / Δz)
k₂ = ceil(Int, 150 / Δz)
space = SpaceIndices(z=k₁:k₂)
mult = length(space.z) / regrid_size[3]

transformation = (
                  b = Transformation(; normalization=ZScore()),
                  u = Transformation(; normalization=ZScore()),
                  v = Transformation(; normalization=ZScore()),
                  e = Transformation(; space, normalization=RescaledZScore(mult)),
                 )

observation_library = Dict()

# Don't optimize u, v for free_convection
observation_library["free_convection"] =
    SyntheticObservations(case_path("free_convection"); transformation, times, regrid_size,
                          field_names = (:b, :e))

# Don't optimize v for non-rotating cases
observation_library["strong_wind_no_rotation"] =
    SyntheticObservations(case_path("strong_wind_no_rotation"); transformation, times, regrid_size,
                          field_names = (:b, :e, :u))

# The rest are standard
for case in ["strong_wind", "strong_wind_weak_cooling", "weak_wind_strong_cooling"]
    observation_library[case] = SyntheticObservations(case_path(case); field_names, transformation, times, regrid_size)
end

@show znodes(Center, observation_library["free_convection"].grid)[space.z]

cases = [
         "free_convection",
         "weak_wind_strong_cooling",
         "strong_wind_weak_cooling",
         "strong_wind",
         "strong_wind_no_rotation",
        ]

observations = [observation_library[case] for case in cases]

#####
##### Simulation
#####

# Constant Ri, no convection
mixing_length = MixingLength(Cᴬu   = 0.0,
                             Cᴬc   = 0.0,
                             Cᴬe   = 0.0,
                             Cᴸᵇ   = 1.36,
                             Cᴷu⁻  = 0.101,
                             Cᴷc⁻  = 0.0574,
                             Cᴷe⁻  = 3.32,
                             Cᵟu   = 0.5,
                             Cᵟc   = 0.5,
                             Cᵟe   = 0.5,
                             CᴷRiᶜ = 2.0,
                             Cᴷuʳ  = 0.0,
                             Cᴷcʳ  = 0.0,
                             Cᴷeʳ  = 0.0)

surface_TKE_flux = SurfaceTKEFlux(CᵂwΔ=4.74, Cᵂu★=2.76)
catke = CATKEVerticalDiffusivity(; Cᴰ=1.78, mixing_length)

#####
##### Calibration
#####

mass = 0.6
prior_library = Dict()
prior_library[:CᵂwΔ]  = ScaledLogitNormal(; bounds=(4, 10)) #, interval=(2, 5), mass)
prior_library[:Cᵂu★]  = ScaledLogitNormal(; bounds=(2, 10)) #, interval=(3, 5), mass)
prior_library[:Cᴸᵇ]   = ScaledLogitNormal(; bounds=(0, 0.5)) #, interval=(0.1, 2), mass)
prior_library[:Cᴰ]    = ScaledLogitNormal(; bounds=(0, 1)) #, interval=(0.5, 2), mass)

prior_library[:Cᴷu⁻]  = ScaledLogitNormal(; bounds=(0, 0.1)) #, interval=(0.01, 0.1), mass)
prior_library[:Cᴷc⁻]  = ScaledLogitNormal(; bounds=(0, 1)) #, interval=(0.5, 1.0), mass)
prior_library[:Cᴷe⁻]  = ScaledLogitNormal(; bounds=(0, 3)) #, interval=(1.5, 3), mass)

prior_library[:Cᴷuʳ]  = ScaledLogitNormal(; bounds=(-1, 1))
prior_library[:Cᴷcʳ]  = ScaledLogitNormal(; bounds=(-1, 1))
prior_library[:Cᴷeʳ]  = ScaledLogitNormal(; bounds=(-1, 1))

prior_library[:CᴷRiᶜ] = ScaledLogitNormal(; bounds=(1, 3))
prior_library[:CᴷRiʷ] = ScaledLogitNormal(; bounds=(0, 0.5))

prior_library[:Cᴬu]   = ScaledLogitNormal(; bounds=(0, 0.1))
prior_library[:Cᴬc]   = ScaledLogitNormal(; bounds=(0, 2))
prior_library[:Cᴬe]   = ScaledLogitNormal(; bounds=(0, 0.1)) #, interval=(1, 3), mass)

prior_library[:Cᵟu]  = ScaledLogitNormal(; bounds=(0, 10))
prior_library[:Cᵟc]  = ScaledLogitNormal(; bounds=(0, 10))
prior_library[:Cᵟe]  = ScaledLogitNormal(; bounds=(0, 10))

# No convective adjustment:
constant_Ri_parameters = (:Cᴰ, :CᵂwΔ, :Cᵂu★, :Cᴸᵇ, :Cᴷu⁻, :Cᴷc⁻, :Cᴷe⁻, :Cᵟu, :Cᵟc, :Cᵟe)
variable_Ri_parameters = (:Cᴷuʳ, :Cᴷcʳ, :Cᴷeʳ, :CᴷRiʷ, :CᴷRiᶜ, :Cᴰ, :Cᴸᵇ, :CᵂwΔ, :Cᵂu★)
convective_adjustment_parameters = (:Cᴬc, :Cᴬe) # Cᴬu

# :Cᴸˢ
parameter_names = (:CᵂwΔ, :Cᵂu★, :Cᴷe⁻, :Cᴸᵇ, :Cᴰ, :Cᴷc⁻, :Cᴷu⁻, :Cᴷuʳ, :Cᴷcʳ, :Cᴷeʳ, :CᴷRiᶜ, :CᴷRiʷ, :Cᴬc, :Cᴬe)
free_parameters = FreeParameters(prior_library, names=parameter_names)

Nensemble = 4000
Δt = 20.0

function build_simulation()
    simulation = ensemble_column_model_simulation(observations;
                                                  Nensemble,
                                                  architecture = GPU(),
                                                  tracers = (:b, :e),
                                                  closure = catke)

    simulation.Δt = Δt

    Qᵘ = simulation.model.velocities.u.boundary_conditions.top.condition
    Qᵇ = simulation.model.tracers.b.boundary_conditions.top.condition
    N² = simulation.model.tracers.b.boundary_conditions.bottom.condition

    for (case, obs) in enumerate(observations)
        f = obs.metadata.parameters.coriolis_parameter
        view(Qᵘ, :, case) .= obs.metadata.parameters.momentum_flux
        view(Qᵇ, :, case) .= obs.metadata.parameters.buoyancy_flux
        view(N², :, case) .= obs.metadata.parameters.N²_deep
        view(simulation.model.coriolis, :, case) .= Ref(FPlane(f=f))
    end

    return simulation
end

simulation = build_simulation()
calibration = InverseProblem(observations, simulation, free_parameters)
resampler = Resampler(resample_failure_fraction=0.5, acceptable_failure_fraction=1.0)
eki = EnsembleKalmanInversion(calibration; resampler, convergence_rate=0.8)

#####
##### Plot utils
#####

Nt = length(times)
observed_data = []

for observation in observations
    time_serieses = observation.field_time_serieses
    names = keys(time_serieses)
    case_data = NamedTuple(n => interior(getproperty(time_serieses, n)[Nt])[1, 1, :] for n in names)
    push!(observed_data, case_data)
end

function get_modeled_case(icase, name, k=1)
    model_time_serieses = calibration.time_series_collector.field_time_serieses
    field = getproperty(model_time_serieses, name)[Nt]
    return view(interior(field), k, icase, :)
end

colorcycle =  [:black, :royalblue1, :darkgreen, :lightsalmon, :seagreen, :magenta2, :red4, :khaki1,   :darkgreen, :bisque4,
               :silver, :lightsalmon, :lightseagreen, :teal, :royalblue1, :darkorchid4]

markercycle = [:rect, :utriangle, :star5, :circle, :cross, :+, :pentagon, :ltriangle, :airplane, :diamond, :star4]
markercycle = repeat(markercycle, inner=2)

function make_axes(fig, row=1, label=nothing)
    ax_b = Axis(fig[row, 1], xlabel = "Buoyancy \n[cm s⁻²]", ylabel = "z [m]")
    ax_u = Axis(fig[row, 2], xlabel = "x-velocity \n[cm s⁻¹]")
    ax_v = Axis(fig[row, 3], xlabel = "y-velocity \n[cm s⁻¹]")
    ax_e = Axis(fig[row, 4], xlabel = "Turbulent kinetic energy \n[cm² s⁻²]")
    if !isnothing(label)
        ax_t = Axis(fig[row, 5])
        xlims!(0, 1)
        ylims!(0, 1)
        hidespines!(ax_t)
        hidedecorations!(ax_t)
        text!(ax_t, label, justification=:left, align=(:left, :center), position=(0, 0.5))
    end
    return (ax_b, ax_u, ax_v, ax_e)
end

function plot_fields!(axs, label, color, b, e, u=zeros(size(b)), v=zeros(size(b)); linewidth=2, linestyle=:solid)
    grid = first(values(observation_library)).grid
    z = znodes(Center, grid)
    b, u, v, e = Tuple(Array(f) for f in (b, u, v, e))

    for (q, name) in zip((b, u, v, e), ("b", "u", "v", "e"))
        any(isnan.(q)) && @warn("NaNs found in $label $(name)!")
    end

    ## Note unit conversions below, eg m s⁻¹ -> cm s⁻¹:cyan
    lines!(axs[1], 1e2 * b, z; color, linestyle, label, linewidth)
    lines!(axs[2], 1e2 * u, z; color, linestyle, label, linewidth)
    lines!(axs[3], 1e2 * v, z; color, linestyle, label, linewidth)
    lines!(axs[4], 1e4 * e, z; color, linestyle, label, linewidth)

    return nothing
end

function min_max_parameters(summary)
    names = keys(summary.ensemble_mean)
    Nens = length(summary.parameters)
    parameter_matrix = [summary.parameters[k][name] for name in names, k = 1:Nens]
    θ_min = minimum(parameter_matrix, dims=2)
    θ_max = maximum(parameter_matrix, dims=2)
    return θ_min, θ_max
end

function finitefind(a, val, find)
    b = deepcopy(a)
    b[.!isfinite.(a)] .= val
    return find(b)
end

finitefindmin(a) = finitefind(a, Inf, findmin)
finitefindmax(a) = finitefind(a, -Inf, findmax)

function visualize_parameter_evolution(eki)
    summaries = eki.iteration_summaries
    Niters = length(summaries)
    names = eki.inverse_problem.free_parameters.names
    θ_mean = NamedTuple(name => map(s -> s.ensemble_mean[name], summaries) for name in names)

    k_best(s) = finitefindmin(s.mean_square_errors)[2]
    θ_best = NamedTuple(name => map(s -> s.parameters[k_best(s)][name], summaries) for name in names)

    θ_min_max = [min_max_parameters(s) for s in summaries]
    θ_min = [[θn[1][i] for θn in θ_min_max] for i in 1:length(names)]
    θ_max = [[θn[2][i] for θn in θ_min_max] for i in 1:length(names)]

    θᵢ = NamedTuple(name => first(θ_mean[name]) for name in names)
    Δθ = NamedTuple(name => (θ_mean[name] .- θᵢ[name]) ./ θᵢ[name] for name in names)
    iterations = 0:length(summaries)-1

    fig = Figure(resolution=(1200, 1200))
    ax1 = Axis(fig[1:3, 1], xlabel = "Iteration", ylabel = "Δθ")
    for (i, name) in enumerate(names)
        label = string(name)
        marker = markercycle[i]
        color = colorcycle[i]
        scatterlines!(ax1, iterations, parent(Δθ[name]); marker, color=(color, 0.8), label, linewidth=4)
    end

    fig[1:3, 2] = Legend(fig, ax1)

    Nparts = 3
    Nθpart = floor(Int, length(names) / Nparts)

    for p in 1:Nparts
        axp = Axis(fig[p+3, 1], xlabel = "Iteration", ylabel = "θ")

        if p == Nparts
            np = UnitRange((p-1) * Nθpart + 1, length(names))
        else
            np = UnitRange((p-1) * Nθpart + 1, p * Nθpart)
        end

        partnames = names[np]

        for (n, name) in enumerate(partnames)
            i = (p - 1) * Nθpart + n
            label = string(name)
            marker = markercycle[i]
            color = colorcycle[i]
            scatterlines!(axp, iterations, parent(θ_mean[name]); marker, color=(color, 0.6), label, linewidth=4)
            lines!(axp, iterations, parent(θ_best[name]); color, linewidth=2)
            band!(axp, iterations, parent(θ_min[i]), parent(θ_max[i]), color=(color, 0.3))
        end

        fig[p+3, 2] = Legend(fig, axp)
    end

    display(fig)

    return nothing
end

function plot_latest(eki)
    latest_summary = eki.iteration_summaries[end]
    min_error, k_min = finitefindmin(latest_summary.mean_square_errors)
    max_error, k_max = finitefindmax(latest_summary.mean_square_errors)

    fig = Figure(resolution=(1200, 1200))

    for (c, case) in enumerate(cases)
        label = replace(case, "_" => "\n")
        axs = make_axes(fig, c, label)
        observed = observed_data[c]
        obs = observations[c]

        min_error_data = NamedTuple(n => get_modeled_case(c, n, k_min) for n in keys(obs.field_time_serieses))
        max_error_data = NamedTuple(n => get_modeled_case(c, n, k_max) for n in keys(obs.field_time_serieses))

        plot_fields!(axs, "observed at t = " * prettytime(times[end]), (:gray23, 0.6), observed...; linewidth=4)
        plot_fields!(axs, "min", :navy, min_error_data...)
        plot_fields!(axs, "max", :orangered3, max_error_data...)

        fig[1, 6] = Legend(fig, axs[1])
    end

    display(fig)

    return nothing
end

#####
##### Visualize observations
#####

fig = Figure(resolution=(1200, 1200))
linestyles = [:solid, :dash, :dot, :dashdot, :dashdotdot]
all_axs = []

for (o, observation) in enumerate(observations)
    axs_label = replace(cases[o], "_" => "\n")
    axs = make_axes(fig, o, axs_label)
    append!(all_axs, axs)
    for (n, t) in enumerate(times)
        linestyle = linestyles[o]
        label = "t = " * prettytime(t)
        names = keys(observation.field_time_serieses)
        data = map(name -> interior(observation.field_time_serieses[name][n])[1, 1, :], names)
        plot_fields!(axs, label, colorcycle[n], data...)
    end
end

display(fig)

#####
##### Calibrate
#####

# Initial state after 0 iterations
plot_latest(eki)

# Continuously update
for i = 1:20
    @info "Iterating..."
    start_time = time_ns()
    iterate!(eki)
    elapsed = 1e-9 * (time_ns() - start_time)
    @info string("   done. (", prettytime(elapsed), ")")
    @show eki.iteration_summaries[end]
    visualize_parameter_evolution(eki)
    plot_latest(eki)
end


text:
"""
```
plot_altpolicies(models, var, class, cond_type; kwargs...)

plot_altpolicies(models, vars, class, cond_type;
    forecast_string = "", altpol_string = "", untrans = false,
    fourquarter = false, plotroot = figurespath(m, \"forecast\"), titles = [],
    start_date = iterate_quarters(date_mainsample_end(models[1], -4)),
    end_date = iterate_quarters(date_forecast_start(models[1], 20)),
    kwargs...)
```

Plot `var` or `vars` forecasts under the alternative policies in `models`.

### Inputs

- `models::Vector{AbstractDSGEModel}`: vector of models, where each model has a
  different alternative policy
- `var::Symbol` or `vars::Vector{Symbol}`: variable(s) to be plotted,
  e.g. `:obs_gdp` or `[:obs_gdp, :obs_nominalrate]`
- `class::Symbol`
- `cond_type::Symbol`

### Keyword Arguments

- `forecast_string::String`
- `altpol_string::String`: identifies the group of alternative policies being
  plotted together. Required if `length(models) > 1`
- `start_date::Date`
- `end_date::Date`
- `untrans::Bool`: whether to plot untransformed forecast
- `fourquarter::Bool`: whether to plot four-quarter forecast
- `plotroot::String`: if nonempty, plots will be saved in that directory
- `title::String` or `titles::Vector{String}`

See `?histforecast` for additional keyword arguments, all of which can be passed
into `plot_altpolicies`.
"""
function plot_altpolicies(models::Vector{T}, var::Symbol, class::Symbol,
                          cond_type::Symbol; title::String = "",
                          kwargs...) where {T<:AbstractDSGEModel}
    plots = plot_altpolicies(models, [var], class, cond_type;
                             title = isempty(title) ? String[] : [title],
                             kwargs...)
    return plots[var]
end

function plot_altpolicies(models::Vector{T}, vars::Vector{Symbol}, class::Symbol,
                          cond_type::Symbol;
                          forecast_string::String = "",
                          altpol_string::String = "",
                          start_date::Date = iterate_quarters(date_mainsample_end(models[1]), -4),
                          end_date::Date = iterate_quarters(date_forecast_start(models[1]), 20),
                          untrans::Bool = false,
                          fourquarter::Bool = false,
                          plotroot::String = figurespath(m, "forecast"),
                          titles::Vector{String} = String[],
                          verbose::Symbol = :low,
                          kwargs...) where {T<:AbstractDSGEModel}
    # Determine output_vars
    if untrans && fourquarter
        error("Only one of untrans or fourquarter can be true")
    elseif untrans
        hist_prod  = :histut
        fcast_prod = :forecastut
    elseif fourquarter
        hist_prod  = :hist4q
        fcast_prod = :forecast4q
    else
        hist_prod  = :hist
        fcast_prod = :forecast
    end

    # Get titles if not provided
    if isempty(titles)
        detexify_title = typeof(Plots.backend()) == Plots.GRBackend
        titles = map(var -> describe_series(models[1], var, class, detexify = detexify_title))
    end

    # Temporarily set models[1] altpolicy key to altpol_string for get_forecast_filename
    n_altpolicies = length(models)
    if n_altpolicies != 1
        if isempty(altpol_string)
            error("Must provide nonempty altpol_key if plotting multiple alternative policies")
        else
            altpol = get_setting(models[1], :alternative_policy)
            old_key = altpol.key
            altpol.key = Symbol(altpol_string)
        end
    end

    # Loop through variables
    plots = OrderedDict{Symbol, Plots.Plot}()
    for (var, title) in zip(vars, titles)
        # Loop through alternative policies
        plots[var] = plot()
        for (i, m) in enumerate(models)
            # Get AltPolicy
            altpolicy = alternative_policy(m)

            # Read in MeansBands
            hist     = read_mb(m, :mode, cond_type, Symbol(hist_prod, class), forecast_string = forecast_string)
            forecast = read_mb(m, :mode, cond_type, Symbol(fcast_prod, class), forecast_string = forecast_string)

            # Call recipe
            plots[var] = histforecast!(var, hist, forecast;
                                       start_date = start_date, end_date = end_date,
                                       hist_label = "", forecast_label = string(altpolicy),
                                       forecast_color = altpolicy.color, linestyle = altpolicy.linestyle,
                                       ylabel = series_ylabel(m, var, class, untrans = untrans,
                                                              fourquarter = fourquarter),
                                       title = title, kwargs...)
        end

        # Save plot
        if !isempty(plotroot)
            output_file = get_forecast_filename(plotroot, filestring_base(models[1]), :mode, cond_type,
                                                Symbol("altpol", fcast_prod, "_", detexify(var)),
                                                forecast_string = forecast_string,
                                                fileformat = plot_extension())
            save_plot(plots[var], output_file, verbose = verbose)
        end
    end

    # Reset models[1] altpolicy key to original value
    if n_altpolicies != 1
        altpol.key = old_key
    end

    return plots
end


text:
@testset "TAI<->TT" begin
    ep = TAIEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TTEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 3
    @test second(Float64, out_ep) ≈ 2.3735247436378273e+01
    in_ep = TAIEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01
end
@testset "TAI<->UTC" begin
    ep = TAIEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = UTCEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 2
    @test second(Float64, out_ep) ≈ 1.4551247436378276e+01
    in_ep = TAIEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01
end
@testset "TAI<->TDB" begin
    ep = TAIEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TDBEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 3
    @test second(Float64, out_ep) ≈ 2.3734205844955390e+01 atol=1e-6
    in_ep = TAIEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-6
end
@testset "TAI<->TCB" begin
    ep = TAIEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TCBEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 3
    @test second(Float64, out_ep) ≈ 4.4097432701198270e+01 atol=1e-6
    in_ep = TAIEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-6
end
@testset "TAI<->TCG" begin
    ep = TAIEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TCGEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 3
    @test second(Float64, out_ep) ≈ 2.4650535580099750e+01
    in_ep = TAIEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01
end
@testset "TAI<->UT1" begin
    ep = TAIEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = UT1Epoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 2
    @test second(Float64, out_ep) ≈ 1.4617328996852756e+01 atol=1e-5
    in_ep = TAIEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-5
end
@testset "TT<->TAI" begin
    ep = TTEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TAIEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 2
    @test second(Float64, out_ep) ≈ 1.9367247436378280e+01
    in_ep = TTEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01
end
@testset "TT<->UTC" begin
    ep = TTEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = UTCEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 1
    @test second(Float64, out_ep) ≈ 4.2367247436378280e+01
    in_ep = TTEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01
end
@testset "TT<->TDB" begin
    ep = TTEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TDBEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 2
    @test second(Float64, out_ep) ≈ 5.1550205853115330e+01 atol=1e-6
    in_ep = TTEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-6
end
@testset "TT<->TCB" begin
    ep = TTEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TCBEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 3
    @test second(Float64, out_ep) ≈ 1.1913432210338932e+01 atol=1e-6
    in_ep = TTEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-6
end
@testset "TT<->TCG" begin
    ep = TTEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TCGEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 2
    @test second(Float64, out_ep) ≈ 5.2466535557669786e+01
    in_ep = TTEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01
end
@testset "TT<->UT1" begin
    ep = TTEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = UT1Epoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 1
    @test second(Float64, out_ep) ≈ 4.2433329223481640e+01 atol=1e-5
    in_ep = TTEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-5
end
@testset "UTC<->TAI" begin
    ep = UTCEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TAIEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 3
    @test second(Float64, out_ep) ≈ 2.8551247436378276e+01
    in_ep = UTCEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01
end
@testset "UTC<->TT" begin
    ep = UTCEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TTEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 4
    @test second(Float64, out_ep) ≈ 7.3524743637827330e-01
    in_ep = UTCEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01
end
@testset "UTC<->TDB" begin
    ep = UTCEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TDBEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 4
    @test second(Float64, out_ep) ≈ 7.3420583557444980e-01 atol=1e-6
    in_ep = UTCEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-6
end
@testset "UTC<->TCB" begin
    ep = UTCEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TCBEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 4
    @test second(Float64, out_ep) ≈ 2.1097433265509650e+01 atol=1e-6
    in_ep = UTCEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-6
end
@testset "UTC<->TCG" begin
    ep = UTCEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TCGEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 4
    @test second(Float64, out_ep) ≈ 1.6505356058861196e+00
    in_ep = UTCEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01
end
@testset "UTC<->UT1" begin
    ep = UTCEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = UT1Epoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 2
    @test second(Float64, out_ep) ≈ 5.1617328736316390e+01 atol=1e-5
    in_ep = UTCEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-5
end
@testset "TDB<->TAI" begin
    ep = TDBEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TAIEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 2
    @test second(Float64, out_ep) ≈ 1.9368289019641487e+01 atol=1e-6
    in_ep = TDBEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-6
end
@testset "TDB<->TT" begin
    ep = TDBEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TTEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 2
    @test second(Float64, out_ep) ≈ 5.1552289019641485e+01 atol=1e-6
    in_ep = TDBEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-6
end
@testset "TDB<->UTC" begin
    ep = TDBEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = UTCEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 1
    @test second(Float64, out_ep) ≈ 4.2368289019641490e+01 atol=1e-6
    in_ep = TDBEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-6
end
@testset "TDB<->TCB" begin
    ep = TDBEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TCBEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 3
    @test second(Float64, out_ep) ≈ 1.1914473793618030e+01 atol=1e-6
    in_ep = TDBEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-6
end
@testset "TDB<->TCG" begin
    ep = TDBEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TCGEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 2
    @test second(Float64, out_ep) ≈ 5.2467577140933720e+01 atol=1e-6
    in_ep = TDBEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-6
end
@testset "TDB<->UT1" begin
    ep = TDBEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = UT1Epoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 1
    @test second(Float64, out_ep) ≈ 4.2434370806737520e+01 atol=1e-5
    in_ep = TDBEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-5
end
@testset "TCB<->TAI" begin
    ep = TCBEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TAIEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 1
    @test second(Float64, out_ep) ≈ 5.9005062972974656e+01 atol=1e-6
    in_ep = TCBEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-6
end
@testset "TCB<->TT" begin
    ep = TCBEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TTEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 2
    @test second(Float64, out_ep) ≈ 3.1189062972974654e+01 atol=1e-6
    in_ep = TCBEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-6
end
@testset "TCB<->UTC" begin
    ep = TCBEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = UTCEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 1
    @test second(Float64, out_ep) ≈ 2.2005062972974656e+01 atol=1e-6
    in_ep = TCBEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-6
end
@testset "TCB<->TDB" begin
    ep = TCBEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TDBEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 2
    @test second(Float64, out_ep) ≈ 3.1188021394874360e+01 atol=1e-6
    in_ep = TCBEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-6
end
@testset "TCB<->TCG" begin
    ep = TCBEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TCGEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 2
    @test second(Float64, out_ep) ≈ 3.2104351080075170e+01 atol=1e-6
    in_ep = TCBEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-6
end
@testset "TCB<->UT1" begin
    ep = TCBEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = UT1Epoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 1
    @test second(Float64, out_ep) ≈ 2.2071144903463825e+01 atol=1e-5
    in_ep = TCBEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-5
end
@testset "TCG<->TAI" begin
    ep = TCGEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TAIEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 2
    @test second(Float64, out_ep) ≈ 1.8451959315724658e+01
    in_ep = TCGEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01
end
@testset "TCG<->TT" begin
    ep = TCGEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TTEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 2
    @test second(Float64, out_ep) ≈ 5.0635959315724655e+01
    in_ep = TCGEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01
end
@testset "TCG<->UTC" begin
    ep = TCGEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = UTCEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 1
    @test second(Float64, out_ep) ≈ 4.1451959315724660e+01
    in_ep = TCGEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01
end
@testset "TCG<->TDB" begin
    ep = TCGEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TDBEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 2
    @test second(Float64, out_ep) ≈ 5.0634917732693770e+01 atol=1e-6
    in_ep = TCGEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-6
end
@testset "TCG<->TCB" begin
    ep = TCGEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TCBEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 3
    @test second(Float64, out_ep) ≈ 1.0998144075725655e+01 atol=1e-6
    in_ep = TCGEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-6
end
@testset "TCG<->UT1" begin
    ep = TCGEpoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = UT1Epoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 1
    @test second(Float64, out_ep) ≈ 4.1518041109273234e+01 atol=1e-5
    in_ep = TCGEpoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-5
end
@testset "UT1<->TAI" begin
    ep = UT1Epoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TAIEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 3
    @test second(Float64, out_ep) ≈ 2.8485166135974850e+01 atol=1e-5
    in_ep = UT1Epoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-5
end
@testset "UT1<->TT" begin
    ep = UT1Epoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TTEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 4
    @test second(Float64, out_ep) ≈ 6.6916613597484800e-01 atol=1e-5
    in_ep = UT1Epoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-5
end
@testset "UT1<->UTC" begin
    ep = UT1Epoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = UTCEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 2
    @test second(Float64, out_ep) ≈ 5.1485166135974850e+01 atol=1e-5
    in_ep = UT1Epoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-5
end
@testset "UT1<->TDB" begin
    ep = UT1Epoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TDBEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 4
    @test second(Float64, out_ep) ≈ 6.6812453518777910e-01 atol=1e-5
    in_ep = UT1Epoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-5
end
@testset "UT1<->TCB" begin
    ep = UT1Epoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TCBEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 4
    @test second(Float64, out_ep) ≈ 2.1031351964098377e+01 atol=1e-5
    in_ep = UT1Epoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-5
end
@testset "UT1<->TCG" begin
    ep = UT1Epoch(2018, 8, 14, 10, 2, 5.1551247436378276e+01)
    out_ep = TCGEpoch(ep)
    @test year(out_ep) == 2018
    @test month(out_ep) == 8
    @test day(out_ep) == 14
    @test hour(out_ep) == 10
    @test minute(out_ep) == 4
    @test second(Float64, out_ep) ≈ 1.5844543054366440e+00 atol=1e-5
    in_ep = UT1Epoch(out_ep)
    @test year(in_ep) == 2018
    @test month(in_ep) == 8
    @test day(in_ep) == 14
    @test hour(in_ep) == 10
    @test minute(in_ep) == 2
    @test second(Float64, in_ep) ≈ 5.1551247436378276e+01 atol=1e-5
end


text:
import TimeZones: UTCOffset, value, isdst

@test value(UTCOffset(0, 0)) == 0
@test value(UTCOffset(3600, 0)) == 3600
@test value(UTCOffset(0, 3600)) == 3600
@test value(UTCOffset(-7200, 3600)) == -3600

@test string(UTCOffset(0, 0)) == "+00:00"
@test string(UTCOffset(3600, 0)) == "+01:00"
@test string(UTCOffset(0, 3600)) == "+01:00"
@test string(UTCOffset(-7200, 3600)) == "-01:00"
@test string(UTCOffset(-3661)) == "-01:01:01"

@test !isdst(UTCOffset(0))
@test !isdst(UTCOffset(0, 0))
@test isdst(UTCOffset(0, 3600))
@test isdst(UTCOffset(0, 7200))

# Arithmetic
dt = DateTime(2015, 1, 1)
@test dt + UTCOffset(7200, -3600) == dt + Hour(1)
@test dt - UTCOffset(7200, -3600) == dt - Hour(1)

# Comparisons
let a = UTCOffset(0, 0), b = UTCOffset(0, 0)
    @test a == b
    @test isequal(a, b)
end

let a = UTCOffset(0, 3600), b = UTCOffset(3600, 0)
    @test a != b
    @test !isequal(a, b)
end

let a = UTCOffset(7200, 3600), b = UTCOffset(3600, 7200)
    @test a != b

    # Treated as equal since some time zone representations adjust the standard
    # offset when dealing with midsummer time.
    @test isequal(a, b)
end

buffer = IOBuffer()
show(buffer, UTCOffset(0, 0))
@test takebuf_string(buffer) == "UTC+0/+0"
show(buffer, UTCOffset(3600, 7200))
@test takebuf_string(buffer) == "UTC+1/+2"


text:
module TestCore

using AnisotropicDiffusion
using AnisotropicDiffusion: generic_step!, hand_tuned_step!
using ImageFiltering
using LinearAlgebra
using Test

@testset "1d edge" begin
    img = zeros(20)
    img[1:10] .= 100
    img_noisy = img + randn(size(img)...)
    img_clean = denoise(img_noisy, PeronaMalik(niter=3, lambda=0.25));
    ker = centered(normalize(ones(3),1))
    img_baseline = imfilter(img_noisy, ker)
    @test sum(abs2, img - img_clean) < sum(abs2, img - img_baseline)
end

@testset "generic vs hand tuned $(dim)d" for dim in [1,2]
    alg1 = PeronaMalik(step = generic_step!, niter=1)
    alg2 = PeronaMalik(step = hand_tuned_step!, niter=1)
    img = randn([10,11,12][1:dim]...)
    ret1 = denoise(img, alg1)
    ret2 = denoise(img, alg2)
    @test ret1 ≈ ret2
    @test ret1 !== ret2
end

end#module


text:
struct Linear{BC<:Union{Throw{OnGrid},Periodic{OnCell}}} <: DegreeBC{1}
    bc::BC
    function Linear{BC}(bc::BC=BC()) where BC<:Union{Throw{OnGrid},Periodic{OnCell}}
        new{BC}(bc)
    end
end

Linear() = Linear(Throw(OnGrid()))
Linear(::Periodic{Nothing}) = Linear(Periodic(OnCell()))
Linear(bc::BC) where BC<:Union{Throw{OnGrid},Periodic{OnCell}} = Linear{BC}(bc)

function Base.show(io::IO, deg::Linear{Throw{OnGrid}})
    print(io, nameof(typeof(deg)), '(', ')')
end

"""
    Linear()

Indicate that the corresponding axis should use linear interpolation.

# Extended help

Assuming uniform knots with spacing 1, the `i`th piece of linear b-spline
implemented here is defined as follows.

    y_i(x) = c p(x) + cp p(1-x)

where

    p(δx) = x

and the values `cX` for `X ∈ {_, p}` are the coefficients.

Linear b-splines are naturally interpolating, and require no prefiltering;
there is therefore no need for boundary conditions to be provided.

Also, although the implementation is slightly different in order to re-use
the framework built for general b-splines, the resulting interpolant is just
a piecewise linear function connecting each pair of neighboring data points.
"""
Linear

function positions(deg::Linear, ax::AbstractUnitRange{<:Integer}, x)
    f = floor(x)
    # When x == last(ax) we want to use the x-1, x pair
    f = ifelse(x == last(ax), f - oneunit(f), f)
    fi = fast_trunc(Int, f)
    expand_index(deg, fi, ax), x-f
end
expand_index(::Linear{Throw{OnGrid}}, fi::Number, ax::AbstractUnitRange) = fi
expand_index(::Linear{Periodic{OnCell}}, fi::Number, ax::AbstractUnitRange) =
    (modrange(fi, ax), modrange(fi+1, ax))

value_weights(::Linear, δx) = (1-δx, δx)
gradient_weights(::Linear, δx) = (-oneunit(δx), oneunit(δx))
hessian_weights(::Linear, δx) = (zero(δx), zero(δx))

padded_axis(ax::AbstractUnitRange, ::BSpline{<:Linear}) = ax


text:
export nyquist

#==#
#
# @brief Compute nyquist diagram.
# @param[in] sys System in state space representation.
# @param[in] w Array of frequencies in which the bode diagram will be computed.
# If empty, the frequencies will be chosen according to the system poles.
# @param[in] plot Plot bode diagram (default = true).
#
# @return Real(G(jw)), Imag(G(jw)), w (freq. array) for each pair input/output.
#
#==#

"""
### function nyquist(sys::StateSpace, w::Array=[]; plot=true)

Compute the Nyquist diagram.

##### Args

* sys: System in state space representation.
* w: Array of frequencies in which the bode diagram will be computed.  If empty,
  the frequencies will be chosen according to the system poles.

##### Returns

Real(G(jw)), Imag(G(jw)), w (freq. array) for each input/output pair.

"""

function nyquist(sys::StateSpace, w::Array=[]; plot=true)
    A, B, C, D = sys.A, sys.B, sys.C, sys.D

    # Number of inputs.
    num_u = size(B,2)

    # Number of outputs.
    num_y = size(C,1)

    # Number of states.
    num_states = size(A,1)

    # Identity matrix.
    eye_s = eye(num_states)

    # Check if the vector of frequencies was provided.
    w_user = !isempty(w)

    w_final = 0.0

    tol = 0.0

    # If the user did not provide a vector with frequencies, then they will be
    # computed automatically.
    if (!w_user)
        # If not, compute the bode diagram from a frequency two decades lower
        # than the minimum frequency and two decades higher than the maximum
        # frequency.
        poles      = pole(sys)
        real_poles = abs.(real(poles[poles .!= 0]))

        # Get maximum frequency of the system.
        highest_freq = (isempty(real_poles)) ? 0.0 : maximum(real_poles)
        w_final = (highest_freq != 0) ? highest_freq*100.0 : 1e+2

        # Compute the tolerance using the DC gain.
        dc_gain = maximum(D + C*pinv(eye(A)-A)*B)
        tol     = dc_gain/50.0

        # Don't let tolerance to be 0.
        (tol == 0) && (tol = 1e-2)

        # Initialize the vector `w` with 100 points between the minimum and
        # maximum frequencies.
        w = collect(linspace(0,w_final,100))
    end

    # Frequency response for each pair input/output.
    real_fr     = Vector{Matrix{Float64}}(0)
    imag_fr     = Vector{Matrix{Float64}}(0)
    mir_imag_fr = Vector{Matrix{Float64}}(0)

    # Compute the Nyquist plot.
    i  = 1
    wi = w[1]

    while i <= length(w)
        wi = w[i]

        s = Complex(0,wi)

        fr = C*pinv(s*eye_s-A)*B+D

        real_fr_i  = Float64[real(f) for f in fr]
        imag_fr_i  = Float64[imag(f) for f in fr]

        # If w == 0 and A is singular, then we have a pole at the origin. In
        # this case, we must not consider the frequencies computed to plot the
        # graph.
        if (wi == 0.0) && (det(A) == 0.0)
            real_fr_i = ones(real_fr_i)*Inf
            imag_fr_i = ones(imag_fr_i)*Inf
        end

        if w_user || (i == 1)
            i += 1
            push!(real_fr,     +real_fr_i)
            push!(imag_fr,     +imag_fr_i)

            (!w_user) && push!(mir_imag_fr, -imag_fr_i)
        else
            # Verify the maximum distance between the current point and the last
            # one.
            Δx = maximum(abs.(real_fr_i - real_fr[i-1]))
            Δy = maximum(abs.(imag_fr_i - imag_fr[i-1]))

            # Accept the point if the tolerance is higher than the distances.
            if (maximum([Δx;Δy]) > tol) && ( (w[i]-w[i-1]) > 1e-2 )
                if i == 2
                    w = [collect(linspace(w[1],w[2],10)); w[3:end]]
                elseif i == length(w)
                    w = [w[1:i-2]; collect(linspace(w[i-1],w[i],10))]
                else
                    w = [w[1:i-2]; collect(linspace(w[i-1],w[i],10)); w[i+1:end]]
                end
            else
                i += 1
                push!(real_fr,     +real_fr_i)
                push!(imag_fr,     +imag_fr_i)
                push!(mir_imag_fr, -imag_fr_i)
            end
        end
    end

    ( plot ) && plot_nyquist(real_fr, imag_fr, num_u, num_y)

    real_fr, imag_fr, w
end

#==#
#
# @brief Compute bode diagram.
# @param[in] G Transfer function of the dynamic system.
# @param[in] w Array of frequencies in which the bode diagram will be computed.
# If empty, the frequencies will be chosen according to the system poles.
# @param[in] plot Plot bode diagram (default = true).
#
# @return Real(G(jw)), Imag(G(jw)), w (freq. array) for each pair input/output.
#
#==#

function nyquist(G::TransferFunction, w::Array=[]; plot=true)
    # Convert the transfer function to state space and call bode.jl.
    nyquist(tf2ss(G), w; plot=plot)
end


text:
include("comm_manager.jl")
include("execute_request.jl")

using .CommManager

# Don't send previous lines to the completions function,
# due to issue #380.  Find the start of the first line
# (if any) where the expression is parseable.  Replace
# with find_parsestart(c,p) = start(c) once julia#9467 is merged.
parseok(s) = !Meta.isexpr(parse(s, raise=false), :error)
function find_parsestart(code, cursorpos)
    s = start(code)
    while s < cursorpos
        parseok(code[s:cursorpos]) && return s
        s = nextind(code, s)
        while s < cursorpos && code[s] ∉ ('\n','\r')
            s = nextind(code, s)
        end
    end
    return start(code) # failed to find parseable lines
end

function complete_request(socket, msg)
    code = msg.content["code"]
    cursor_chr = msg.content["cursor_pos"]
    cursorpos = cursor_chr <= 0 ? 0 : chr2ind(code, cursor_chr)
    if isspace(code[1:cursorpos])
        send_ipython(requests[], msg_reply(msg, "complete_reply",
                                 Dict("status" => "ok",
                                              "matches" => String[],
                                              "cursor_start" => cursor_chr,
                                              "cursor_end" => cursor_chr)))
        return
    end

    codestart = find_parsestart(code, cursorpos)
    comps, positions = Base.REPLCompletions.completions(code[codestart:end], cursorpos-codestart+1)
    positions += codestart-1
    if isempty(positions) # true if comps to be inserted without replacement
        cursor_start = (cursor_end = ind2chr(code, last(positions)))
    else
        cursor_start = ind2chr(code, first(positions)) - 1
        cursor_end = ind2chr(code, last(positions))
    end
    send_ipython(requests[], msg_reply(msg, "complete_reply",
                                     Dict("status" => "ok",
                                                  "matches" => comps,
                                                  "cursor_start" => cursor_start,
                                                  "cursor_end" => cursor_end)))
end

function kernel_info_request(socket, msg)
    send_ipython(requests[],
                 msg_reply(msg, "kernel_info_reply",
                           Dict("protocol_version" => "5.0",
                                        "implementation" => "ijulia",
                                        # TODO: "implementation_version" => IJulia version string from Pkg
                                        "language_info" =>
                                        Dict("name" => "julia",
                                             "version" =>
                                             string(VERSION.major, '.',
                                                    VERSION.minor, '.',
                                                    VERSION.patch),
                                             "mimetype" => "application/julia",
                                             "file_extension" => ".jl"),
                                        "banner" => "Julia: A fresh approach to technical computing.",
                                        "help_links" => [
                                                         Dict("text"=>"Julia Home Page",
                                                              "url"=>"http://julialang.org/"),
                                                         Dict("text"=>"Julia Documentation",
                                                              "url"=>"http://docs.julialang.org/"),
                                                         Dict("text"=>"Julia Packages",
                                                              "url"=>"http://pkg.julialang.org/")
                                                        ])))
end

function connect_request(socket, msg)
    send_ipython(requests[],
                 msg_reply(msg, "connect_reply",
                           Dict("shell_port" => profile["shell_port"],
                                        "iopub_port" => profile["iopub_port"],
                                        "stdin_port" => profile["stdin_port"],
                                        "hb_port" => profile["hb_port"])))
end

function shutdown_request(socket, msg)
    send_ipython(requests[], msg_reply(msg, "shutdown_reply",
                                     msg.content))
    sleep(0.1) # short delay (like in ipykernel), to hopefully ensure shutdown_reply is sent
    exit()
end

stripdots(ex) = :_
stripdots(ex::Symbol) = ex
stripdots(ex::Expr) = Meta.isexpr(ex, :.) ? stripdots(ex.args[2]) : :_
stripdots(ex::QuoteNode) = ex.value
stripdots(ex::GlobalRef) = ex.name
rm_sideeffects(ex) = ex
function rm_sideeffects(ex::Expr)
    if Meta.isexpr(ex, :call)
        name = stripdots(ex.args[1])
        if name == :repl_search || name == :repl_corrections
            return nothing
        else
            return ex
        end
    else
        return Expr(ex.head, map(rm_sideeffects, ex.args)...)
    end
end
function docdict(s::AbstractString)
    ex = macroexpand(parse(helpcode(s)))
    # unfortunately, the REPL help macros sometimes have
    # expressions with side effects (I/O), so we need to
    # remove these.
    display_dict(eval(Main, rm_sideeffects(ex)))
end
import Base: is_id_char, is_id_start_char
function get_token(code, pos)
    # given a string and a cursor position, find substring to request
    # help on by:
    #   1) searching backwards, skipping invalid identifier chars
    #        ... search forward for end of identifier
    #   2) search backwards to find the biggest identifier (including .)
    #   3) if nothing found, do return empty string
    # TODO: detect operators?

    startpos = pos
    while startpos > start(code)
        if is_id_char(code[startpos])
            break
        else
            startpos = prevind(code, startpos)
        end
    end
    endpos = startpos
    while startpos >= start(code) && (is_id_char(code[startpos]) || code[startpos] == '.')
        startpos = prevind(code, startpos)
    end
    startpos = startpos < pos ? nextind(code, startpos) : pos
    if !is_id_start_char(code[startpos])
        return ""
    end
    while endpos < endof(code) && is_id_char(code[endpos])
        endpos = nextind(code, endpos)
    end
    if !is_id_char(code[endpos])
        endpos = prevind(code, endpos)
    end
    return code[startpos:endpos]
end

function inspect_request(socket, msg)
    try
        code = msg.content["code"]
        s = get_token(code, chr2ind(code, msg.content["cursor_pos"]))
        if isempty(s)
            content = Dict("status" => "ok", "found" => false)
        else
            d = docdict(s)
            content = Dict("status" => "ok",
                           "found" => !isempty(d),
                           "data" => d)
        end
        send_ipython(requests[], msg_reply(msg, "inspect_reply", content))
    catch e
        content = error_content(e, backtrace_top=:inspect_request);
        content["status"] = "error"
        send_ipython(requests[],
                     msg_reply(msg, "inspect_reply", content))
    end
end

function history_request(socket, msg)
    # we will just send back empty history for now, pending clarification
    # as requested in ipython/ipython#3806
    send_ipython(requests[],
                 msg_reply(msg, "history_reply",
                           Dict("history" => [])))
end

function is_complete_request(socket, msg)
    ex = parse(msg.content["code"], raise=false)
    status = Meta.isexpr(ex, :incomplete) ? "incomplete" : Meta.isexpr(ex, :error) ? "invalid" : "complete"
    send_ipython(requests[],
                 msg_reply(msg, "is_complete_reply",
                           Dict("status"=>status, "indent"=>"")))
end

const handlers = Dict{String,Function}(
    "execute_request" => execute_request,
    "complete_request" => complete_request,
    "kernel_info_request" => kernel_info_request,
    "inspect_request" => inspect_request,
    "connect_request" => connect_request,
    "shutdown_request" => shutdown_request,
    "history_request" => history_request,
    "is_complete_request" => is_complete_request,
    "comm_open" => comm_open,
    "comm_info_request" => comm_info_request,
    "comm_msg" => comm_msg,
    "comm_close" => comm_close
)