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

using Markdown
using InteractiveUtils

# ╔═╡ 88d2a9b8-d8e4-47a2-b1de-369d8fb44c20
using PlutoUI; PlutoUI.TableOfContents(title = "Benchmarking Procesos")

# ╔═╡ 07ddd182-470c-11ec-3704-93171dbc510e
begin
	using CSV, DataFrames, DataFramesMeta, Dates, IndexedTables
	using FileIO, JLD2, JuliaDB, JuliaDBMeta, Statistics, StatsBase
	using Plots, StatsPlots
	#cd(raw"\\srvfs02\Departamento de Investigacion Economica\Departamento\BD Internacional\Remesas")
	cd(raw"D:\Departamento\BD Internacional\Remesas")
	#include("functions/fn_all.jl");
	Plots.plotly()
	#ENV["COLUMNS"]=1000
	#ENV["LINES"] = 60

	## Cargar base desde archivo, como dataframe
	#directorio = "//srvfs02/Departamento de Investigacion Economica/Departamento/BD Internacional/Remesas/"
	directorio = "D:/Departamento/BD Internacional/Remesas/"
end

# ╔═╡ d1d0d4b2-92c1-4faa-8e42-7b7241fa1cc2
begin
	# 801197710224
	x = round.(rand(5)*10^12)#00000000000)
	using Formatting
	format.(x)
	#x
	round.(x./10^11)
end

# ╔═╡ 454637cc-589c-459d-93fd-7da3ba6fd67c
md"# Leer CSVs"

# ╔═╡ 53266a68-82e6-4024-af9d-8cefe89399f9
global function mergeCSV(fecha_proceso)
	n1 = Dates.now()
	searchdir(path, key) = Base.filter(
		x -> Base.contains(x, key), Base.readdir(path))
	key = fecha_proceso * ".csv"
	path = "csv/"
	csv = searchdir(path, key)
	NArchivo_csv = Base.size(csv)[1];

	# Importar archivos, seleccionar y renombrar variables
	#path * csv[1]
	bases_csv = []

	for i in 1:NArchivo_csv
		Base.push!(bases_csv,
			CSV.read(path * csv[i],
				delim = ',',
				DataFrames.DataFrame))
		DataFrames.rename!(bases_csv[i],
			:Fecha_Contable => :Fecha,
			:Departamento => :DeptoTransacc,
			:Municipio => :MunicTransacc,
			:RTN_Identidad => :ID)
		bases_csv[i] = DataFramesMeta.select(bases_csv[i],
			:Fecha,
			#:Agente_Cambiario,
			#:Tipo_Transaccion,
			#:Actividad_Economica,
			#:Medio_Pago,
			:DeptoTransacc,
			:MunicTransacc,
			:Tasa_Cambio,
			:Monto,
			:ID,
			:Tipo_Cliente);
	end

	# Unir en una sola tabla
	database = Base.copy(bases_csv[1])
	for i in 2:NArchivo_csv
		database = Base.vcat(
			database, bases_csv[i])
	end
	database
end

# ╔═╡ 4a5b02b5-83d5-46b9-88ec-41f4d0df54d1
database = mergeCSV("2021")

# ╔═╡ ae2c7408-7ca2-4166-98be-f9ba3d8476ed
md"# Cambiar Formatos"

# ╔═╡ 38b7bde3-7c34-4930-8582-bfd067a05c09
md"## Fecha"

# ╔═╡ cd85e83f-f380-480a-9c78-5fdbbba068ab
begin
	database[!, :Fechas] .= Dates.Date("1-1-1")
	database[!, :Days] .= 0
    database[!, :Weeks] .= 0
    database[!, :Months] .= 0
    database[!, :Years] .= 0
	database
end

# ╔═╡ 5cb827be-5ad3-49a3-b920-b1df8a8cb861

# ╔═╡ df699c60-d076-4319-9d95-13ea2afe955c

# ╔═╡ fbd5cb7c-e722-4a2e-84ef-4c6702134234
# Benchmarking
#=
@time database[!, :Fechas] = Dates.Date.(
		database[!, :Fecha],
		Dates.DateFormat("d/m/Y H:M:S"))
@time for i in 1:Base.size(database)[1]
	database[i, :Fechas] = Dates.Date.(
		database[i, :Fecha],
		Dates.DateFormat("d/m/Y H:M:S"))
	end
=#

# ╔═╡ e5552459-e2ad-4f44-9760-757ea37fa233
@time for i in 1:Base.size(database)[1]
	database[i, :Fechas] = Dates.Date.(
		database[i, :Fecha],
		Dates.DateFormat("d/m/Y H:M:S"))
end

# ╔═╡ 50ca420b-2a11-4366-9541-555585eeb5ea

# ╔═╡ b45cf1f6-c70c-450b-a8e4-29f7eb44e791
database

# ╔═╡ e3c415b1-0964-4337-9f9f-c1a9616a8b5d

# ╔═╡ 029a721d-4a55-4cc0-8898-61bdfb80c10c
md"## Depurar IDs, solo enteros"

# ╔═╡ 7ffa2874-0690-4587-91db-f0e55570bcb1
global function StrToInt(x)
	try
		Base.parse.(Int64, x)
	catch
		return 0
	end
end

# ╔═╡ 11a22c40-1da4-4f9c-a650-000230546bfa
md"## Tipo de Cliente"

# ╔═╡ 3e06e476-3c9d-42c4-b167-c2ce926b2c3c
begin
	database[!, :Tipo_Clientes] .= 0
	for i in 1:Base.size(database)[1]
		if database[i, :Tipo_Cliente] == "Natural Residente"
			database[i, :Tipo_Clientes] = 1
		else
			database[i, :Tipo_Clientes] = 0
		end
	end
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
DataFramesMeta = "1313f7d8-7da2-5740-9ea0-a2ca25f37964"
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
Formatting = "59287772-0a20-5a39-b81b-1366585eb4c0"
IndexedTables = "6deec6e2-d858-57c5-ab9b-e6ca5bd20e43"
JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
JuliaDB = "a93385a2-3734-596a-9a66-3cfbb77141e6"
JuliaDBMeta = "2c06ca41-a429-545c-b8f0-5ca7dd64ba19"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"

[compat]
CSV = "~0.8.5"
DataFrames = "~0.22.7"
DataFramesMeta = "~0.8.0"
FileIO = "~1.11.2"
Formatting = "~0.4.2"
IndexedTables = "~0.8.1"
JLD2 = "~0.4.3"
JuliaDB = "~0.9.0"
JuliaDBMeta = "~0.4.0"
Plots = "~0.28.4"
PlutoUI = "~0.7.1"
StatsBase = "~0.32.2"
StatsPlots = "~0.14.5"
"""

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

[[AbstractFFTs]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "485ee0867925449198280d4af84bdb46a2a404d0"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.0.1"

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

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

[[Arpack]]
deps = ["Arpack_jll", "Libdl", "LinearAlgebra"]
git-tree-sha1 = "2ff92b71ba1747c5fdd541f8fc87736d82f40ec9"
uuid = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
version = "0.4.0"

[[Arpack_jll]]
deps = ["Libdl", "OpenBLAS_jll", "Pkg"]
git-tree-sha1 = "e214a9b9bd1b4e1b4f15b22c0994862b66af7ff7"
uuid = "68821587-b530-5797-8361-c406ea357684"
version = "3.5.0+3"

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

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

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

[[BinaryProvider]]
deps = ["Libdl", "Logging", "SHA"]
git-tree-sha1 = "ecdec412a9abc8db54c0efc5548c64dfce072058"
uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
version = "0.5.10"

[[CSV]]
deps = ["Dates", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode"]
git-tree-sha1 = "b83aa3f513be680454437a0eee21001607e5d983"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.8.5"

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

[[CategoricalArrays]]
deps = ["DataAPI", "Future", "JSON", "Missings", "Printf", "Statistics", "StructTypes", "Unicode"]
git-tree-sha1 = "18d7f3e82c1a80dd38c16453b8fd3f0a7db92f23"
uuid = "324d7699-5711-5eae-9e2f-1d82baa6b597"
version = "0.9.7"

[[Chain]]
git-tree-sha1 = "cac464e71767e8a04ceee82a889ca56502795705"
uuid = "8be319e6-bccf-4806-a6f7-6fae938471bc"
version = "0.4.8"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "f885e7e7c124f8c92650d61b9477b9ac2ee607dd"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.11.1"

[[ChangesOfVariables]]
deps = ["LinearAlgebra", "Test"]
git-tree-sha1 = "9a1d594397670492219635b35a3d830b04730d62"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.1"

[[Clustering]]
deps = ["Distances", "LinearAlgebra", "NearestNeighbors", "Printf", "SparseArrays", "Statistics", "StatsBase"]
git-tree-sha1 = "75479b7df4167267d75294d14b58244695beb2ac"
uuid = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5"
version = "0.14.2"

[[CodecZlib]]
deps = ["BinaryProvider", "Libdl", "TranscodingStreams"]
git-tree-sha1 = "05916673a2627dd91b4969ff8ba6941bc85a960e"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.6.0"

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

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "InteractiveUtils", "Printf", "Reexport"]
git-tree-sha1 = "177d8b959d3c103a6d57574c38ee79c81059c31b"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.11.2"

[[CommonSubexpressions]]
deps = ["MacroTools", "Test"]
git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
version = "0.3.0"

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

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

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

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[Dagger]]
deps = ["Distributed", "LinearAlgebra", "MemPool", "Profile", "Random", "Serialization", "SharedArrays", "SparseArrays", "Statistics", "StatsBase", "Test"]
git-tree-sha1 = "39dd9351c6637a71d1925e26e49adcd8a25ebf0b"
uuid = "d58978e5-989f-55fb-8d15-ea34adc7bf54"
version = "0.8.0"

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

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

[[DataFramesMeta]]
deps = ["Chain", "DataFrames", "MacroTools", "Reexport"]
git-tree-sha1 = "55be907a531471de062f321147006c04b8c1e75f"
uuid = "1313f7d8-7da2-5740-9ea0-a2ca25f37964"
version = "0.8.0"

[[DataStructures]]
deps = ["InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "88d48e133e6d3dd68183309877eac74393daa7eb"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.17.20"

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

[[DataValues]]
deps = ["DataValueInterfaces", "Dates"]
git-tree-sha1 = "d88a19299eba280a6d062e135a43f00323ae70bf"
uuid = "e7dc6d0d-1eca-5fa6-8ad6-5aecde8b7ea5"
version = "0.4.13"

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

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

[[DiffEqDiffTools]]
deps = ["LinearAlgebra", "SparseArrays", "StaticArrays"]
git-tree-sha1 = "b992345a39b4d9681342ae795a8dacc100730182"
uuid = "01453d9d-ee7c-5054-8395-0335cb756afa"
version = "0.14.0"

[[DiffResults]]
deps = ["StaticArrays"]
git-tree-sha1 = "c18e98cba888c6c25d1c3b048e4b3380ca956805"
uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
version = "1.0.3"

[[DiffRules]]
deps = ["LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
git-tree-sha1 = "3287dacf67c3652d3fed09f4c12c187ae4dbb89a"
uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
version = "1.4.0"

[[Distances]]
deps = ["LinearAlgebra", "Statistics", "StatsAPI"]
git-tree-sha1 = "837c83e5574582e07662bbbba733964ff7c26b9d"
uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
version = "0.10.6"

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

[[Distributions]]
deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"]
git-tree-sha1 = "9c41285c57c6e0d73a21ed4b65f6eec34805f937"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.23.8"

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

[[DoubleFloats]]
deps = ["GenericLinearAlgebra", "LinearAlgebra", "Polynomials", "Printf", "Quadmath", "Random", "Requires", "SpecialFunctions"]
git-tree-sha1 = "79ce76892ef75c65b210fc531891e2427152d4fd"
uuid = "497a8b3b-efae-58df-a0af-a86822472b78"
version = "1.1.24"

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

[[ExprTools]]
git-tree-sha1 = "b7e3d17636b348f005f11040025ae8c6f645fe92"
uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
version = "0.1.6"

[[FFMPEG]]
deps = ["BinaryProvider", "Libdl"]
git-tree-sha1 = "9143266ba77d3313a4cf61d8333a1970e8c5d8b6"
uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
version = "0.2.4"

[[FFTW]]
deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"]
git-tree-sha1 = "463cb335fa22c4ebacfd1faba5fde14edb80d96c"
uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
version = "1.4.5"

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

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

[[FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays"]
git-tree-sha1 = "4863cbb7910079369e258dee4add9d06ead5063a"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.8.14"

[[FixedPointNumbers]]
git-tree-sha1 = "d14a6fa5890ea3a7e5dcab6811114f132fec2b4b"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.6.1"

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

[[ForwardDiff]]
deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
git-tree-sha1 = "6406b5112809c08b1baa5703ad274e1dded0652f"
uuid = "f6369f11-7733-5829-9624-2563aa707210"
version = "0.10.23"

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

[[GR]]
deps = ["Base64", "DelimitedFiles", "LinearAlgebra", "Printf", "Random", "Serialization", "Sockets", "Test"]
git-tree-sha1 = "c690c2ab22ac9ee323d9966deae61a089362b25c"
uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
version = "0.44.0"

[[GenericLinearAlgebra]]
deps = ["LinearAlgebra", "Printf", "Random"]
git-tree-sha1 = "ac44f4f51ffee9ff1ea50bd3fbb5677ea568d33d"
uuid = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
version = "0.2.7"

[[GeometryTypes]]
deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "StaticArrays"]
git-tree-sha1 = "78f0ce9d01993b637a8f28d84537d75dc0ce8eef"
uuid = "4d00f742-c7ba-57c2-abde-4428a4b178cb"
version = "0.7.10"

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

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

[[IndexedTables]]
deps = ["DataValues", "Dates", "IteratorInterfaceExtensions", "LinearAlgebra", "OnlineStats", "PooledArrays", "Random", "Serialization", "SparseArrays", "Statistics", "TableTraits", "TableTraitsUtils", "Test", "WeakRefStrings"]
git-tree-sha1 = "cb4406eeed8729f8179940d39e6dec5598971341"
uuid = "6deec6e2-d858-57c5-ab9b-e6ca5bd20e43"
version = "0.8.1"

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

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

[[Interpolations]]
deps = ["AxisAlgorithms", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"]
git-tree-sha1 = "2b7d4e9be8b74f03115e64cf36ed2f48ae83d946"
uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
version = "0.12.10"

[[Intervals]]
deps = ["Dates", "Printf", "RecipesBase", "Serialization", "TimeZones"]
git-tree-sha1 = "323a38ed1952d30586d0fe03412cde9399d3618b"
uuid = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
version = "1.5.0"

[[InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "a7254c0acd8e62f1ac75ad24d5db43f5f19f3c65"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.2"

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

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

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

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

[[JLD2]]
deps = ["CodecZlib", "DataStructures", "MacroTools", "Mmap", "Pkg", "Printf", "Requires", "UUIDs"]
git-tree-sha1 = "9f2f2f24e60305feb6ae293a617ddf06f429efc3"
uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
version = "0.4.3"

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

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

[[JuliaDB]]
deps = ["Dagger", "DataValues", "Dates", "Distributed", "Glob", "IndexedTables", "MemPool", "Nullables", "OnlineStats", "OnlineStatsBase", "PooledArrays", "Printf", "Random", "RecipesBase", "Serialization", "Statistics", "StatsBase", "Test", "TextParse", "WeakRefStrings"]
git-tree-sha1 = "83cfc7566e079721ab9c2b19efa6a26bc5ba9d20"
uuid = "a93385a2-3734-596a-9a66-3cfbb77141e6"
version = "0.9.0"

[[JuliaDBMeta]]
deps = ["Dagger", "Distributed", "IndexedTables", "IterTools", "JuliaDB", "MacroTools", "Reexport", "Test"]
git-tree-sha1 = "ca0c4c156d4a3b8055dcabd36eacb3d7469b8583"
uuid = "2c06ca41-a429-545c-b8f0-5ca7dd64ba19"
version = "0.4.0"

[[KernelDensity]]
deps = ["Distributions", "FFTW", "Interpolations", "Optim", "StatsBase", "Test"]
git-tree-sha1 = "c1048817fe5711f699abc8fabd47b1ac6ba4db04"
uuid = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
version = "0.5.1"

[[LazyArtifacts]]
deps = ["Artifacts", "Pkg"]
uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"

[[LearnBase]]
deps = ["LinearAlgebra", "SparseArrays", "StatsBase", "Test"]
git-tree-sha1 = "c4b5da6d68517f46f70ed5157b28336b56cd2ff3"
uuid = "7f8f8fb0-2700-5f03-b4bd-41f8cfc144b6"
version = "0.2.2"

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

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

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

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

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

[[LineSearches]]
deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"]
git-tree-sha1 = "f27132e551e959b3667d8c93eae90973225032dd"
uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
version = "7.1.1"

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

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

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

[[LossFunctions]]
deps = ["InteractiveUtils", "LearnBase", "Markdown", "Random", "RecipesBase", "SparseArrays", "Statistics", "StatsBase", "Test"]
git-tree-sha1 = "08d87fec43e7d335811dfae5b55dbfc5690e915b"
uuid = "30fc2ffe-d236-52d8-8643-a9d8f7c094a7"
version = "0.5.1"

[[MKL_jll]]
deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
git-tree-sha1 = "5455aef09b40e5020e1520f551fa3135040d4ed0"
uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
version = "2021.1.1+2"

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

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

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

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

[[MemPool]]
deps = ["DataStructures", "Distributed", "Mmap", "Random", "Serialization", "Sockets", "Test"]
git-tree-sha1 = "d52799152697059353a8eac1000d32ba8d92aa25"
uuid = "f9f48841-c794-520a-933b-121f7ba6ed94"
version = "0.2.0"

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

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

[[Mocking]]
deps = ["Compat", "ExprTools"]
git-tree-sha1 = "29714d0a7a8083bba8427a4fbfb00a540c681ce7"
uuid = "78c3b35d-d492-501b-9361-3d52fe80e533"
version = "0.7.3"

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

[[MultivariateStats]]
deps = ["Arpack", "LinearAlgebra", "SparseArrays", "Statistics", "StatsBase"]
git-tree-sha1 = "352fae519b447bf52e6de627b89f448bcd469e4e"
uuid = "6f286f6a-111f-5878-ab1e-185364afe411"
version = "0.7.0"

[[MutableArithmetics]]
deps = ["LinearAlgebra", "SparseArrays", "Test"]
git-tree-sha1 = "7bb6853d9afec54019c1397c6eb610b9b9a19525"
uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
version = "0.3.1"

[[NLSolversBase]]
deps = ["Calculus", "DiffEqDiffTools", "DiffResults", "Distributed", "ForwardDiff"]
git-tree-sha1 = "f1b8ed89fa332f410cfc7c937682eb4d0b361521"
uuid = "d41bc354-129a-5804-8e4c-c37616107c6c"
version = "7.5.0"

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

[[NearestNeighbors]]
deps = ["Distances", "StaticArrays"]
git-tree-sha1 = "16baacfdc8758bc374882566c9187e785e85c2f0"
uuid = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
version = "0.4.9"

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

[[Nullables]]
git-tree-sha1 = "8f87854cc8f3685a60689d8edecaa29d2251979b"
uuid = "4d1e1d77-625e-5b40-9113-a560ec7a8ecd"
version = "1.0.0"

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

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

[[OnlineStats]]
deps = ["DataStructures", "Dates", "LearnBase", "LinearAlgebra", "LossFunctions", "OnlineStatsBase", "PenaltyFunctions", "Random", "RecipesBase", "Reexport", "Statistics", "StatsBase", "SweepOperator", "Test"]
git-tree-sha1 = "a783fe25075b9d09d6dd2ba67dafd374d835acfa"
uuid = "a15396b6-48d5-5d58-9928-6d29437db91e"
version = "0.23.0"

[[OnlineStatsBase]]
deps = ["Dates", "LearnBase", "Statistics", "Test"]
git-tree-sha1 = "2d62fb8d9b1b407f9f80e99f4a63f32dea90423c"
uuid = "925886fa-5bf2-5e8e-b522-a9147a512338"
version = "0.10.2"

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

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

[[Optim]]
deps = ["Compat", "FillArrays", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"]
git-tree-sha1 = "c05aa6b694d426df87ff493306c1c5b4b215e148"
uuid = "429524aa-4258-5aef-a3af-852621145aeb"
version = "0.22.0"

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

[[PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse", "Test"]
git-tree-sha1 = "95a4038d1011dfdbde7cecd2ad0ac411e53ab1bc"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.10.1"

[[Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.3"

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

[[PenaltyFunctions]]
deps = ["InteractiveUtils", "LearnBase", "LinearAlgebra", "RecipesBase", "Reexport", "Test"]
git-tree-sha1 = "b0baaa5218ca0ffd6a8ae37ef0b58e0df688ac8b"
uuid = "06bb1623-fdd5-5ca2-a01c-88eae3ea319e"
version = "0.1.2"

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

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

[[PlotUtils]]
deps = ["Colors", "Dates", "Printf", "Random", "Reexport"]
git-tree-sha1 = "51e742162c97d35f714f9611619db6975e19384b"
uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
version = "0.6.5"

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

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "Logging", "Markdown", "Random", "Suppressor"]
git-tree-sha1 = "45ce174d36d3931cd4e37a47f93e07d1455f038d"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.1"

[[Polynomials]]
deps = ["Intervals", "LinearAlgebra", "MutableArithmetics", "RecipesBase"]
git-tree-sha1 = "79bcbb379205f1c62913fa9ebecb413c7a35f8b0"
uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
version = "2.0.18"

[[PooledArrays]]
deps = ["DataAPI"]
git-tree-sha1 = "b1333d4eced1826e15adbdf01a4ecaccca9d353c"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "0.5.3"

[[PositiveFactorizations]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "17275485f373e6673f7e7f97051f703ed5b15b20"
uuid = "85a6dd25-e78a-55b7-8502-1745935b8125"
version = "0.2.4"

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

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

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

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

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

[[Quadmath]]
deps = ["Printf", "Random", "Requires"]
git-tree-sha1 = "5a8f74af8eae654086a1d058b4ec94ff192e3de0"
uuid = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd"
version = "0.5.5"

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

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

[[Ratios]]
deps = ["Requires"]
git-tree-sha1 = "01d341f502250e81f6fec0afe662aa861392a3aa"
uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439"
version = "0.4.2"

[[RecipesBase]]
git-tree-sha1 = "7bdce29bc9b2f5660a6e5e64d64d91ec941f6aa2"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "0.7.0"

[[Reexport]]
deps = ["Pkg"]
git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "0.2.0"

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

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

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

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

[[SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "f45b34656397a1f6e729901dc9ef679610bd12b5"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.8"

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

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

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

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

[[SortingAlgorithms]]
deps = ["DataStructures", "Random", "Test"]
git-tree-sha1 = "03f5898c9959f8115e30bc7226ada7d0df554ddd"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "0.3.1"

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

[[SpecialFunctions]]
deps = ["OpenSpecFun_jll"]
git-tree-sha1 = "d8d8b8a9f4119829410ecd706da4cc8594a1e020"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "0.10.3"

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

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

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

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

[[StatsFuns]]
deps = ["Rmath", "SpecialFunctions"]
git-tree-sha1 = "ced55fd4bae008a8ea12508314e725df61f0ba45"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.7"

[[StatsPlots]]
deps = ["Clustering", "DataStructures", "DataValues", "Distributions", "Interpolations", "KernelDensity", "MultivariateStats", "Observables", "Plots", "RecipesBase", "Reexport", "StatsBase", "TableOperations", "Tables", "Widgets"]
git-tree-sha1 = "388e6588a1e09e763b0f38c498aeb350b5f76828"
uuid = "f3b207a7-027a-5e70-b257-86293d7955fd"
version = "0.14.5"

[[StructTypes]]
deps = ["Dates", "UUIDs"]
git-tree-sha1 = "d24a825a95a6d98c385001212dc9020d609f2d4f"
uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
version = "1.8.1"

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

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

[[SweepOperator]]
deps = ["LinearAlgebra", "Test"]
git-tree-sha1 = "2039aaa96f7b21634a1b3246cb2699dd1d26d473"
uuid = "7522ee7d-7047-56d0-94d9-4bc626e7058d"
version = "0.2.0"

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

[[TableOperations]]
deps = ["Tables", "Test"]
git-tree-sha1 = "208630a14884abd110a8f8008b0882f0d0f5632c"
uuid = "ab02a1b2-a7df-11e8-156e-fb1833f50b87"
version = "0.2.1"

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

[[TableTraitsUtils]]
deps = ["DataValues", "IteratorInterfaceExtensions", "Missings", "TableTraits"]
git-tree-sha1 = "78fecfe140d7abb480b53a44f3f85b6aa373c293"
uuid = "382cd787-c1b6-5bf2-a167-d5b971a19bda"
version = "1.0.2"

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

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

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

[[TextParse]]
deps = ["CodecZlib", "DataStructures", "Dates", "DoubleFloats", "Mmap", "Nullables", "PooledArrays", "WeakRefStrings"]
git-tree-sha1 = "26b43d6746b52cca13c4cdef90f89652273b413e"
uuid = "e0df1984-e451-5cb5-8b61-797a481e67e3"
version = "0.9.1"

[[TimeZones]]
deps = ["Dates", "Future", "LazyArtifacts", "Mocking", "Pkg", "Printf", "RecipesBase", "Serialization", "Unicode"]
git-tree-sha1 = "a5688ffdbd849a98503c6650effe79fe89a41252"
uuid = "f269a46b-ccf7-5d73-abea-4c690281aa53"
version = "1.5.9"

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

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

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

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

[[WeakRefStrings]]
deps = ["Missings", "Random", "Test"]
git-tree-sha1 = "960639a12ffd223ee463e93392aeb260fa325566"
uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
version = "0.5.8"

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

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

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

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

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

# ╔═╡ Cell order:
# ╠═88d2a9b8-d8e4-47a2-b1de-369d8fb44c20
# ╠═07ddd182-470c-11ec-3704-93171dbc510e
# ╠═d1d0d4b2-92c1-4faa-8e42-7b7241fa1cc2
# ╟─454637cc-589c-459d-93fd-7da3ba6fd67c
# ╠═53266a68-82e6-4024-af9d-8cefe89399f9
# ╠═4a5b02b5-83d5-46b9-88ec-41f4d0df54d1
# ╟─ae2c7408-7ca2-4166-98be-f9ba3d8476ed
# ╟─38b7bde3-7c34-4930-8582-bfd067a05c09
# ╠═cd85e83f-f380-480a-9c78-5fdbbba068ab
# ╠═5cb827be-5ad3-49a3-b920-b1df8a8cb861
# ╠═df699c60-d076-4319-9d95-13ea2afe955c
# ╠═fbd5cb7c-e722-4a2e-84ef-4c6702134234
# ╠═e5552459-e2ad-4f44-9760-757ea37fa233
# ╠═50ca420b-2a11-4366-9541-555585eeb5ea
# ╠═b45cf1f6-c70c-450b-a8e4-29f7eb44e791
# ╠═e3c415b1-0964-4337-9f9f-c1a9616a8b5d
# ╟─029a721d-4a55-4cc0-8898-61bdfb80c10c
# ╠═7ffa2874-0690-4587-91db-f0e55570bcb1
# ╟─11a22c40-1da4-4f9c-a650-000230546bfa
# ╠═3e06e476-3c9d-42c4-b167-c2ce926b2c3c
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: ElvisCasco/Remesas

max_issues_repo_name: ElvisCasco/Remesas

max_forks_repo_name: ElvisCasco/Remesas


text:
export TR_Nwt

function TR_Nwt(h :: AbstractNLPModel,
				nlpstop :: AbstractStopping;
				kwargs...)

	optimal, stop = TR_generic(h, nlpstop; direction = :Nwt, kwargs...)
	return optimal, stop
end

ext: jl


text:
using GeoPhyInv
using Test
using SparseArrays

# # Intro
#md # ```@docs
#md # AGeom
#md # ```

# # Examples

# Lets create a `mgrid` for the experiment.
mgrid = [range(0.0, stop=10.,step=0.1), range(0.0, stop=30.,step=0.2)];

# Then initialize an acquisition on `mgrid`.

ageom=AGeom(mgrid, SSrcs(10), Srcs(10), Recs(10));

# Otherwise, initialize with one of the predefined acquisitions.

ageom=AGeom(mgrid, :xwell, SSrcs(10), Recs(10));

# By default, the sources and receivers are randomly placed
# on the grid.
# If unsure, test it.
@test (ageom ∈ mgrid)

# The source and receiver positions can be updated as desired.

update!(ageom[1], Srcs(), [0,1], [10,20],);
update!(ageom[1], Recs(), [0,0], [10,20],);
update!(ageom, SSrcs(), [0,1], [10,20], );

# It is easy to combine supersources. Now `ageom2` has 20 supersources.
ageom2=vcat(ageom, ageom);

# # Methods
#md # ```@docs
#md # update!(::AGeomss, ::Srcs)
#md # update!(::AGeomss, ::Recs)
#md # update!(::AGeom, ::SSrcs)
#md # update!(::AGeom, ::Recs)
#md # Base.in(::AGeom, ::AbstractVector{StepRangeLen})
#md # Base.isequal(::AGeom, ::AGeom)
#md # SparseArrays.SparseMatrixCSC(::AGeomss,::AbstractVector{StepRangeLen})
#md # ```


text:
using InteractionNetworkModels, BenchmarkTools, Distributions
using Plots, StatsBase

mcmc = SisMcmcSplitMerge(
    β=0.3, p=0.7, η=0.5,
    ν_td=1
)

model_mode = Hollywood(-3.0, Poisson(7), 10)
S = [[1,1,1,1], [2,2,2,2], [3,3,3,3]]
V = 1:20
# d = FastEditDistance(FastLCS(51),51)
d = DTW(FastLCS(100))
model = SIS(S, 3.9, d, V, 50, 50)

S2 = [[1,1,1,1], [2,2,3,11,2,1,2], [2,2,2,2], [3,3,3,3]]
d(S,S2)

x = mcmc(model, desired_samples=100)
plot(x)
model.mode

S_prop = deepcopy(S)
inds = [2]
multiple_adj_split_noisy!(
    S_prop,
    [Int[] for i in 1:100],
    inds,
    Geometric(0.7),
    0.7, 1:10,
    1,10
)
S_prop

# Base dimension range performance
# --------------------------------
xl,xu  = (1,100)
r = DimensionRange(xl,xu)
y=10
@btime notin(r,y)

function test(x::Real, l::Real, u::Real)
    return (x >= l) & (x<=u)
end
@btime test(y, xl, xu)

# But do this and we get performace back
rxl = r.l
rxu = r.u
@btime test(y, rxl, rxu)

# Undo a merge

function test_undo_merge!(
    x::Path{Int},
    ind::Vector{Int}
    )
    live_index = 0
    for i in ind
        x[i+live_index] = 0
        insert!(x, i+live_index+1, 0)
        live_index += 1
    end
end

x = [1,2,3,4,5]
test_undo_merge!(x, [1,5])
x

function test_undo_split!(
    x::Path{Int},
    ind::Vector{Int}
    )

    for i in ind
        tmp = popat!(x, i+1)
        x[i] = 0
    end

end

x = [1,1,2,3,4,5,5]
ind=[1,5]

test_undo_split!(x, [1,5])
x

function get_cond_dist(
    S::InteractionSequence{Int};
    α::Float64=0.0,
    V::Int=length(unique(vcat(S...)))
    )
    tot_len = 0
    c = fill(α, V)
    for path in S
        for v in path
            c[v] += 1.0
            tot_len += 1
        end
    end
    c ./= (tot_len + α * V)
    return c
end

get_cond_dist(S, V=10)

S = [[1,2,1,2,1], [2,2,2,2], [3,3,3,3,10]]

@time get_cond_dist(S, V=10, α=10.0)

d = Multinomial(10, [0.9,0.1])

@time rand(Multinomial(1, [0.9,0.1]))

@time rand(Categorical([0.1,0.9]))


text:
include("maximum_step_size.jl")

function merit_function_derivative(res::class_homogeneous_residuals, theta::class_theta)
		return (theta.mu - 1.0) * res.mu  + (theta.dual - 1.0) * (res.r_G_norm + res.r_D_norm) + (theta.primal - 1.0) * res.r_P_norm;
end

function is_pos(arr::Array{Float64,1})
		for el in arr
				if el <= 0
						return false
				end
		end

		return true
end

function sx(vars::class_variables)
		sx_vec = [s(vars) .* x(vars); kappa(vars) * tau(vars)]
		return minimum(sx_vec)/mean(sx_vec) > 1e-2
end

function safe_guard(value::Float64, vars::class_variables)
		if is_pos(x(vars)) && is_pos(s(vars)) && tau(vars) > 0 && kappa(vars) > 0
			if sx(vars)
					return value
			else
					#print("SX!!")
					return Inf
			end
		else
			return Inf
		end
end

function current_merit_function(res::class_homogeneous_residuals, vars::class_variables)
		return safe_guard(res.mu + res.r_norm, vars)
end

function proximal_merit_function!(res::class_homogeneous_residuals, nlp_eval::internal_AbstractNLPEvaluator, vars::class_variables, newt::abstract_newton_direction, intial_point::class_variables)
		update_residuals!(res, nlp_eval, vars, newt, intial_point)
		return current_merit_function(res, vars)
end

function objective_merit_function(vars::class_variables, nlp::internal_AbstractNLPEvaluator, τ_k::Float64, μ_k::Float64)
		# (τ_k)^2 * c(𝑥/τ) - μ * ( sum(log(𝑥)) + log(τ) )
		obj = internal_eval_c(nlp, x_scaled(vars))
		value = (τ_k)^2 * obj - μ_k * ( sum(log(x(vars))) + log(tau(vars)) )
		return value
end

#function derivative_merit_function

function weighted_objective_constraints_merit_function(vars::class_variables, intial_point::class_variables, nlp::internal_AbstractNLPEvaluator, τ_k::Float64, μ_k::Float64)
		 y_mod = 1e-8 * (y(vars) - y(intial_point))
		 primal = norm(tau(vars) * internal_eval_a(nlp, x_scaled(vars)) + y_mod,1)
		 obj = objective_merit_function(vars, nlp, τ_k, μ_k)
		 #@show primal, obj
		 value = primal
		 return safe_guard(value,vars)
end

function step(newton_solver::abstract_newton_direction, nlp::internal_AbstractNLPEvaluator, intial_point::class_variables, direction::class_variables, theta::class_theta)
	try
		alpha_max = maximum_step(intial_point, direction);
		alpha = 0.5 * min(1.0, 0.8 * alpha_max);
		new_vars = deepcopy(intial_point)
		move!(new_vars, alpha, direction);

		res = newton_solver.residuals
		orginal_merit_func = proximal_merit_function!(res, nlp, intial_point, newton_solver, intial_point)
		new_merit_func = proximal_merit_function!(deepcopy(res), nlp, new_vars, newton_solver, intial_point)
		v = new_merit_func/orginal_merit_func

		return new_vars, alpha, true, v
	catch e
		println("ERROR in step")
		throw(e)
	end
end

function line_search(newton_solver::abstract_newton_direction, nlp::internal_AbstractNLPEvaluator, intial_point::class_variables, direction::class_variables, theta::class_theta)
	try
		alpha_max = maximum_step(intial_point, direction);
		alpha = min(1.0, 0.9 * alpha_max);

		res = newton_solver.residuals

		merit_function_diff = merit_function_derivative(res, theta)

		if merit_function_diff >= 0.0
				println(merit_function_diff)
				error("non-decreasing direction")
		end

		#finite_difference_merit_function(res, theta, nlp, newton_solver, intial_point, direction)

		#orginal_merit_func = current_merit_function(res)
		orginal_merit_func = proximal_merit_function!(res, nlp, intial_point, newton_solver, intial_point)

		actual_gain = 0.0
		new_merit_func = 0.0
		new_vars = intial_point
		for i = 1:20
				new_vars = deepcopy(intial_point)
				move!(new_vars, alpha, direction);

				# why do I need to deep copy the residuals ???????
				new_res = deepcopy(res)
				new_merit_func = proximal_merit_function!(new_res, nlp, new_vars, newton_solver, intial_point)

				# try non-linear update
				#
				nonlinear_vars = deepcopy(new_vars)
				non_linear_update(new_res, res, nonlinear_vars, (1 - alpha * (1 - theta.mu)))
				nonlinear_merit_func = proximal_merit_function!(new_res, nlp, nonlinear_vars, newton_solver, intial_point)

				is_non_linear_update = false
				if nonlinear_merit_func < new_merit_func
						new_vars = nonlinear_vars
						is_non_linear_update = true
						new_merit_func = nonlinear_merit_func
				end
				#

				# what is the gain?
				actual_gain = new_merit_func - orginal_merit_func
				v = new_merit_func/orginal_merit_func

				if actual_gain <= 1e-4 * merit_function_diff * alpha || alpha < 1e-4
					return new_vars, alpha, true, v
				end

				alpha = 0.5 * alpha
		end

		@show alpha
		@show actual_gain, merit_function_diff * alpha
		@show orginal_merit_func, new_merit_func
		@show norm(direction._v)
		@show norm(intial_point._v)
		@show norm(new_vars._v)
		@show norm(y(direction))
		@show norm(s(direction))
		@show tau(direction), kappa(direction)

		error("failure of line search")
	catch e
		println("ERROR in line_search")
		throw(e)
	end
end

function exact_line_search(newton_solver::abstract_newton_direction, nlp::internal_AbstractNLPEvaluator, intial_point::class_variables, direction::class_variables, theta::class_theta)
	try
		alpha_max = maximum_step(intial_point, direction);
		alpha = min(1.0, 0.9 * alpha_max);

		res = newton_solver.residuals

		merit_function_diff = merit_function_derivative(res, theta)

		if merit_function_diff >= 0.0
				error("non-decreasing direction")
		end

		#finite_difference_merit_function(res, theta, nlp, newton_solver, intial_point, direction)

		#orginal_merit_func = current_merit_function(res)
		orginal_merit_func = proximal_merit_function!(res, nlp, intial_point, newton_solver, intial_point)

		actual_gain = 0.0

		best_gain = actual_gain
		best_point = intial_point;
		best_alpha = 0.0
		best_v = NaN
		for i = 1:30
				new_vars = deepcopy(intial_point)
				move!(new_vars, alpha, direction);

				# why do I need to deep copy the residuals ???????
				new_res = deepcopy(res)
				new_merit_func = proximal_merit_function!(new_res, nlp, new_vars, newton_solver, intial_point)

				# try non-linear update
				#
					nonlinear_vars = deepcopy(new_vars)
					non_linear_update(new_res, res, nonlinear_vars, (1 - alpha * (1 - theta.mu)))
					nonlinear_merit_func = proximal_merit_function!(new_res, nlp, nonlinear_vars, newton_solver, intial_point)

					is_non_linear_update = false
					if nonlinear_merit_func < new_merit_func
							new_vars = nonlinear_vars
							is_non_linear_update = true
							new_merit_func = nonlinear_merit_func
					end
				#

				# what is the gain?
				actual_gain = new_merit_func - orginal_merit_func
				v = new_merit_func/orginal_merit_func

				if actual_gain <= best_gain
					best_gain = actual_gain
					best_point = new_vars;
					best_alpha = alpha
					best_v = v
				end

				alpha = 0.6 * alpha
		end
		if best_gain >= 0.0
				@show alpha
				@show actual_gain, merit_function_diff * alpha
				error("failure of line search")
		end

		return best_point, best_alpha, true, best_v
	catch e
		println("ERROR in line_search")
		throw(e)
	end
end

####################
#	TO DO!!!!
####################
function weighted_objective_constraint_line_search(newton_solver::abstract_newton_direction, nlp::internal_AbstractNLPEvaluator, intial_point::class_variables, direction::class_variables, theta::class_theta)
	try
		alpha_max = maximum_step(intial_point, direction);
		alpha = min(1.0, 0.9 * alpha_max);

		res = newton_solver.residuals

		merit_function_diff = -(1.0 - theta.primal) * res.r_P_norm

		if merit_function_diff >= 0.0
				@show merit_function_diff
				error("non-decreasing direction")
		end

		#finite_difference_merit_function(res, theta, nlp, newton_solver, intial_point, direction)

		τ_k = tau(intial_point)
		μ_k = theta.mu * mu(intial_point)
		orginal_merit_func = weighted_objective_constraints_merit_function(intial_point, intial_point, nlp, τ_k, μ_k)

		new_vars = deepcopy(intial_point)
		move!(new_vars, 1e-8, direction);

		inc_merit_func = weighted_objective_constraints_merit_function(new_vars, intial_point, nlp, τ_k, μ_k)

		@show merit_function_diff, (inc_merit_func - orginal_merit_func)/1e-8
		@show -μ_k * (sum(x(direction) ./ x(intial_point)) + tau(direction)/tau(intial_point))# + sum(1 ./ s(direction)) + 1./kappa(direction))
		@show dot(y(direction), res.r_P)
		#@show (1 - theta.mu) * dot(-y(direction), res.r_P) + dot(x(direction), res.r_D) + tau(direction) * res.r_G
		actual_gain = 1.0

		best_gain = actual_gain
		best_point = intial_point;
		best_alpha = 0.0
		best_v = NaN
		for i = 1:4
				new_vars = deepcopy(intial_point)
				move!(new_vars, alpha, direction);

				# why do I need to deep copy the residuals ???????
				new_res = deepcopy(res)
				new_merit_func = weighted_objective_constraints_merit_function(new_vars, intial_point, nlp, τ_k, μ_k)

				# what is the gain?
				actual_gain = new_merit_func - orginal_merit_func
				v = new_merit_func/orginal_merit_func

				if actual_gain <= best_gain
					best_gain = actual_gain
					best_point = new_vars;
					best_alpha = alpha
					best_v = v
				end

				alpha = 0.6 * alpha
		end
		if best_gain >= 1.0
				@show alpha
				@show actual_gain #, merit_function_diff * alpha
				@show norm(x(direction))
				@show norm(y(direction))
				@show norm(s(direction))
				@show tau(direction), kappa(direction)
				error("failure of weighted line search")
		end

		return best_point, best_alpha, true, best_v
	catch e
		println("ERROR in weighted_objective_constraint_line_search")
		throw(e)
	end
end

function oli_line_search(newton_solver::abstract_newton_direction, nlp::internal_AbstractNLPEvaluator, intial_point::class_variables, direction::class_variables, theta::class_theta)
	try
		# make sufficient progress on constraints
	catch

	end
end

#non_linear_update(previous_res, newton_solver.nlp_vals, deepcopy(new_vars), 0.99)
function non_linear_update(new_res::class_homogeneous_residuals, previous_res::class_homogeneous_residuals, point::class_variables, reduction_factor::Float64)
		s( point, reduction_factor * previous_res.r_D + (s(point) - new_res.r_D) )
		kappa( point, reduction_factor * previous_res.r_G + (kappa(point) - new_res.r_G) )
end

function s_heuristic()

end

#

max_forks_repo_name: ohinder/Everest


text:
using
  PyPlot,
  JLD2,
  Printf,
  FourierFlows

using FFTW: ifft

import GeophysicalFlows.BarotropicQG
import GeophysicalFlows.BarotropicQG: energy, meanenergy, enstrophy, meanenstrophy

# Numerical parameters and time-stepping parameters
nx  = 512      # 2D resolution = nx^2
stepper = "FilteredETDRK4"   # timestepper
dt  = 2e-2     # timestep
nsteps = 20000 # total number of time-steps
nsubs  = 500   # number of time-steps for plotting
               # (nsteps must be multiple of nsubs)

# Physical parameters
Lx  = 2π       # domain size
nu  = 8.0e-10  # viscosity
nnu = 2        # viscosity order
f0  = -1.0     # Coriolis parameter
beta = 1.4015  # the y-gradient of planetary PV
mu   = 1.0e-2  # linear drag
   F = 0.0012  # normalized wind stress forcing on domain-averaged
               # zonal flow U(t) flow

# Topographic PV
topoPV(x, y) = @. 2*cos(10x)*cos(10y)

# Forcing on the domain-averaged U equation
calcFU(t) = F

# Initialize problem
prob = BarotropicQG.ForcedProblem(nx=nx, Lx=Lx, f0=f0, beta=beta, eta=topoPV,
                  calcFU=calcFU, nu=nu, nnu=nnu, mu=mu, dt=dt, stepper=stepper)
sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid

x, y = gridpoints(g)

# Files
filepath = "."
plotpath = "./plots_acctopo"
plotname = "snapshots"
filename = joinpath(filepath, "acctopo.jl.jld2")

# File management
if isfile(filename); rm(filename); end
if !isdir(plotpath); mkdir(plotpath); end

# Initialize with zeros
BarotropicQG.set_zeta!(prob, 0*x)

# Create Diagnostics
E = Diagnostic(energy, prob; nsteps=nsteps)
Q = Diagnostic(enstrophy, prob; nsteps=nsteps)
Emean = Diagnostic(meanenergy, prob; nsteps=nsteps)
Qmean = Diagnostic(meanenergy, prob; nsteps=nsteps)
diags = [E, Emean, Q, Qmean]

# Create Output
get_sol(prob) = sol # extracts the Fourier-transformed solution
get_u(prob) = irfft(im*g.lr.*g.invKrsq.*sol, g.nx)
out = Output(prob, filename, (:sol, get_sol), (:u, get_u))

function plot_output(prob, fig, axs; drawcolorbar=false)

  # Plot the PV field and the evolution of energy and enstrophy.

  sol, v, p, g = prob.sol, prob.vars, prob.params, prob.grid
  BarotropicQG.updatevars!(prob)

  sca(axs[1])
  pcolormesh(x, y, v.q)
  axis("square")
  xlim(0, 2)
  xticks(0:0.5:2)
  ylim(0, 2)
  yticks(0:0.5:2)
  title(L"$\nabla^2\psi + \eta$ (part of the domain)")
  if drawcolorbar==true
    colorbar()
  end

  sca(axs[2])
  cla()
  plot(mu*E.t[1:E.i], E.data[1:E.i], label=L"$E_{\psi}$")
  plot(mu*E.t[1:Emean.i], Emean.data[1:Emean.i], label=L"$E_U$")

  xlabel(L"\mu t")
  ylabel(L"E")
  legend()

  sca(axs[3])
  cla()
  plot(mu*Q.t[1:Q.i], Q.data[1:Q.i], label=L"$Q_{\psi}$")
  plot(mu*Qmean.t[1:Qmean.i], Qmean.data[1:Qmean.i], label=L"$Q_U$")
  xlabel(L"\mu t")
  ylabel(L"Q")
  legend()
  pause(0.001)
end

fig, axs = subplots(ncols=3, nrows=1, figsize=(15, 4))
plot_output(prob, fig, axs; drawcolorbar=true)

# Step forward
startwalltime = time()

while cl.step < nsteps
  stepforward!(prob, diags, nsubs)

  # Message
  log = @sprintf("step: %04d, t: %d, E: %.4f, Q: %.4f, τ: %.2f min",
    cl.step, cl.t, E.data[E.i], Q.data[Q.i],
    (time()-startwalltime)/60)

  println(log)

  plot_output(prob, fig, axs; drawcolorbar=false)

end
println((time()-startwalltime))

plot_output(prob, fig, axs; drawcolorbar=false)

savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), cl.step)
savefig(savename, dpi=240)


text:
import Mercury as Hg
using ZMQ
using Sockets
using Test
using BenchmarkTools
if !isdefined(@__MODULE__, :TestMsg)
    include("jlout/test_msg_pb.jl")
end

mutable struct PubNode <: Hg.Node
    ctx::ZMQ.Context
    nodeio::Hg.NodeIO
    test_msg::TestMsg
    should_finish::Bool
end

function PubNode(ctx)
    test_msg = TestMsg(x = 10, y = 10, z = 10)
    nodedata = Hg.NodeIO()
    PubNode(ctx, nodedata, test_msg, false)
end

function Hg.setupIO!(node::PubNode, nodeio::Hg.NodeIO)
    ctx = node.ctx
    addr = ip"127.0.0.1"
    port = 5555
    pub = Hg.ZmqPublisher(ctx, addr, port, name = "test_pub")
    empty!(nodeio.pubs)
    Hg.add_publisher!(nodeio, node.test_msg, pub)
    return nothing
end

function Hg.compute(node::PubNode)
    node.test_msg.x += 1
    node.test_msg.y += 2
    node.test_msg.z += 3
    Hg.publish(Hg.getIO(node).pubs[1])
    if (node.test_msg.x % 100) == 0
        println("Sent x value of ", node.test_msg.x)
    end
end
Hg.getrate(::PubNode) = 10

mutable struct SubNode <: Hg.Node
    ctx::ZMQ.Context
    nodeio::Hg.NodeIO
    test_msg::TestMsg
    should_finish::Bool
end

function SubNode(ctx = ZMQ.Context())
    test_msg = TestMsg(x = 0, y = 0, z = 0)
    nodedata = Hg.NodeIO()
    SubNode(ctx, nodedata, test_msg, false)
end

function Hg.setupIO!(node::SubNode, nodeio::Hg.NodeIO)
    ctx = node.ctx
    addr = ip"127.0.0.1"
    port = 5555
    sub = Hg.ZmqSubscriber(ctx, addr, port)
    empty!(nodeio.subs)
    Hg.add_subscriber!(nodeio, node.test_msg, sub)
end

Hg.getrate(::SubNode) = 10

function Hg.compute(node::SubNode)
    x = node.test_msg.x
    y = node.test_msg.y
    z = node.test_msg.z
    submsg = Hg.getIO(node).subs[1]

    x = node.test_msg.x
    y = node.test_msg.y
    z = node.test_msg.z

    println("Received x = $x")
end

@testset "Node" begin

    println("######## NODE TESTS #############")
    ## Initialize publisher
    ctx = ZMQ.Context()
    node = PubNode(ctx)
    Hg.setupIO!(node, Hg.getIO(node))
    @test Hg.getname(node) == "PubNode"
    @test Hg.getname(Hg.getpublisher(node, 1)) == "test_pub"
    @test Hg.getname(Hg.getpublisher(node, 1).pub) == "test_pub"
    @test Hg.getname(Hg.getpublisher(node, "test_pub")) == "test_pub"
    @test Hg.getpublisher(node, "test_pub").msg isa TestMsg
    @test Hg.numpublishers(node) == 1
    @test Hg.numsubscribers(node) == 0

    Hg.printstatus(node)

    ## Initialize subscriber
    using Base.Threads
    Hg.reset_sub_count()
    subnode = SubNode()
    Hg.setupIO!(subnode, Hg.getIO(subnode))
    sub = Hg.getsubscriber(subnode, 1)
    @test isopen(sub.sub)
    @test sub === Hg.getsubscriber(subnode, "subscriber_1")
    @test isnothing(Hg.getsubscriber(subnode, "sub2"))
    @test Hg.numsubscribers(subnode) == 1
    @test Hg.numpublishers(subnode) == 0

    Hg.printstatus(subnode)

    ## Launch tasks
    task = @async Hg.launch(node)
    @test !istaskdone(task)
    @test !istaskfailed(task)
    @test Hg.node_sockets_are_open(node)

    Hg.getflags(subnode).should_finish[] = false
    subtask = Threads.@spawn Hg.launch(subnode)
    @test !istaskdone(subtask)
    @test !istaskfailed(subtask)
    @test Hg.node_sockets_are_open(subnode)

    ##
    sleep(1)
    Hg.stopnode(node)     # closing publisher first should be fine
    Hg.stopnode(subnode)
    sleep(0.1)  # give some time for the tasks to finish
    @test !Hg.getflags(node).is_running[]
    @test !Hg.getflags(subnode).is_running[]

    @test !Hg.getflags(node).did_error[]
    @test !Hg.getflags(subnode).did_error[]
    Hg.printstatus(node)
    Hg.printstatus(subnode)

    sleep(0.5)
    pub = Hg.getpublisher(node, 1)
    @test !isopen(pub.pub)
    sub = Hg.getsubscriber(subnode, 1)
    @test !isopen(sub.sub)

    dx = node.test_msg.x - subnode.test_msg.x
    dy = node.test_msg.y - subnode.test_msg.y
    dz = node.test_msg.z - subnode.test_msg.z
    @test abs(dx) <= 1
    @test abs(dy) <= 2
    @test abs(dz) <= 3
    @show node.test_msg.x
end


text:
"""
    coordGeoRegion(geo::PolyRegion) ->
        blon::Vector{<:Real}, blat::Vector{<:Real},
        slon::Vector{<:Real}, slat::Vector{<:Real},

For a given RectRegion, extract the [N,S,E,W] bounds and create a longitude and latitude vectors for the bound and the shape of the GeoRegion

Arguments
=========

- `geo` : A PolyRegion (i.e. a polygonal GeoRegion)

Returns
=======

- `blon` : A vector of longitude points for the bound of the GeoRegion
- `blat` : A vector of latitude points for the bound of the GeoRegion
- `slon` : A vector of longitude points for the shape of the GeoRegion
- `slat` : A vector of latitude points for the shape of the GeoRegion
"""
function coordGeoRegion(GeoReg::PolyRegion)

    shape = GeoReg.shape
    npnt  = length(shape)
    lon   = zeros(21,npnt-1)
    lat   = zeros(21,npnt-1)

    for ipnt = 1 : (npnt-1)
        lon[:,ipnt] .= collect(range(shape[ipnt][1],shape[ipnt+1][1],21))
        lat[:,ipnt] .= collect(range(shape[ipnt][2],shape[ipnt+1][2],21))
    end

    lon = lon[:]
    lat = lat[:]

    N = GeoReg.N
    S = GeoReg.S
    E = GeoReg.E
    W = GeoReg.W

    blon = vcat(range(W,E,21),range(E,E,21),range(E,W,21),range(W,W,21))
    blat = vcat(range(N,N,21),range(N,S,21),range(S,S,21),range(S,N,21))

    return blon,blat,lon,lat

end

"""
    coordGeoRegion(geo::RectRegion) -> blon::Vector{<:Real}, blat::Vector{<:Real}

For a given RectRegion, extract the [N,S,E,W] bounds and create a longitude and latitude vector.

Arguments
=========

- `geo` : A RectRegion (i.e. a rectilinear GeoRegion)

Returns
=======

- `blon` : A vector of longitude points for the bound of the GeoRegion
- `blat` : A vector of latitude points for the bound of the GeoRegion
"""
function coordGeoRegion(GeoReg::RectRegion)

    N = GeoReg.N
    S = GeoReg.S
    E = GeoReg.E
    W = GeoReg.W

    blon = vcat(range(W,E,21),range(E,E,21),range(E,W,21),range(W,W,21))
    blat = vcat(range(N,N,21),range(N,S,21),range(S,S,21),range(S,N,21))

    return blon,blat

end


text:
function _extract_dense(model, fields)
    extracted = []
    for f in fields
        dat = getfield(model, f)
        if dat isa AbstractArray{<:Number} && issparse(dat)
            dat = Array(dat)
        end
        push!(extracted, dat)
    end
    return extracted
end


text:
function Base.copy(Dadj::Adjoint{T,<:DArray{T,2}}) where T
    D = parent(Dadj)
    DArray(reverse(size(D)), procs(D)) do I
        lp = Array{T}(undef, map(length, I))
        rp = convert(Array, D[reverse(I)...])
        adjoint!(lp, rp)
    end
end

function Base.copy(Dtr::Transpose{T,<:DArray{T,2}}) where T
    D = parent(Dtr)
    DArray(reverse(size(D)), procs(D)) do I
        lp = Array{T}(undef, map(length, I))
        rp = convert(Array, D[reverse(I)...])
        transpose!(lp, rp)
    end
end

const DVector{T,A} = DArray{T,1,A}
const DMatrix{T,A} = DArray{T,2,A}

# Level 1

function axpy!(α, x::DArray, y::DArray)
    if length(x) != length(y)
        throw(DimensionMismatch("vectors must have same length"))
    end
    asyncmap(procs(y)) do p
        @async remotecall_fetch(p) do
            axpy!(α, localpart(x), localpart(y))
            return nothing
        end
    end
    return y
end

function dot(x::DVector, y::DVector)
    if length(x) != length(y)
        throw(DimensionMismatch(""))
    end

    results=Any[]
    asyncmap(procs(x)) do p
            push!(results, remotecall_fetch((x, y) -> dot(localpart(x), makelocal(y, localindices(x)...)), p, x, y))
    end
    return reduce(+, results)
end

function norm(x::DArray, p::Real = 2)
    results = []
    @sync begin
        for pp in procs(x)
            @async push!(results, remotecall_fetch(() -> norm(localpart(x), p), pp))
        end
    end
    return norm(results, p)
end

function LinearAlgebra.rmul!(A::DArray, x::Number)
    @sync for p in procs(A)
        @async remotecall_fetch((A,x)->(rmul!(localpart(A), x); nothing), p, A, x)
    end
    return A
end

# Level 2
function add!(dest, src, scale = one(dest[1]))
    if length(dest) != length(src)
        throw(DimensionMismatch("source and destination arrays must have same number of elements"))
    end
    if scale == one(scale)
        @simd for i = eachindex(dest)
            @inbounds dest[i] += src[i]
        end
    else
        @simd for i = eachindex(dest)
            @inbounds dest[i] += scale*src[i]
        end
    end
    return dest
end

function LinearAlgebra.mul!(y::DVector, A::DMatrix, x::AbstractVector, α::Number = 1, β::Number = 0)

    # error checks
    if size(A, 2) != length(x)
        throw(DimensionMismatch(""))
    end
    if y.cuts[1] != A.cuts[1]
        throw(ArgumentError("cuts of output vector must match cuts of first dimension of matrix"))
    end

    # Multiply on each tile of A
    R = Array{Future}(undef, size(A.pids))
    for j = 1:size(A.pids, 2)
        xj = x[A.cuts[2][j]:A.cuts[2][j + 1] - 1]
        for i = 1:size(A.pids, 1)
            R[i,j] = remotecall(procs(A)[i,j]) do
                localpart(A)*convert(localtype(x), xj)
            end
        end
    end

    # Scale y if necessary
    if β != one(β)
        asyncmap(procs(y)) do p
            remotecall_fetch(p) do
                if !iszero(β)
                    rmul!(localpart(y), β)
                else
                    fill!(localpart(y), 0)
                end
                return nothing
            end
        end
    end

    # Update y
    @sync for i = 1:size(R, 1)
        p = y.pids[i]
        for j = 1:size(R, 2)
            rij = R[i,j]
            @async remotecall_fetch(() -> (add!(localpart(y), fetch(rij), α); nothing), p)
        end
    end

    return y
end

function LinearAlgebra.mul!(y::DVector, adjA::Adjoint{<:Number,<:DMatrix}, x::AbstractVector, α::Number = 1, β::Number = 0)

    A = parent(adjA)

    # error checks
    if size(A, 1) != length(x)
        throw(DimensionMismatch(""))
    end
    if y.cuts[1] != A.cuts[2]
        throw(ArgumentError("cuts of output vector must match cuts of second dimension of matrix"))
    end

    # Multiply on each tile of A
    R = Array{Future}(undef, reverse(size(A.pids)))
    for j = 1:size(A.pids, 1)
        xj = x[A.cuts[1][j]:A.cuts[1][j + 1] - 1]
        for i = 1:size(A.pids, 2)
            R[i,j] = remotecall(() -> localpart(A)'*convert(localtype(x), xj), procs(A)[j,i])
        end
    end

    # Scale y if necessary
    if β != one(β)
        asyncmap(procs(y)) do p
            remotecall_fetch(p) do
                if !iszero(β)
                    rmul!(localpart(y), β)
                else
                    fill!(localpart(y), 0)
                end
                return nothing
            end
        end
    end

    # Update y
    @sync for i = 1:size(R, 1)
        p = y.pids[i]
        for j = 1:size(R, 2)
            rij = R[i,j]
            @async remotecall_fetch(() -> (add!(localpart(y), fetch(rij), α); nothing), p)
        end
    end
    return y
end

function LinearAlgebra.lmul!(D::Diagonal, DA::DMatrix)
    d = D.diag
    s = verified_destination_serializer(procs(DA), size(DA.indices)) do pididx
        d[DA.indices[pididx][1]]
    end
    map_localparts!(DA) do lDA
        lmul!(Diagonal(localpart(s)), lDA)
    end
end

function LinearAlgebra.rmul!(DA::DMatrix, D::Diagonal)
    d = D.diag
    s = verified_destination_serializer(procs(DA), size(DA.indices)) do pididx
        d[DA.indices[pididx][2]]
    end
    map_localparts!(DA) do lDA
        rmul!(lDA, Diagonal(localpart(s)))
    end
end

# Level 3
function _matmatmul!(C::DMatrix, A::DMatrix, B::AbstractMatrix, α::Number, β::Number, tA)
    # error checks
    Ad1, Ad2 = (tA == 'N') ? (1,2) : (2,1)
    mA, nA = (size(A, Ad1), size(A, Ad2))
    mB, nB = size(B)
    if mB != nA
        throw(DimensionMismatch("matrix A has dimensions ($mA, $nA), matrix B has dimensions ($mB, $nB)"))
    end
    if size(C,1) != mA || size(C,2) != nB
        throw(DimensionMismatch("result C has dimensions $(size(C)), needs ($mA, $nB)"))
    end
    if C.cuts[1] != A.cuts[Ad1]
        throw(ArgumentError("cuts of the first dimension of the output matrix must match cuts of dimension $Ad1 of the first input matrix"))
    end

    # Multiply on each tile of A
    if tA == 'N'
        R = Array{Future}(undef, size(procs(A))..., size(procs(C), 2))
    else
        R = Array{Future}(undef, reverse(size(procs(A)))..., size(procs(C), 2))
    end
    for j = 1:size(A.pids, Ad2)
        for k = 1:size(C.pids, 2)
            Acuts = A.cuts[Ad2]
            Ccuts = C.cuts[2]
            Bjk = B[Acuts[j]:Acuts[j + 1] - 1, Ccuts[k]:Ccuts[k + 1] - 1]
            for i = 1:size(A.pids, Ad1)
                p = (tA == 'N') ? procs(A)[i,j] : procs(A)[j,i]
                R[i,j,k] = remotecall(p) do
                    if tA == 'T'
                        return transpose(localpart(A))*convert(localtype(B), Bjk)
                    elseif tA == 'C'
                        return adjoint(localpart(A))*convert(localtype(B), Bjk)
                    else
                        return localpart(A)*convert(localtype(B), Bjk)
                    end
                end
            end
        end
    end

    # Scale C if necessary
    if β != one(β)
        @sync for p in C.pids
            if β != zero(β)
                @async remotecall_fetch(() -> (rmul!(localpart(C), β); nothing), p)
            else
                @async remotecall_fetch(() -> (fill!(localpart(C), 0); nothing), p)
            end
        end
    end

    # Update C
    @sync for i = 1:size(R, 1)
        for k = 1:size(C.pids, 2)
            p = C.pids[i,k]
            for j = 1:size(R, 2)
                rijk = R[i,j,k]
                @async remotecall_fetch(d -> (add!(localpart(d), fetch(rijk), α); nothing), p, C)
            end
        end
    end
    return C
end

LinearAlgebra.mul!(C::DMatrix, A::DMatrix, B::AbstractMatrix, α::Number = 1, β::Number = 0) = _matmatmul!(C, A, B, α, β, 'N')
LinearAlgebra.mul!(C::DMatrix, A::Adjoint{<:Number,<:DMatrix}, B::AbstractMatrix, α::Number = 1, β::Number = 0) = _matmatmul!(C, parent(A), B, α, β, 'C')
LinearAlgebra.mul!(C::DMatrix, A::Transpose{<:Number,<:DMatrix}, B::AbstractMatrix, α::Number = 1, β::Number = 0) = _matmatmul!(C, parent(A), B, α, β, 'T')

_matmul_op = (t,s) -> t*s + t*s

function Base.:*(A::DMatrix, x::AbstractVector)
    T = Base.promote_op(_matmul_op, eltype(A), eltype(x))
    y = DArray(I -> Array{T}(undef, map(length, I)), (size(A, 1),), procs(A)[:,1], (size(procs(A), 1),))
    return mul!(y, A, x)
end
function Base.:*(A::DMatrix, B::AbstractMatrix)
    T = Base.promote_op(_matmul_op, eltype(A), eltype(B))
    C = DArray(I -> Array{T}(undef, map(length, I)),
            (size(A, 1), size(B, 2)),
            procs(A)[:,1:min(size(procs(A), 2), size(procs(B), 2))],
            (size(procs(A), 1), min(size(procs(A), 2), size(procs(B), 2))))
    return mul!(C, A, B)
end

function Base.:*(adjA::Adjoint{<:Any,<:DMatrix}, x::AbstractVector)
    A = parent(adjA)
    T = Base.promote_op(_matmul_op, eltype(A), eltype(x))
    y = DArray(I -> Array{T}(undef, map(length, I)),
            (size(A, 2),),
            procs(A)[1,:],
            (size(procs(A), 2),))
    return mul!(y, adjA, x)
end
function Base.:*(adjA::Adjoint{<:Any,<:DMatrix}, B::AbstractMatrix)
    A = parent(adjA)
    T = Base.promote_op(_matmul_op, eltype(A), eltype(B))
    C = DArray(I -> Array{T}(undef, map(length, I)), (size(A, 2),
        size(B, 2)),
        procs(A)[1:min(size(procs(A), 1), size(procs(B), 2)),:],
        (size(procs(A), 2), min(size(procs(A), 1), size(procs(B), 2))))
    return mul!(C, adjA, B)
end

function Base.:*(trA::Transpose{<:Any,<:DMatrix}, x::AbstractVector)
    A = parent(trA)
    T = Base.promote_op(_matmul_op, eltype(A), eltype(x))
    y = DArray(I -> Array{T}(undef, map(length, I)),
            (size(A, 2),),
            procs(A)[1,:],
            (size(procs(A), 2),))
    return mul!(y, trA, x)
end
function Base.:*(trA::Transpose{<:Any,<:DMatrix}, B::AbstractMatrix)
    A = parent(trA)
    T = Base.promote_op(_matmul_op, eltype(A), eltype(B))
    C = DArray(I -> Array{T}(undef, map(length, I)), (size(A, 2),
        size(B, 2)),
        procs(A)[1:min(size(procs(A), 1), size(procs(B), 2)),:],
        (size(procs(A), 2), min(size(procs(A), 1), size(procs(B), 2))))
    return mul!(C, trA, B)
end


text:
using FastGaussQuadrature: gausshermite
using SpecialFunctions: loggamma
using ChainRulesCore: ChainRulesCore
using IrrationalConstants: sqrt2, invsqrtπ

struct DefaultExpectationMethod end

struct AnalyticExpectation end

struct GaussHermiteExpectation
    xs::Vector{Float64}
    ws::Vector{Float64}
end
GaussHermiteExpectation(n::Integer) = GaussHermiteExpectation(gausshermite(n)...)

ChainRulesCore.@non_differentiable gausshermite(n)

struct MonteCarloExpectation
    n_samples::Int
end

default_expectation_method(_) = GaussHermiteExpectation(20)

"""
    expected_loglikelihood(
        quadrature,
        lik,
        q_f::AbstractVector{<:Normal},
        y::AbstractVector,
    )

This function computes the expected log likelihood:

```math
    ∫ q(f) log p(y | f) df
```
where `p(y | f)` is the process likelihood. This is described by `lik`, which should be a
callable that takes `f` as input and returns a Distribution over `y` that supports
`loglikelihood(lik(f), y)`.

`q(f)` is an approximation to the latent function values `f` given by:
```math
    q(f) = ∫ p(f | u) q(u) du
```
where `q(u)` is the variational distribution over inducing points (see
[`elbo`](@ref)). The marginal distributions of `q(f)` are given by `q_f`.

`quadrature` determines which method is used to calculate the expected log
likelihood - see [`elbo`](@ref) for more details.

# Extended help

`q(f)` is assumed to be an `MvNormal` distribution and `p(y | f)` is assumed to
have independent marginals such that only the marginals of `q(f)` are required.
"""
expected_loglikelihood(quadrature, lik, q_f, y)

"""
    expected_loglikelihood(::DefaultExpectationMethod, lik, q_f::AbstractVector{<:Normal}, y::AbstractVector)

The expected log likelihood, using the default quadrature method for the given likelihood.
(The default quadrature method is defined by `default_expectation_method(lik)`, and should
be the closed form solution if it exists, but otherwise defaults to Gauss-Hermite
quadrature.)
"""
function expected_loglikelihood(
    ::DefaultExpectationMethod, lik, q_f::AbstractVector{<:Normal}, y::AbstractVector
)
    quadrature = default_expectation_method(lik)
    return expected_loglikelihood(quadrature, lik, q_f, y)
end

function expected_loglikelihood(
    mc::MonteCarloExpectation, lik, q_f::AbstractVector{<:Normal}, y::AbstractVector
)
    # take `n_samples` reparameterised samples
    f_μ = mean.(q_f)
    fs = f_μ .+ std.(q_f) .* randn(eltype(f_μ), length(q_f), mc.n_samples)
    lls = loglikelihood.(lik.(fs), y)
    return sum(lls) / mc.n_samples
end

# Compute the expected_loglikelihood over a collection of observations and marginal distributions
function expected_loglikelihood(
    gh::GaussHermiteExpectation, lik, q_f::AbstractVector{<:Normal}, y::AbstractVector
)
    # Compute the expectation via Gauss-Hermite quadrature
    # using a reparameterisation by change of variable
    # (see e.g. en.wikipedia.org/wiki/Gauss%E2%80%93Hermite_quadrature)
    return sum(Broadcast.instantiate(
        Broadcast.broadcasted(y, q_f) do yᵢ, q_fᵢ  # Loop over every pair
            # of marginal distribution q(fᵢ) and observation yᵢ
            expected_loglikelihood(gh, lik, q_fᵢ, yᵢ)
        end,
    ))
end

# Compute the expected_loglikelihood for one observation and a marginal distributions
function expected_loglikelihood(gh::GaussHermiteExpectation, lik, q_f::Normal, y)
    μ = mean(q_f)
    σ̃ = sqrt2 * std(q_f)
    return invsqrtπ * sum(Broadcast.instantiate(
        Broadcast.broadcasted(gh.xs, gh.ws) do x, w # Loop over every
            # pair of Gauss-Hermite point x with weight w
            f = σ̃ * x + μ
            loglikelihood(lik(f), y) * w
        end,
    ))
end

function expected_loglikelihood(
    ::AnalyticExpectation, lik, q_f::AbstractVector{<:Normal}, y::AbstractVector
)
    return error(
        "No analytic solution exists for $(typeof(lik)). Use `DefaultExpectationMethod`, `GaussHermiteExpectation` or `MonteCarloExpectation` instead.",
    )
end


text:
module Parallax

import Gloria: onevent!, render!, update!
using Gloria: Gloria, AbstractObject, Event, Layer, Scene, Texture, Window, Resources, iskey

using Colors: @colorant_str

struct Controls <: AbstractObject end

mutable struct Object <: AbstractObject
    texture::Texture
    x::Float64
    y::Float64
    vx::Float64
    vy::Float64
    θ::Float64
    ω::Float64
end
Object(fname::String, args...) = Object(Texture(Resources(window), fname), args...)

function onevent!(::Controls, e::Event{:mousemotion})
    # Pan
    if Gloria.getmousestate().left
        for layer in object_layers
            layer.x += e.rel_x * layer.scale
            layer.y += e.rel_y * layer.scale
        end
    end
    # Zoom
    if Gloria.getmousestate().right
        for layer in object_layers
            layer.scale *= 1.001^e.rel_y
            layer.x = (layer.x - width/2)*1.001^e.rel_y + width/2
            layer.y = (layer.y - height/2)*1.001^e.rel_y + height/2
            # truck/pedestal (moving the scene)
            # layer.x += e.rel_x
            # layer.y += e.rel_y
        end
    end
end

function onevent!(::Controls, e::Event{:mousebutton_down})
    (e.button == 1 || e.button == 3) && Gloria.setrelativemousemode!(true)
end

function onevent!(::Controls, e::Event{:mousebutton_up})
    (e.button == 1 || e.button == 3) && Gloria.setrelativemousemode!(false)
end

function onevent!(::Controls, e::Event{:key_down})
    iskey(e, "escape") && Gloria.quit!(window, e)
end

function update!(obj::Object; t::Float64, dt::Float64)
    obj.θ += obj.ω*dt
    obj.x += obj.vx*dt
    obj.y += obj.vy*dt
end

function render!(layer::Layer, obj::Object, frame::Int, fps::Float64)
    render!(layer, obj.texture, obj.x, obj.y, obj.θ)
end

# Setup

function main(; keepalive=true)
    @eval begin
        # const width, height = 1920, 1080
        const width, height = 800, 600
        const controls_layer = Layer([Controls()], show=false)
        const object_layers = [Layer(Object[], width/2, height/2, scale=1.5^-n) for n in 10:-1:0]
        const scene = Scene(controls_layer, object_layers...)
        const window = Window("Parallax", width, height, scene, fullscreen=false)

        for layer in object_layers
            push!(layer, [
                Object(
                    abspath(@__DIR__, "..", "assets", "sample.svg"),
                    (rand() - 0.5)*0.5width, (rand() - 0.5)*0.5height,
                    (rand() - 0.5)*width, (rand() - 0.5)*height,
                    rand()*360, (rand() - 0.5)*1080
                ) for _ in 1:50
            ]...)
        end
    end

    Gloria.run!(window)
    keepalive && wait(window)
end

# # precompile
# const dir = abspath(@__DIR__, "..", "precompile")
# const blacklist_import = [:Parallax, :unknown]
# const fnames = collect(filter(x->occursin(r"^precompile_.*\.jl$", x), readdir(dir)))
# const names = (fname->Symbol(match(r"^precompile_(.*)\.jl$", fname)[1])).(fnames)
# for name in names
#     name in blacklist_import && continue
#     try
#         @eval import $name
#     catch e
#         @warn "Failed import of: $name ($e)"
#     end
# end
# for fname in fnames
#     try
#         include(joinpath(dir, fname))
#         _precompile_()
#         catch e
#         @warn "Failed additional precompilation of: $fname ($e)"
#     end
# end

end # module


text:
function test_error_handler(f, e::Exception, args...)
    try
        throw(e)
    catch err
        f(err, args...)
    end
end

@testset "Error handling" begin
    # checkbodyforerrors
    resp = GraphQLClient.GQLResponse{Any}(nothing, nothing)
    @test GraphQLClient.checkbodyforerrors(resp) == false
    resp = GraphQLClient.GQLResponse{Any}(GraphQLClient.GQLError[], nothing)
    @test GraphQLClient.checkbodyforerrors(resp) == false
    resp = GraphQLClient.GQLResponse{Any}([GraphQLClient.GQLError("message", nothing)], nothing)
    @test_throws GraphQLClient.GraphQLError GraphQLClient.checkbodyforerrors(resp)

    # handle_error
    @test_throws ArgumentError test_error_handler(GraphQLClient.handle_error, ArgumentError("msg"))
    @test_throws HTTP.StatusError test_error_handler(GraphQLClient.handle_error, HTTP.StatusError(404, HTTP.Response(404;request=HTTP.Request(), body="{}")))
    @test_throws GraphQLClient.GraphQLError test_error_handler(GraphQLClient.handle_error, HTTP.StatusError(400, HTTP.Response(400;request=HTTP.Request(), body="{}")))

    # handle_deserialisation_error
    @test_throws MethodError test_error_handler(GraphQLClient.handle_deserialisation_error, MethodError(""), "", "")
    # Argument error without "invalid JSON"
    @test_throws ArgumentError test_error_handler(GraphQLClient.handle_deserialisation_error, ArgumentError(""), "", "")
    # Argument error with "invalid JSON" but default type
    @test_throws ArgumentError test_error_handler(GraphQLClient.handle_deserialisation_error, ArgumentError(""), "", Any)
    # Actual deserialisation error
    resp = HTTP.Response(200;body="{\"data\": {\"query\": 1}}")
    @test_throws ArgumentError test_error_handler(
        GraphQLClient.handle_deserialisation_error,
        ArgumentError("invalid JSON at byte"),
        resp,
        Dict, # anything apart from Any
    )
    @test_logs (:warn, r"Deserialisation of GraphQL response failed, trying to access execution errors") match_mode=:any try
        test_error_handler(GraphQLClient.handle_deserialisation_error, ArgumentError("invalid JSON at byte"), resp, Dict)
    catch
    end
    @test_logs (:error, r"No errors in GraphQL response") match_mode=:any try
        test_error_handler(GraphQLClient.handle_deserialisation_error, ArgumentError("invalid JSON at byte"), resp, Dict)
    catch
    end

    # Actual error that resulted in a deserialisation error
    resp = HTTP.Response(400;body="{\"errors\": [{\"message\": \"I stopped deserialisation!\"}]}")
    @test_throws GraphQLClient.GraphQLError test_error_handler(
        GraphQLClient.handle_deserialisation_error,
        ArgumentError("invalid JSON at byte"),
        resp,
        Dict, # anything apart from Any
    )
    @test_logs (:warn, r"Deserialisation of GraphQL response failed, trying to access execution errors") match_mode=:any try
        test_error_handler(GraphQLClient.handle_deserialisation_error, ArgumentError("invalid JSON at byte"), resp, Dict)
    catch
    end
end

function local_server_success(port)
    @async HTTP.serve(HTTP.Sockets.localhost, port) do req
        execution_string = String(req.body)
        return HTTP.Response("""
            {
                "data": {
                    "queryName": {
                        "field": \"$execution_string\"
                    }
                }
            }
        """)
    end
end

function local_server_success_json(port)
    @async HTTP.serve(HTTP.Sockets.localhost, port) do req
        return HTTP.Response(JSON3.write(
            Dict(
                "data" => Dict(
                    "queryName" => Dict(
                        "field" => JSON3.read(req.body)
                    )
                )
            )
        ))
    end
end

function local_server_error(port)
    @async HTTP.serve(HTTP.Sockets.localhost, port) do req
        str = """
            {
                "data": {
                    "queryName": {
                        "field": null
                    }
                },
                "errors": [
                    {
                        "message": "Error Text"
                    }
                ]
            }
        """
        return HTTP.Response(str)
    end
end

@testset "execute" begin
    # Successful query
    port = 7999
    local_server_success(7999)
    client = Client("http://$(HTTP.Sockets.localhost):$port";introspect=false)

    execution_string = "execute this"
    response = GraphQLClient._execute(client.endpoint, execution_string, Dict())

    # Test execution string passed through correctly
    @test response.data["queryName"]["field"] == execution_string
    @inferred GraphQLClient._execute(client.endpoint, execution_string, Dict())

    # Test struct type inference
    struct S
        field::String
    end
    StructTypes.StructType(::Type{S}) = StructTypes.Struct()
    response = GraphQLClient._execute(client.endpoint, execution_string, Dict(), S)
    @test response.data["queryName"].field == execution_string
    @inferred GraphQLClient._execute(client.endpoint, execution_string, Dict(), S)

    # Test error in response
    port = 7996
    local_server_error(7996)
    client = Client("http://$(HTTP.Sockets.localhost):$port";introspect=false)
    execution_string = "execute this"
    response = GraphQLClient._execute(client.endpoint, execution_string, Dict())
    @test !isnothing(response.errors)
    @test !isempty(response.errors)
    @test response.errors[1].message == "Error Text"
    @test_throws GraphQLClient.GraphQLError GraphQLClient._execute(client.endpoint, execution_string, Dict(), throw_on_execution_error=true)
    @test_throws GraphQLClient.GraphQLError GraphQLClient._execute(client.endpoint, execution_string, Dict(), S) # deserialisation fails

    # Test error before response - wrong input to retries causes method error
    @test_throws MethodError GraphQLClient._execute(client.endpoint, execution_string, Dict(), retries=Dict())

    # Test execute methods
    port = 7990
    local_server_success_json(7990)
    struct S1; query::String; end
    struct S2; field::S1; end
    StructTypes.StructType(::Type{S1}) = StructTypes.Struct()
    StructTypes.StructType(::Type{S2}) = StructTypes.Struct()
    client = Client("http://$(HTTP.Sockets.localhost):$port";introspect=false)
    @inferred GraphQLClient.execute(client.endpoint, Dict("query" => "val"))
    @inferred GraphQLClient.execute(client.endpoint, Dict("query" => "val"), Dict(), S2)
    response = GraphQLClient.execute(client.endpoint, Dict("query" => "val"))
    @test response.data["queryName"]["field"]["query"] == "val"
    @inferred GraphQLClient.execute(client.endpoint, "val")
    @inferred GraphQLClient.execute(client.endpoint, "val", Dict(), S2)
    response = GraphQLClient.execute(client.endpoint, "val")
    @test response.data["queryName"]["field"]["query"] == "val"
    @inferred GraphQLClient.execute(client, Dict("query" => "val"))
    @inferred GraphQLClient.execute(client, Dict("query" => "val"), S2)
    response = GraphQLClient.execute(client, Dict("query" => "val"))
    @test response.data["queryName"]["field"]["query"] == "val"
    @inferred GraphQLClient.execute(client, "val")
    @inferred GraphQLClient.execute(client, "val", S2)
    response = GraphQLClient.execute(client, "val")
    @test response.data["queryName"]["field"]["query"] == "val"

    global_graphql_client(client)
    @inferred GraphQLClient.execute(Dict("query" => "val"))
    @inferred GraphQLClient.execute(Dict("query" => "val"), S2)
    response = GraphQLClient.execute(Dict("query" => "val"))
    @test response.data["queryName"]["field"]["query"] == "val"
    @inferred GraphQLClient.execute("val")
    @inferred GraphQLClient.execute("val", S2)
    response = GraphQLClient.execute("val")
    @test response.data["queryName"]["field"]["query"] == "val"
end


text:
abstract type NoDecision <: Exception end
Base.showerror(io::IO, nd::NoDecision) = print(io, """
    Planner failed to choose an action because the following exception was thrown:
    $nd

    To specify an action for this case, use the default_action solver parameter.
    """)

struct AllSamplesTerminal <: NoDecision
    belief
end
Base.showerror(io::IO, ast::AllSamplesTerminal) = print(io, """
    Planner failed to choose an action because all states sampled from the belief were terminal.

    To see the belief, catch this exception as ex and see ex.belief.

    To specify an action for this case, use the default_action solver parameter.
    """)


text:
@testset "Internal Methods" begin
    m = InfiniteModel()
    @infinite_parameter(m, x ~ Normal())
    @infinite_parameter(m, xi[1:2] ~ MvNormal(ones(2), [1. 0.; 0. 2.]))
    domain1 = UniDistributionDomain(Normal())
    domain2 = MultiDistributionDomain(MvNormal(ones(2), [1. 0.; 0. 2.]))
    domain3 = CollectionDomain([domain1, domain1])
    domain4 = BadDomain()
    domain5 = IntervalDomain(0, 1)
    domain6 = CollectionDomain([domain1, domain5])
    domain7 = IntervalDomain(0, Inf)
    domain8 = CollectionDomain([domain5, domain5])
    # test _expect_coeffs
    @testset "_expect_coeffs" begin
        @test IOMT._expect_coeffs([1, 2]) == [0.5, 0.5]
    end
    # test _default_pdf
    @testset "_default_pdf" begin
        @test IOMT._default_pdf(0, 1, 3) == 0.5
    end
    # test generate_expect_data
    @testset "generate_expect_data" begin
        @test_throws ErrorException IOMT.generate_expect_data(domain4, x, 2)
        @test_throws ErrorException IOMT.generate_expect_data(domain1, x, 2, pdf = sin)
        @test IOMT.generate_expect_data(domain1, x, 2) isa FunctionalDiscreteMeasureData
        @test_throws ErrorException IOMT.generate_expect_data(domain2, xi, 2, pdf = sin)
        @test IOMT.generate_expect_data(domain2, xi, 2) isa FunctionalDiscreteMeasureData
        @test_throws ErrorException IOMT.generate_expect_data(domain6, xi, 2)
        @test_throws ErrorException IOMT.generate_expect_data(domain3, xi, 2, pdf = sin)
        @test IOMT.generate_expect_data(domain3, xi, 2) isa FunctionalDiscreteMeasureData
        @test_throws ErrorException IOMT.generate_expect_data(domain7, x, 2)
        @test_throws ErrorException IOMT.generate_expect_data(domain5, x, 2, d= 2)
        @test IOMT.generate_expect_data(domain5, x, 2) isa FunctionalDiscreteMeasureData
        @test IOMT.generate_expect_data(domain5, x, 2, pdf = (s) -> 0.5) isa FunctionalDiscreteMeasureData
        @test IOMT.generate_expect_data(domain8, xi[1], 2) isa FunctionalDiscreteMeasureData
    end
end

@testset "Expect" begin
    # Setup the model
    m = InfiniteModel()
    @infinite_parameter(m, x ~ Normal())
    @infinite_parameter(m, xi[1:2] ~ MvNormal(ones(2), [1. 0.; 0. 2.]))
    @infinite_parameter(m, y in [0, 1])
    @variable(m, inf, Infinite(x, xi, y))

    # Test normal usage
    @test InfiniteOpt._index_type(expect(inf, x)) == MeasureIndex
    warn = "Cannot specify a nonzero `num_supports` for individual " *
           "dependent parameters."
    @test_logs (:warn, warn) expect(inf, xi[1], num_supports = 3)
    @test !InfiniteOpt._is_expect(InfiniteOpt._core_variable_object(expect(inf, y)).data)
    @test 𝔼(inf, x) isa GeneralVariableRef

    # Test with JuMP container
    @infinite_parameter(m, z[2:3] in [-1, 1])
    @test_throws ErrorException expect(z[2], z)
    @test expect(z[2], z[2]) isa GeneralVariableRef
end

@testset "Macro" begin
    # Setup the model
    m = InfiniteModel()
    @infinite_parameter(m, x ~ Normal())
    @infinite_parameter(m, xi[1:2] ~ MvNormal(ones(2), [1. 0.; 0. 2.]))
    @variable(m, inf, Infinite(x, xi))

    # Test normal usage
    @test InfiniteOpt._index_type(@expect(inf, x)) == MeasureIndex
    @test InfiniteOpt._index_type(@expect(inf, xi, num_supports = 20)) == MeasureIndex
    @test_macro_throws ErrorException @expect(inf, x, 1)
    @test_macro_throws ErrorException @𝔼(inf, x, 5)
    @test @𝔼(inf, x) isa GeneralVariableRef
end


text:
#!/usr/bin/env julia

using Pkg, SHA, BinaryBuilderBase, Pkg.BinaryPlatforms, Pkg.Artifacts

include("common.jl")

function find_old_shards(name::String)
    all_shards = BinaryBuilderBase.all_compiler_shards()
    max_version = maximum([s.version for s in all_shards if s.name == name])
    return [s for s in all_shards if s.name == name && s.version != max_version]
end

artifacts_toml = joinpath(dirname(dirname(pathof(BinaryBuilderBase))), "Artifacts.toml")
function unbind_shard(cs::CompilerShard)
    @info("Unbinding $(repr(cs))")
    unbind_artifact!(artifacts_toml, BinaryBuilderBase.artifact_name(cs))
end

unbind_shard.(find_old_shards("Rootfs"))
unbind_shard.(find_old_shards("PlatformSupport"))
unbind_shard.(find_old_shards("RustBase"))
unbind_shard.(find_old_shards("RustToolchain"))

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
    repSTRdiff("runparams.txt", [hapsize_1, hapsize_2, ... hapsize_N];
        keepall=false, showhaps=false, min_uss_score=10)

Run getSTRdiff for multiple user-specified haplotype sizes.
Retrieve the best estimates from all runs and output a final table.

See getSTRdiff for additional details on input format.

Options:

    keepall        Keep all intermediate files [false]

    showhaps       Show the offspring and de novo trio haplotypes [false]

    min_uss_score  Minimum value for the unique solution score.  For values
                   below this number, both chromosomes from each of the
                   grandparents will be tested, and the minimum difference among
                   probabilities for both male and both female grandparental
                   chromosomes will be shown parenthetically next to the uss
                   score (e.g. 3,(17)).  This allows the sex of the transmitter
                   to be determined in (rare) cases where there is high similarity
                   between the transmitting parent's chromosomes and no one single
                   chromosome solution produced a high uss score.  For example,
                   [0.91, 0.90, 0.40, 0.73] would return (17), indicating support for
                   de novo transmission by the first grandparent (male)
                   even though the single chromosome solution is low.
                   This value is shown only when the uss score for the best estimate
                   falls below the min_uss_score setting [default = 10].

----------------------------------------------------------
Example sequence of commands for automation:
----------------------------------------------------------

1.  Run check\\_vcf\\_sites to clean input sites.

    check\\_vcf\\_sites("infile.vcf.gz", "datafile.txt"; c=4)

2.  Run repSTRdiff at a single haplotype length:

    repSTRdiff("runparams.txt", [20000])

    Remove any STR sites that fail from the input file.

3.  Run repSTRdiff over a range of haplotype lengths:

    repSTRdiff("runparams.txt", [20000:10000:300000...])

Examine the results in runparams.txt.out.final.

-----------------------------------------------------------

"""
function repSTRdiff(infile::String, hapsizes::Array; keepall=false, showhaps=false, min_uss_score=10)

    sites = []
    n = splitext(infile)
    fn = n[1]

    global s_haps = showhaps

    open(infile, "r") do file   #get all loci to run
        for i in eachline(file)

            if occursin(r"^#", i)
                continue
            end

            s = split(i, '\t')
            l = s[4]
            push!(sites, l)

        end
    end

    for i in hapsizes
        getSTRdiff(infile, i, showhaps=s_haps)
    end

    for i in sites
        run(`bash -c "grep $i $fn.*.out |sed 's/^.*\.[0-9]*\.out://' > $fn.$i.all.out" `)
    end

    header = "#ped\tiid\tloc\ttransmitting_gp\ttransmitting_gp_sex\ttransmitting_gp_str\ttransmitting_gp_str_size\tparent_str\tparent_str_size\tstr_allele_change\tsize_change_bp\tphased_offspring_with_de_novo\tfour_gp_match_probabilities\tuniq_soln_score\thapsize\tstr_type"

    OUT = open("$fn.out.final", "w")

    println(OUT, header)

    for i in sites
        best = ""
        open("$fn.$i.all.out") do file

            uss = 0; off = 0;

            for i in eachline(file)

                l = []
                l = split(i, '\t')
                f = String(l[4])
                global ped = l[1]
                global iid = l[2]
                global loc = l[3]

                if f == "na"
                    continue
                end

                if length(l) ≥ 14

                    t = parse(Int, l[14])
                    o = parse(Int, l[12])
                    loc = l[3]

                    if (t + o) > (uss + off)
                        best = i
                        uss = t
                        off = o
                    end

                else

                    best = "$ped\t$iid\t$loc\toptimal solution not found"

                end

            end

            if length(best) < 1
                best = "$ped\t$iid\t$loc\toptimal solution not found"
            end

            function uss_check(best, buss, min_uss_score)
                if (buss ≤ min_uss_score) && (length(best) > 70) && (!occursin(r"\t-9\t", best))
                    b = split(best, '\t')
                    s = parse(Int64, b[14])
                    c = replace(b[13], "[" => "")
                    c = replace(c, "]" => "")
                    p = parse.(Float64, split(c, r", "))
                    d1 = abs(p[1] - p[3])
                    d2 = abs(p[1] - p[4])
                    d3 = abs(p[2] - p[3])
                    d4 = abs(p[2] - p[4])
#                    println("$d1  $d2  $d3  $d4")
                    uss2 = floor(Int64,(minimum([d1, d2, d3, d4]) * 100))
                    b[14] = join([ s, ",", "(", uss2, ")" ],"")
                    best2 = ( join(b, "\t") )
                else
                    best2 = best
                end
                return best2
            end

            if length(best) > 70             #Retest grandparental chromosomes together
                g = split(best, '\t')
                buss = parse(Int64, g[14])
                best = uss_check(best, buss, min_uss_score)
            end

            println(OUT, best)

        end
    end

    close(OUT)

    function clean_dir()
        rt1 = run(`bash -c "rm -f *.out"`);
        rt2 = run(`bash -c "rm -f *_fam* *_trio*"`);
        rt3 = run(`bash -c "rm -f log* error*"`);
        println()
    end

    if keepall == false
        clean_dir()
    elseif keepall == true
        println("INFO: intermediate files for STR sites and haplotype lengths retained.")
    else
        clean_dir()
    end

end

max_stars_repo_name: ScottWatkins/STRdiff

max_issues_repo_name: ScottWatkins/STRdiff

max_forks_repo_name: ScottWatkins/STRdiff


text:
using Test
using OMEinsum
using SimpleTensorNetworks

@testset "contract" begin
    for (iAs, iBs, iOut) in [
        [(1,2,3), (2,3,4), (1,4)],
        [(1,3,2), (2,3,4), (4,1)],
        [(), (2,4), (4,2)],
        [(2,4), (), (4,2)],
        [(2,4), (1,), (4,2,1)],
        [(2,4), (2,4), ()],
        [(2,4), (4,), (2,)],
        ]
        A = asarray(randn(rand(4:12, length(iAs))...))
        B = asarray(randn([(iBs[i] in iAs) ? size(A, indexin(iBs[i], collect(iAs))[]) : rand(4:12) for i=1:length(iBs)]...))
        out = tensorcontract(iAs, A, iBs, B, iOut)
        eincode = EinCode{(iAs, iBs), iOut}()
        eout = eincode(A, B)
        @test eout ≈ out
    end
end

@testset "tensor contract" begin
    A = zeros(Float64, 10, 32, 21);
    B = zeros(Float64, 32, 11, 5, 2, 41, 10);
    tA = LabeledTensor(A, [1,2,3])
    tB = LabeledTensor(B, [2,4,5,6,7,1])

    tOut = tA * tB
    @test tOut.array == ein"abc,bdefga->cdefg"(A, B)
    @test tOut.labels == [3,4,5,6,7]

    tnet1 = TensorNetwork([tA])
    @test tnet1 isa TensorNetwork
    tnet = TensorNetwork([tA, tB])

    tOut2, contracted_labels = contract_label!(tnet, 1)
    @test tnet.tensors[] ≈ tOut
    @test tnet.tensors[] ≈ tOut2
    @test contracted_labels == [1, 2]
end

@testset "mul_dim" begin
    A = zeros(Float64, 10, 32, 21);
    B = zeros(Float64, 32, 11);
    tA = LabeledTensor(A, [1,2,3])
    tB = LabeledTensor(B, [2,4])
    tA1 = SimpleTensorNetworks.mul_dim(tA, B; dim=2)
    tA2 = tA * tB
    @test tA1.array ≈ permutedims(tA2.array, (1,3,2))
    @test tA1.labels == [1, 2, 3]
    @test tA2.labels == [1, 3, 4]
end


text:
using Test
using LinearAlgebra
using SummationByPartsOperators

# check construction of interior part of upwind operators
@testset "Check interior parts" begin
  N = 21
  xmin = 0.
  xmax = Float64(N + 1)
  interior = 10:N-10

  for acc_order in 2:7
    Dp_bounded = derivative_operator(Mattsson2017(:plus   ), 1, acc_order, xmin, xmax, N)
    Dm_bounded = derivative_operator(Mattsson2017(:minus  ), 1, acc_order, xmin, xmax, N)
    Dc_bounded = derivative_operator(Mattsson2017(:central), 1, acc_order, xmin, xmax, N)
    for compact in (true, false)
      show(IOContext(devnull, :compact=>false), Dp_bounded)
      show(IOContext(devnull, :compact=>false), Dm_bounded)
      show(IOContext(devnull, :compact=>false), Dc_bounded)
    end
    M = mass_matrix(Dp_bounded)
    @test M == mass_matrix(Dm_bounded)
    @test M == mass_matrix(Dc_bounded)
    Dp_periodic = periodic_derivative_operator(1, acc_order, xmin, xmax, N-1, -(acc_order - 1) ÷ 2)
    Dm_periodic = periodic_derivative_operator(1, acc_order, xmin, xmax, N-1, -acc_order + (acc_order - 1) ÷ 2)
    Dp = Matrix(Dp_bounded)
    Dm = Matrix(Dm_bounded)
    @test Dp[interior,interior] ≈ Matrix(Dp_periodic)[interior,interior]
    @test Dm[interior,interior] ≈ Matrix(Dm_periodic)[interior,interior]
    res = M * Dp + Dm' * M
    res[1,1] += 1
    res[end,end] -= 1
    @test norm(res) < N * eps()
    x = grid(Dp_bounded)
    for D in (Dp_bounded, Dm_bounded, Dc_bounded)
      @test norm(D * x.^0) < N * eps()
      for k in 1:acc_order÷2
        @test D * x.^k ≈ k .* x.^(k-1)
      end
      for k in acc_order÷2+1:acc_order
        @test (D * x.^k)[interior] ≈ (k .* x.^(k-1))[interior]
      end
    end
    diss = M * (Dp - Dm)
    @test diss ≈ diss'
    @test maximum(eigvals(Symmetric(diss))) < N * eps()
  end
end