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

using Markdown
using InteractiveUtils

# ╔═╡ 61c31884-4e1a-4f71-aa47-a0573d3e341a
begin
	using VegaLite
	using DataFrames
	using PlutoUI
	using CSV
	using Downloads
end

# ╔═╡ f4588c88-48d4-473e-b0d7-daa9ba1f0a00
html"""

Data Visualization

Basic diagrams with VegaLite

&copy  Dr. Roland Schätzle

"""

# ╔═╡ d8a90b71-b2ce-4766-83a7-f134ebd91298
PlutoUI.TableOfContents(title = "VegaLite")

# ╔═╡ 8d71d310-4db1-4921-8c86-26221f2ceb5a
md"""
# Data: Countries, Population and GDP

Load and prepare basic data for use in the following diagrams. For more information see the Pluto notebook "Basic diagrams with Gadfly".
"""

# ╔═╡ 5932694a-acae-4a1c-a2bb-ba0ec538c467
begin
	Downloads.download("https://raw.githubusercontent.com/roland-KA/StatisticalPlotsWithJulia/main/data/countries.csv", "countries.csv")
	countries = CSV.read("countries.csv", DataFrame)
	dropmissing!(countries)
	countries.GDPperCapita = countries.GDP ./ countries.Pop2019

	# group and aggregate regions
	regions = groupby(select(
			select(countries, Not(:Country)), Not(:Subregion)), :Region)
	regions_cum = combine(regions,
		:Pop2018 => sum, :Pop2019 => sum, :PopChangeAbs => sum, :GDP => sum,
		renamecols = false)

	# group and aggregate subregions
	subregions = groupby(select(countries, Not(:Country)), :Subregion)
	subregions_cum = combine(subregions, :Region => first,
		:Pop2018 => sum, :Pop2019 => sum, :PopChangeAbs => sum, :GDP => sum,
		renamecols = false)
end

# ╔═╡ d1bbe6fa-1f02-4158-a91f-7f4b3d641256
regions_cum

# ╔═╡ f791a7b9-549c-4f2d-b493-8bd3da683f3b
md"""
# Bar Plots
"""

# ╔═╡ 157e4e4a-a7d6-4c98-8843-d837a2232f36
md"""
## Population by Region
"""

# ╔═╡ ed0b8989-8f4e-4890-909e-218761ba4ea2
md"""
A bar plot to compare the population of the different regions in 2019.

First, a simple version using default values for several aspects of the diagram.
"""

# ╔═╡ 17ea91b1-db35-43b1-8cf0-f91f070b3049
let
regions_cum |>
	@vlplot(
		width = 600, height = 300,
		:bar,
		x = :Region, y = :Pop2019, color = :Region
	)
end

# ╔═╡ 17ebb9cc-1080-41f6-a914-f73afbb8e957
md"""
The second version uses serveral formatting options. It has different labels (x-axis, y-axis, title) as well as a more readable number format on the y-axis.
"""

# ╔═╡ afe3969f-6ff9-4c12-88ce-62dc2c3b2b3e
let
regions_cum |>
	@vlplot(
		width = 600, height = 300,
		title = "Population by Region, 2019",
		:bar,
		x = {:Region, title = "Region", axis = {labelAngle = 0}},
		y = {:Pop2019, title = "Population [millions]"},
		color = :Region,
		config = {background = "ghostwhite"}
	)
end

# ╔═╡ 9f9924c4-a937-4dea-b045-ba8101c2971b
md"""
## Population by Subregion
"""

# ╔═╡ a2fa5c60-aaf1-4a12-bb17-d260c41b85a6
md"""
Next we have a look at the population of the subregions.
"""

# ╔═╡ 55b9aa0b-da52-4a8c-a7b8-4bf25da2ca2f
let
subregions_cum |>
	@vlplot(
		width = 600, height = 300,
		:bar,
		x = :Subregion,	y = :Pop2019, color = :Region
	)
end

# ╔═╡ ea4f0f45-7fb3-4faa-bdb9-5b4a1df04280
md"""
As there are quite a few subregions, a horizontal bar diagram might be more readable. Apart from that we adapt the labels.
"""

# ╔═╡ 94bfd9c5-2877-417e-a18c-4983f243fc85
let
subregions_cum |>
	@vlplot(
		title = "Population by Subregion, 2019",
		width = 600, height = 300,
		:bar,
		x = {:Pop2019, title = "Population [millions]"},
		y = {:Subregion, title = "Subregion"},
		color = :Region,
		config = {background = "ghostwhite"}
	)
end

# ╔═╡ 0d84e074-af24-4823-a355-57a3485edc36
md"""
It get's even more readable, if we sort the subregions by population size before rendering the diagram.
"""

# ╔═╡ 9f72a8e5-8777-40c0-a6af-75b1a2aa7115
let
subregions_cum |>
	@vlplot(
		title = "Population by Subregion, 2019",
		width = 600, height = 300,
		:bar,
		x = {:Pop2019, title = "Population [millions]"},
		y = {:Subregion, sort = "-x", title = "Subregion"},
		color = :Region,
		config = {background = "ghostwhite"}
	)
end

# ╔═╡ 371dad4d-963c-42c9-a785-c477e124bcb4
md"""
The sorting can be done within VegaLite. `sort = "-x"` means: sort by x-values in descending order.
"""

# ╔═╡ 305f5666-2880-43dd-ae0b-7778ca4be212
md"""
# Scatter Plots
"""

# ╔═╡ f7e54982-88f4-4210-9343-ed187f65fa33
md"""
In the next step we have a look at the population at country level in relation to the growth rate. For this purpose a scatter plot is good way to show this relation.
"""

# ╔═╡ b38a4a8a-b655-4756-b396-bd6b2f97579c
let
	countries |>
		@vlplot(
			width = 600, height = 300,
			:point,
			x = :Pop2019, y = :PopChangePct, color = :Region
		)
end

# ╔═╡ e0b07b43-4446-43f5-814a-0642ece326e7
md"""
The distribution of the data is quite skewed. So a logarithmic scale on the x-axis might give a better insight of the data. And again, we add some labels.
"""

# ╔═╡ 1dc24d9c-8d8d-479c-b9e1-79792ee2ab0c
let
	countries |>
		@vlplot(
			title = "Population vs. Growth Rate, 2019",
			width = 600, height = 300,
			:point,
			x = {:Pop2019, title = "Population [millions]",
				scale = {type = :log, base = 10}},
			y = {:PopChangePct, title = "Growth Rate [%]"},
			color = :Region,
			config = {background = "ghostwhite"}
		)
end

# ╔═╡ ec791080-6c52-4c80-b7ca-43b604d4ede1
md"""
# Histograms
"""

# ╔═╡ 389e95fc-d410-481e-a753-f36ec70fa931
md"""
## Distribution of GDP per Capita
"""

# ╔═╡ f2753a4d-7bcd-4866-b59c-5568385e88c3
let
	countries |>
		@vlplot(
			width = 600, height = 300,
			:bar,
			x = {:GDPperCapita, bin = true}, y = "count()"
		)
end

# ╔═╡ 5782433d-5ab6-4a9c-a6fd-04610edb1dd8
md"""
The labels on the x-axis are set automatically to match class boundaries. And a reasonable bin size has been chosen by default.
"""

# ╔═╡ e46e377d-a637-4ba8-a9ba-789d6e6415f3
md"""
And again with labels (and the same number of bins as in the Gadfly example).
"""

# ╔═╡ 56892be5-2c1c-40a8-b0cd-d958d71a8845
let
	countries |>
		@vlplot(
			title = "Distribution of GDP per Capita, 2019",
			width = 600, height = 300,
			:bar,
			x = {:GDPperCapita, bin = {maxbins = 20},
				title = "GDP per Capita [USD]"},
			y = {"count()", title = "Number of countries"},
			config = {background = "ghostwhite"}
		)
end

# ╔═╡ a6ce4860-949b-4e84-b884-6090b40c3b96
md"""
# Box Plots and Violin Plots
"""

# ╔═╡ 11e4f95f-3849-4260-aed4-e0ac29cad61f
md"""
In the next step we have a look at the distribution of GDP per Capita among the different regions.
"""

# ╔═╡ 8e15d66c-0e8f-4842-bc6c-65f0c49934b7
md"""
## Box Plot: GDP per Capita by Region, 2019
"""

# ╔═╡ 8580f7b7-564d-4469-a022-4bdd3a1b4453
let
	countries |>
		@vlplot(
			title = "GDP per Capita by Region, 2019",
			width = 600, height = 300,
			:boxplot,
			transform = [{filter = "datum.GDPperCapita <= 100000"}],
			x = {:Region, title = "Region", axis = {labelAngle = 0}},
			y = {:GDPperCapita, title = "GDP per Capita [USD]"},
			color = :Region,
			config = {background = "ghostwhite"}
		)
end

# ╔═╡ f1777ba0-7414-4963-803e-d7335fecc3a9
md"""
In order to restrict the GDP values to values <= 100.000, VegaLite uses a different logic than Gadfly. In Gadfly this is done by restricting the coordinate system on the y-axis to values between 0 and 100.000 (via `Coord` and `ymax = 100000`), but a `scale = {domain = [0, 100000]}` on the y-axis doesn't produce the desired result here.

In VegaLite the range of *data* used (i.e. the `GDPperCapita` values themselfes) have to be restricted to that range. This is done using a `filter`-command within a `transform`.

"""

# ╔═╡ 64bdab07-2d6e-4517-b581-35370390fed9
md"""
## Violin Plot: GDP per Capita by Region, 2019
"""

# ╔═╡ da36a333-1639-4a4a-a05b-75e0d04784cc
let
	countries |>
		@vlplot(
			title = "GDP per Capita by Region, 2019",
			height = 300,
			mark = {:area, orient = "horizontal"},

			transform = [
				{filter = "datum.GDPperCapita <= 100000"},
				{density = "GDPperCapita", groupby = ["Region"]}
			],

			x = {"density:q", stack = "center", impute = nothing, title = nothing,
				axis = {labels = false}},
			y = {"value:q", title = "GDP per Capita [USD]", axis = {labelAngle = 0}},

			column = {"Region",
				header = {titleOrient = "bottom", labelOrient = "bottom"}},
			color = :Region,

			config = {background = "ghostwhite"},
			width = 120, spacing = 0
		)
end

# ╔═╡ 72b9a687-30ef-47ea-81d7-ba563eaeb714
md"""
As VegaLite doesn't support violin plots as a `mark` on it's own, they have to be constructed using density plots and listed in a row using the `column` dimension (here: five density plots, one for each region).
"""

# ╔═╡ 7bb96854-2752-488f-b5fd-55dcfd3077ce
md"""
*Note*: The GDP values on the y-axis in this diagram are restricted to values <= 100.000. To get the same result in the box plot above, the line starting with `transform` has to be uncommented.
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
VegaLite = "112f6efa-9a02-5b7d-90c0-432ed331239a"

[compat]
CSV = "~0.8.5"
DataFrames = "~1.2.2"
PlutoUI = "~0.7.9"
VegaLite = "~2.6.0"
"""

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

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

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

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

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

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

[[ConstructionBase]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f74e9d5388b8620b4cee35d4c5a618dd4dc547f4"
uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
version = "1.3.0"

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

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

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

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

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

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

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

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

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

[[FilePaths]]
deps = ["FilePathsBase", "MacroTools", "Reexport", "Requires"]
git-tree-sha1 = "919d9412dbf53a2e6fe74af62a73ceed0bce0629"
uuid = "8fc22ac5-c921-52a6-82fd-178b2807b824"
version = "0.8.3"

[[FilePathsBase]]
deps = ["Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "0f5e8d0cb91a6386ba47bd1527b240bd5725fbae"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.10"

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

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

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "MbedTLS", "Sockets"]
git-tree-sha1 = "c7ec02c4c6a039a98a15f955462cd7aea5df4508"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.8.19"

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

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

[[InvertedIndices]]
deps = ["Test"]
git-tree-sha1 = "15732c475062348b0165684ffe28e85ea8396afc"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.0.0"

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

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

[[JSONSchema]]
deps = ["HTTP", "JSON", "ZipFile"]
git-tree-sha1 = "b84ab8139afde82c7c65ba2b792fe12e01dd7307"
uuid = "7d188eb4-7ad8-530c-ae41-71a32a6d4692"
version = "0.3.3"

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

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

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

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

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

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

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

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

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

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

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

[[NodeJS]]
deps = ["Pkg"]
git-tree-sha1 = "905224bbdd4b555c69bb964514cfa387616f0d3a"
uuid = "2bd173c7-0d6d-553b-b6af-13a54713934c"
version = "1.3.0"

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

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "bfd7d8c7fd87f04543810d9cbd3995972236ba1b"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.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"

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

[[PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "cde4ce9d6f33219465b55162811d8de8139c0414"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.2.1"

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

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

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

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

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

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

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

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

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

[[Setfield]]
deps = ["ConstructionBase", "Future", "MacroTools", "Requires"]
git-tree-sha1 = "fca29e68c5062722b5b4435594c3d1ba557072a3"
uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
version = "0.7.1"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[Vega]]
deps = ["DataStructures", "DataValues", "Dates", "FileIO", "FilePaths", "IteratorInterfaceExtensions", "JSON", "JSONSchema", "MacroTools", "NodeJS", "Pkg", "REPL", "Random", "Setfield", "TableTraits", "TableTraitsUtils", "URIParser"]
git-tree-sha1 = "43f83d3119a868874d18da6bca0f4b5b6aae53f7"
uuid = "239c3e63-733f-47ad-beb7-a12fde22c578"
version = "2.3.0"

[[VegaLite]]
deps = ["Base64", "DataStructures", "DataValues", "Dates", "FileIO", "FilePaths", "IteratorInterfaceExtensions", "JSON", "MacroTools", "NodeJS", "Pkg", "REPL", "Random", "TableTraits", "TableTraitsUtils", "URIParser", "Vega"]
git-tree-sha1 = "3e23f28af36da21bfb4acef08b144f92ad205660"
uuid = "112f6efa-9a02-5b7d-90c0-432ed331239a"
version = "2.6.0"

[[ZipFile]]
deps = ["Libdl", "Printf", "Zlib_jll"]
git-tree-sha1 = "c3a5637e27e914a7a445b8d0ad063d701931e9f7"
uuid = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea"
version = "0.9.3"

[[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:
# ╟─f4588c88-48d4-473e-b0d7-daa9ba1f0a00
# ╟─61c31884-4e1a-4f71-aa47-a0573d3e341a
# ╟─d8a90b71-b2ce-4766-83a7-f134ebd91298
# ╟─8d71d310-4db1-4921-8c86-26221f2ceb5a
# ╠═5932694a-acae-4a1c-a2bb-ba0ec538c467
# ╠═d1bbe6fa-1f02-4158-a91f-7f4b3d641256
# ╟─f791a7b9-549c-4f2d-b493-8bd3da683f3b
# ╟─157e4e4a-a7d6-4c98-8843-d837a2232f36
# ╟─ed0b8989-8f4e-4890-909e-218761ba4ea2
# ╠═17ea91b1-db35-43b1-8cf0-f91f070b3049
# ╟─17ebb9cc-1080-41f6-a914-f73afbb8e957
# ╠═afe3969f-6ff9-4c12-88ce-62dc2c3b2b3e
# ╟─9f9924c4-a937-4dea-b045-ba8101c2971b
# ╟─a2fa5c60-aaf1-4a12-bb17-d260c41b85a6
# ╠═55b9aa0b-da52-4a8c-a7b8-4bf25da2ca2f
# ╟─ea4f0f45-7fb3-4faa-bdb9-5b4a1df04280
# ╠═94bfd9c5-2877-417e-a18c-4983f243fc85
# ╟─0d84e074-af24-4823-a355-57a3485edc36
# ╠═9f72a8e5-8777-40c0-a6af-75b1a2aa7115
# ╟─371dad4d-963c-42c9-a785-c477e124bcb4
# ╟─305f5666-2880-43dd-ae0b-7778ca4be212
# ╟─f7e54982-88f4-4210-9343-ed187f65fa33
# ╠═b38a4a8a-b655-4756-b396-bd6b2f97579c
# ╟─e0b07b43-4446-43f5-814a-0642ece326e7
# ╠═1dc24d9c-8d8d-479c-b9e1-79792ee2ab0c
# ╟─ec791080-6c52-4c80-b7ca-43b604d4ede1
# ╟─389e95fc-d410-481e-a753-f36ec70fa931
# ╠═f2753a4d-7bcd-4866-b59c-5568385e88c3
# ╟─5782433d-5ab6-4a9c-a6fd-04610edb1dd8
# ╟─e46e377d-a637-4ba8-a9ba-789d6e6415f3
# ╠═56892be5-2c1c-40a8-b0cd-d958d71a8845
# ╟─a6ce4860-949b-4e84-b884-6090b40c3b96
# ╟─11e4f95f-3849-4260-aed4-e0ac29cad61f
# ╟─8e15d66c-0e8f-4842-bc6c-65f0c49934b7
# ╠═8580f7b7-564d-4469-a022-4bdd3a1b4453
# ╟─f1777ba0-7414-4963-803e-d7335fecc3a9
# ╟─64bdab07-2d6e-4517-b581-35370390fed9
# ╠═da36a333-1639-4a4a-a05b-75e0d04784cc
# ╟─72b9a687-30ef-47ea-81d7-ba563eaeb714
# ╟─7bb96854-2752-488f-b5fd-55dcfd3077ce
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: roland-KA/StatisticalPlotsWithJulia

max_issues_repo_name: roland-KA/StatisticalPlotsWithJulia

max_forks_repo_name: roland-KA/StatisticalPlotsWithJulia


text:
#=
Tools for matching algorithms.

Author: Daisuke Oyama

=#

# random_prefs

"""
    random_prefs([rng, ]m, n, allow_unmatched=true)

Generate random preference order lists for two groups, say, m males and n
females.

Each male has a preference order over femals [1, ..., n] and "unmatched" which
is represented by 0, while each female has a preference order over males
[1, ..., m] and "unmatched" which is again represented by 0.

The optional `rng` argument specifies a random number generator.

# Arguments

- `m::Integer` : Number of males.
- `n::Integer` : Number of females.
- `allow_unmatched::Bool(true)` : If false, return preference order lists of
males and females where 0 is always placed in the last place, (i.e.,
"unmatched" is least preferred by every individual).

# Returns

- `m_prefs::Matrix{Int}` :  Array of shape (n+1, m), where each column contains
a random permutation of 0, 1, ..., n.
- `f_prefs::Matrix{Int}` :  Array of shape (m+1, n), where each column contains
a random permutation of 0, 1, ..., m.

# Examples

```julia
julia> m_prefs, f_prefs = random_prefs(4, 3);

julia> m_prefs
4x4 Array{Int64,2}:
 3  3  1  3
 0  2  3  1
 2  1  2  0
 1  0  0  2

julia> f_prefs
5x3 Array{Int64,2}:
 1  2  4
 4  3  1
 3  4  2
 0  0  0
 2  1  3

julia> m_prefs, f_prefs = random_prefs(4, 3, allow_unmatched=false);

julia> m_prefs
4x4 Array{Int64,2}:
 1  3  1  2
 2  1  3  3
 3  2  2  1
 0  0  0  0

julia> f_prefs
5x3 Array{Int64,2}:
 1  2  3
 2  3  4
 4  1  1
 3  4  2
 0  0  0
```
"""
function random_prefs(rng::AbstractRNG,
	                  m::Integer, n::Integer; allow_unmatched::Bool=true)
    m_prefs = _random_prefs(rng, m, n, allow_unmatched)
    f_prefs = _random_prefs(rng, n, m, allow_unmatched)
    return m_prefs, f_prefs
end

random_prefs(m::Integer, n::Integer; allow_unmatched::Bool=true) =
    random_prefs(Base.GLOBAL_RNG, m, n, allow_unmatched=allow_unmatched)

function _random_prefs(r::AbstractRNG,
	                   m::Integer, n::Integer, allow_unmatched::Bool)
    prefs = Array(Int, n+1, m)
    for j in 1:m
        prefs[end, j] = 0
    end

    _randperm2d!(r, sub(prefs, 1:n, :))

    if allow_unmatched
        unmatched_rankings = rand(r, 2:n+1, m)
        for j in 1:m
            prefs[end, j] = prefs[unmatched_rankings[j], j]
            prefs[unmatched_rankings[j], j] = 0
        end
    end

    return prefs
end

# Copied from base/random.jl
"Return a random `Int` (masked with `mask`) in ``[0, n)``, when `n <= 2^52`."
@inline function rand_lt(r::AbstractRNG, n::Int, mask::Int=nextpow2(n)-1)
    # this duplicates the functionality of RangeGenerator objects,
    # to optimize this special case
    while true
        x = (Base.Random.rand_ui52_raw(r) % Int) & mask
        x < n && return x
    end
end

# In-place version of randperm in base/random.jl
"""
Given a vector `a` of length n, generate a random permutation of length n and
store it in `a`.
"""
function _randperm!{T<:Integer}(r::AbstractRNG, a::AbstractVector{T})
    # a = Array{typeof(n)}(n)
    n = length(a)
    @assert n <= Int64(2)^52
    if n == 0
       return a
    end
    a[1] = 1
    mask = 3
    @inbounds for i = 2:Int(n)
        j = 1 + rand_lt(r, i, mask)
        if i != j # a[i] is uninitialized (and could be #undef)
            a[i] = a[j]
        end
        a[j] = i
        i == 1+mask && (mask = 2mask + 1)
    end
    return a
end

"""
Given an m x n matrix `a`, generate n random permutations of length m and store
them in columns of `a`.
"""
function _randperm2d!{T<:Integer}(r::AbstractRNG, a::AbstractMatrix{T})
    m, n = size(a)
    for j in :1:n
        _randperm!(r, sub(a, :, j))
    end
    a
end

ext: jl


text:
@testset "Array math" begin
    @testset "AbstractArray-of-StaticArray with scalar math" begin
        v = SVector{2,Float64}[SVector{2,Float64}(1,1)]
        @test v .* 1.0 == v
        @test typeof(v .* 1.0) == typeof(v)
        @test 1 .- v == v .- v
        @test typeof(1 .- v) == typeof(v)
        v2 = SVector{2,Int}[SVector{2,Int}(1,1)]
        @test v2 .* 1.0 == v
        @test typeof(v2 .* 1.0) == typeof(v)
    end

    @testset "Array-scalar math" begin
        m = @SMatrix [1 2; 3 4]

        @test m .+ 1 === @SMatrix [2 3; 4 5]
        @test 1 .+ m === @SMatrix [2 3; 4 5]
        @test m .* 2 === @SMatrix [2 4; 6 8]
        @test 2 .* m === @SMatrix [2 4; 6 8]
        @test m .- 1 === @SMatrix [0 1; 2 3]
        @test 1 .- m === @SMatrix [0 -1; -2 -3]
        @test m ./ 2 === @SMatrix [0.5 1.0; 1.5 2.0]
        @test 12 ./ m === @SMatrix [12.0 6.0; 4.0 3.0]

    end

    @testset "Elementwise array math" begin
        m1 = @SMatrix [1 2; 3 4]
        m2 = @SMatrix [4 3; 2 1]

        @test .-(m1) === @SMatrix [-1 -2; -3 -4]

        @test m1 .+ m2 === @SMatrix [5 5; 5 5]
        @test m1 .* m2 === @SMatrix [4 6; 6 4]
        @test m1 .- m2 === @SMatrix [-3 -1; 1 3]
        @test m1 ./ m2 === @SMatrix [0.25 2/3; 1.5 4.0]
    end

    @testset "zeros() and ones()" begin
        @test zeros(SVector{3,Float64}) === @SVector [0.0, 0.0, 0.0]
        @test zeros(SVector{3,Int}) === @SVector [0, 0, 0]
        @test ones(SVector{3,Float64}) === @SVector [1.0, 1.0, 1.0]
        @test ones(SVector{3,Int}) === @SVector [1, 1, 1]

        @test zeros(SVector{3}) === @SVector [0.0, 0.0, 0.0]
        @test zeros(SMatrix{2,2}) === @SMatrix [0.0 0.0; 0.0 0.0]
        @test zeros(SArray{(1,1,1)}) === SArray{(1,1,1)}((0.0,))
        @test zeros(MVector{3})::MVector == @MVector [0.0, 0.0, 0.0]
        @test zeros(MMatrix{2,2})::MMatrix == @MMatrix [0.0 0.0; 0.0 0.0]
        @test zeros(MArray{(1,1,1)})::MArray == MArray{(1,1,1)}((0.0,))

        @test ones(SVector{3}) === @SVector [1.0, 1.0, 1.0]
        @test ones(SMatrix{2,2}) === @SMatrix [1.0 1.0; 1.0 1.0]
        @test ones(SArray{(1,1,1)}) === SArray{(1,1,1)}((1.0,))
        @test ones(MVector{3})::MVector == @MVector [1.0, 1.0, 1.0]
        @test ones(MMatrix{2,2})::MMatrix == @MMatrix [1.0 1.0; 1.0 1.0]
        @test ones(MArray{(1,1,1)})::MArray == MArray{(1,1,1)}((1.0,))
    end
end


text:
function readvar(path::String, name::String)
    raw_data = ncread(path, name)
    scale_factor = ncgetatt(path, name, "scale_factor")
    if scale_factor != Void()
        raw_data *= scale_factor
    end

    raw_data[:]
end

# # TODO: generalize(?) ===>> this code is slower than it should
# # need to make the cartesian product of the dimensions
# function cartesianproduct(long::Array{Float32}, lat::Array{Float32}, time::Array{DateTime})
#     longlat = [repeat(long, inner=length(lat)) repeat(lat, outer=length(long))]
#
#     [repeat(view(longlat, :, 1), inner=length(time)) repeat(view(longlat, :, 2), inner=length(time)) repeat(time, outer=length(view(longlat, :, 1)))]
# end
#
# # returns a dict of dimension => data
# function readdims(path::String)
#     long_raw = readvar(path, "longitude")
#     lat_raw = readvar(path, "latitude")
#     time_raw = map(x -> Base.Dates.DateTime(1950, 1, 1) + Base.Dates.Day(x), readvar(path, "time"))
#     dims = cartesianproduct(long_raw, lat_raw, time_raw)
#     d = OrderedDict()
#     d["longitude"] = Array{Float32}(view(dims, :, 1))
#     d["latitude"] = Array{Float32}(view(dims, :, 2))
#     d["time"] = Array{DateTime}(view(dims, :, 3))
#
#     d
# end

# looped version, significantly faster for some reason
function cartesianproduct(long::Array{Float32}, lat::Array{Float32}, time::Array{Int32})
    length_long = length(long)
    length_lat = length(lat)
    length_time = length(time)
    total_size = length_long * length_lat * length_time

    d = OrderedDict()

    # inner x2
    new_long = zeros(Float32, total_size)
    times = length_lat * length_time
    index = 1
    for i = 1:length_long
        for j = 1:times
            new_long[index] = long[i]
            index += 1
        end
    end

    d["longitude"] = new_long

    # outer + inner
    new_lat_temp = zeros(Float32, length_long * length_lat)
    index = 1
    for i = 1:length_long
        for j = 1:length_lat
            new_lat_temp[index] = lat[j]
            index += 1
        end
    end
    length_new_lat_temp = length(new_lat_temp)
    new_lat = zeros(Float32, total_size)
    index = 1
    for i = 1:length_new_lat_temp
        for j = 1:length_time
            new_lat[index] = new_lat_temp[i]
            index += 1
        end
    end

    d["latitude"] = new_lat

    # outer x2
    new_time = zeros(Int32, total_size)
    times = length_long * length_lat
    index = 1
    for i = 1:times
        for j = 1:length_time
            new_time[index] = time[j]
            index += 1
        end
    end

    d["time"] = new_time

    d
end

# returns a dict of dimension => data
function readdims(path::String)
    long_raw = readvar(path, "longitude")
    lat_raw = readvar(path, "latitude")
    time_raw = readvar(path, "time")
    d = cartesianproduct(long_raw, lat_raw, time_raw)
    d["time"] = map(x -> Base.Dates.DateTime(1950, 1, 1) + Base.Dates.Day(x), d["time"])

    d
end

# returns a dict of variable => data, for all variables not dimensions
function readdata(path::String)
    info = ncinfo(path)

    data_names = [k for k in keys(info.vars) if !(k in keys(info.dim))]

    d = Dict()
    for data_name in data_names
        d[data_name] = readvar(path, data_name)
    end

    d
end

# read the entire netcdf file as a DataFrame
function readnetcdf(path::String)
    dict = readdims(path)
    merge!(dict, readdata(path))

    DataFrame(dict)
end

max_stars_repo_name: radujica/data-analysis-pipelines

max_issues_repo_name: radujica/data-analysis-pipelines

max_forks_repo_name: radujica/data-analysis-pipelines


text:
using Neo4j
using Base.Test

@test isdefined(:Neo4j) == true
@test typeof(Neo4j) == Module

graph = nothing
try
  print("[TEST] Creating a Neo4j connection to localhost:7474 with no credentials...");
  graph = getgraph()
catch
  print("[TEST] Anonymous connection failed! Creating a Neo4j connection to localhost:7474 with neo4j:neo5j credentials...");
  #Trying with security.
  conn = Neo4j.Connection("localhost"; user="neo4j", password="neo5j");
  graph = getgraph(conn);
end
conn = graph.connection;
println("Success!");

print("[TEST] Checking version of connected graph = Neo4j ", ascii(graph.version), "...")
# Have to account for newer Neo4j! Using version text - ref: http://docs.julialang.org/en/release-0.4/manual/strings/
@test convert(VersionNumber, graph.version) >  v"2.0.0"
# Check that
@test graph.node == "http://localhost:7474/db/data/node"
println("Success!");

print("[TEST] Creating a node...");
barenode = Neo4j.createnode(graph)
@test barenode.self == "http://localhost:7474/db/data/node/$(barenode.id)"
println("Success!");

print("[TEST] Creating a node with properties...");
propnode = Neo4j.createnode(graph, Dict{AbstractString,Any}("a" => "A", "b" => 1))  #UTF8String
@test propnode.data["a"] == "A"
@test propnode.data["b"] == 1
println("Success!");

print("[TEST] Retrieving the created node...");
gotnode = getnode(graph, propnode.id)
@test gotnode.id == propnode.id
@test gotnode.data["a"] == "A"
@test gotnode.data["b"] == 1
println("Success!");

print("[TEST] Setting node properties...");
setnodeproperty(barenode, "a", "A")
barenode = getnode(barenode)
@test barenode.data["a"] == "A"
println("Success!");

print("[TEST] Getting node properties...");
props = getnodeproperties(propnode)
@test props["a"] == "A"
@test props["b"] == 1
@test length(props) == 2
println("Success!");

print("[TEST] Updating node properties...");
updatenodeproperties(barenode, Dict{AbstractString,Any}("a" => 1, "b" => "A"))  #UTF8String
barenode = getnode(barenode)
@test barenode.data["a"] == 1
@test barenode.data["b"] == "A"
println("Success!");

print("[TEST] Deleting node properties...");
deletenodeproperties(barenode)
barenode = getnode(barenode)
@test length(barenode.data) == 0
println("Success!");

print("[TEST] Deleting a specific property...");
deletenodeproperty(propnode, "b")
propnode = getnode(propnode)
@test length(propnode.data) == 1
@test propnode.data["a"] == "A"
println("Success!");

print("[TEST] Adding a node label...")
addnodelabel(barenode, "A")
barenode = getnode(barenode)
@test getnodelabels(barenode) == ["A"]
println("Success!");

print("[TEST] Adding multiple node labels...")
addnodelabels(barenode, ["B", "C"])
barenode = getnode(barenode)
labels = getnodelabels(barenode)
@test "A" in labels
@test "B" in labels
@test "C" in labels
@test length(labels) == 3
println("Success!");

print("[TEST] Updating node labels...")
updatenodelabels(barenode, ["D", "E", "F"])
barenode = getnode(barenode)
labels = getnodelabels(barenode)
@test "D" in labels
@test "E" in labels
@test "F" in labels
@test length(labels) == 3
println("Success!");

print("[TEST] Deleting a node label...")
deletenodelabel(barenode, "D")
barenode = getnode(barenode)
labels = getnodelabels(barenode)
@test "E" in labels
@test "F" in labels
@test length(labels) == 2
println("Success!");

print("[TEST] Getting nodes for a given label...")
nodes = getnodesforlabel(graph, "E")
@test length(nodes) > 0
@test barenode.id in [n.id for n = nodes]
println("Success!");

print("[TEST] Getting all labels...")
labels = getlabels(graph)
# TODO Can't really test this because there might be other crap in the local DB
println("Success!");

print("[TEST] Creating a relationship...")
rel1 = createrel(barenode, propnode, "test"; props=Dict{AbstractString,Any}("a" => "A", "b" => 1)); #UTF8String
rel1alt = getrel(graph, rel1.id);
@test rel1.reltype == "TEST"
@test rel1.data["a"] == "A"
@test rel1.data["b"] == 1
@test rel1.id == rel1alt.id
println("Success!");

print("[TEST] Getting relationships from nodes...")
endnode = Neo4j.createnode(graph, Dict{AbstractString,Any}("a" => "A", "b" => 1))  # UTF8String
rel2 = createrel(propnode, endnode, "test"; props=Dict{AbstractString,Any}("a" => "A", "b" => 1));  # UTF8String
@test length(Neo4j.getrels(endnode)) == 1
@test length(Neo4j.getrels(propnode)) == 2
@test length(Neo4j.getrels(barenode)) == 1
@test length(Neo4j.getrels(endnode, incoming=true, outgoing=false)) == 1
@test length(Neo4j.getrels(endnode, incoming=false, outgoing=true)) == 0
@test length(Neo4j.getrels(propnode, incoming=true, outgoing=false)) == 1
@test length(Neo4j.getrels(propnode, incoming=false, outgoing=true)) == 1
println("Success!")

print("[TEST] Getting neighbors...")
neighbors = Neo4j.getneighbors(propnode)
@test length(neighbors) == 2
neighbors = Neo4j.getneighbors(propnode, incoming=true, outgoing=false)
@test length(neighbors) == 1
@test neighbors[1].metadata["id"] == barenode.metadata["id"]
neighbors = Neo4j.getneighbors(propnode, incoming=false, outgoing=true)
@test length(neighbors) == 1
@test neighbors[1].metadata["id"] == endnode.metadata["id"]
println("Success!")

print("[TEST] Getting relationship properties...")
rel1prop = getrelproperties(rel1);
@test rel1prop["a"] == "A"
@test rel1prop["b"] == 1
@test length(rel1prop) == 2
@test getrelproperty(rel1, "a") == "A"
@test getrelproperty(rel1, "b") == 1
println("Success!");

print("[TEST] Deleting a relationship...")
deleterel(rel1)
deleterel(rel2)
@test_throws ErrorException getrel(graph, rel1.id)
@test_throws ErrorException getrel(graph, rel2.id)
println("Success!");

print("[TEST] Deleting a node...")
deletenode(graph, barenode.id)
deletenode(graph, propnode.id)
@test_throws ErrorException getnode(graph, barenode.id)
@test_throws ErrorException getnode(graph, propnode.id)
println("Success!");

# --- New transaction code from Glesica source ---

print("[TEST] Creating a transaction and nodes with the new Glesica transaction framework...")
loadtx = transaction(conn)

function createnode(txn, name, age; submit=false)
  q = "CREATE (n:Neo4jjl) SET n.name = {name}, n.age = {age}"
  txn(q, "name" => name, "age" => age; submit=submit)
end

@test length(loadtx.statements) == 0

createnode(loadtx, "John Doe", 20)

@test length(loadtx.statements) == 1

createnode(loadtx, "Jane Doe", 20)

@test length(loadtx.statements) == 2
println("Success!")

query = "MATCH (n:Neo4jjl) WHERE n.age = {age} RETURN n.name";
print("[TEST] Doing a match query '", query, "'...")
people = loadtx(query, "age" => 20; submit=true)

@test length(loadtx.statements) == 0
@test length(people.results) == 3
@test length(people.errors) == 0

matchresult = people.results[3]
@test matchresult["columns"][1] == "n.name"
@test "John Doe" in [row["row"][1] for row = matchresult["data"]]
@test "Jane Doe" in [row["row"][1] for row = matchresult["data"]]

loadresult = commit(loadtx)

@test length(loadresult.results) == 0
@test length(loadresult.errors) == 0
println("Success!")

query = "MATCH (n:Neo4jjl) WHERE n.age = {age} DELETE n"
print("[TEST] Deleting nodes '", query, "'...")

deletetx = transaction(conn)
deletetx(query, "age" => 20)

deleteresult = commit(deletetx)

@test length(deleteresult.results) == 1
@test length(deleteresult.results[1]["columns"]) == 0
@test length(deleteresult.results[1]["data"]) == 0
@test length(deleteresult.errors) == 0
println("Success!")

print("[TEST] Rolling back transactions...")

rolltx = transaction(conn)

person = createnode(rolltx, "John Doe", 20; submit=true)

@test length(rolltx.statements) == 0
@test length(person.results) == 1
@test length(person.errors) == 0

rollback(rolltx)

rolltx = transaction(conn)
rollresult = rolltx("MATCH (n:Neo4jjl) WHERE n.name = 'John Doe' RETURN n"; submit=true)

@test length(rollresult.results) == 1
@test length(rollresult.results[1]["columns"]) == 1
@test length(rollresult.results[1]["data"]) == 0
@test length(rollresult.errors) == 0

println("Success!");

# --- New cypherQuery using transaction/commit endpoint ---

print("[TEST] MATCH node and return DataFrame using cypherQuery()...")

# Open transaction and create node
loadtx = transaction(conn)
createnode(loadtx, "John Doe", 20; submit=true)
Neo4j.commit(loadtx)

matchresult = cypherQuery(conn,
                  "MATCH (n:Neo4jjl {name: {name}}) RETURN n.name AS Name, n.age AS Age;",
                  "name" => "John Doe")
@test DataFrames.DataFrame(Name = "John Doe", Age = 20) == matchresult

# Cleanup
deletetx = transaction(conn)
deletetx(query, "age" => 20)
deleteresult = commit(deletetx)

println("Success!");
println("--- All tests passed!");


text:
global charNum = parse(Int, ARGS[1])
include("loadData.jl")

##

char = @pipe d |>
      zip(_[:, :taxon], _[:, charNum+1]) |>
      Dict
##

nnodes = [length(post_order(t)) for t in ttrees[1]]

data = zeros(nbase, nsites, maximum(nnodes), length(ttrees[1]))

for i in 1:length(taxa)
      for l in taxa[i]
            s = char[l]
            t_ind = find_by_name(ttrees[1][i],l).num
            if s in states
                  data[:,1,t_ind, i] = (states .== s)
            else
                  data[:,1,t_ind, i] .= 1
            end
      end
end

##

nLineages = length(lineages)

my_data = Dict{Symbol, Any}(
      :data => data,
      :nLineages => nLineages
)

##

model = Model(
      ind = Stochastic(0, () -> Logistic(), true),
      data = Stochastic(
            4,
            (fullSrates, nLineages, ind) -> MultiplePhyloDist(
                  ttrees[Integer(ceil(1000 * invlogit(ind[])))],
                  fullSrates,
                  ones(1, nLineages),
                  freeK,
            ),
            false,
      ),
      s_rates = Stochastic(
            2,
            nLineages -> UnivariateDistribution[
                  Normal(0,1) for i in 1:12, j in 1:nLineages
            ],
            true,
      ),
      fullSrates = Logical(
            2,
            s_rates -> exp.(s_rates),
            false,
      ),
)

##

inits = [Dict{Symbol, Any}(
      :data => data,
      :nLineages => nLineages,
      :s_rates => rand(Normal(0,1), (12, nLineages)),
      :ind => (a = zeros(); a[]=rand(Logistic()); a),
) for i in 1:2]

##

scheme = [Slice(:s_rates, 1.), Empirical(:ind, 1) ]

setsamplers!(model, scheme)

##

try
      mkdir("output")
catch e
end

##

sim = mcmc(
      model,
      my_data,
      inits,
      1000,
      burnin = 0,
      thin = 10,
      chains = 2,
      trees = false,
)
bi = 1 + size(sim)[1] ÷ 2
gd = gelmandiag(sim[bi:end,:,:])
psrf = maximum(gd.value[:,1])
write("output/lineage_$(lpad(charNum, 2, "0")).jls", sim[bi:end,:,:])
write("output/lineage_$(lpad(charNum, 2, "0")).log", string(psrf)*"\n")

while psrf > 1.1
      global sim, psrf, bi, gd
      @show psrf
      sim = mcmc(sim, 1000)
      bi = 1 + size(sim)[1] ÷ 2
      gd = gelmandiag(sim[bi:end,:,:])
      psrf = maximum(gd.value[:,1])
      write("output/lineage_$(lpad(charNum, 2, "0")).jls", sim[bi:end,:,:])
      write("output/lineage_$(lpad(charNum, 2, "0")).log", string(psrf)*"\n")
end

max_stars_repo_name: erathorn/phylogeneticTypology

max_issues_repo_name: erathorn/phylogeneticTypology

max_forks_repo_name: erathorn/phylogeneticTypology


text:
# The code in this file is a small port from scikit-learn's and numpy's
# library which is distributed under the 3-Clause BSD license.
# The rest of DecisionTree.jl is released under the MIT license.

# written by Poom Chiarawongse <eight1911@gmail.com>

module treeregressor
    include("../util.jl")

    import Random
    export fit

    mutable struct NodeMeta{S}
        l           :: NodeMeta{S}  # right child
        r           :: NodeMeta{S}  # left child
        label       :: Float64      # most likely label
        feature     :: Int          # feature used for splitting
        threshold   :: S            # threshold value
        is_leaf     :: Bool
        depth       :: Int
        region      :: UnitRange{Int} # a slice of the samples used to decide the split of the node
        features    :: Vector{Int}    # a list of features not known to be constant
        split_at    :: Int            # index of samples
        function NodeMeta{S}(
            features :: Vector{Int},
            region   :: UnitRange{Int},
            depth    :: Int) where S
        node = new{S}()
        node.depth = depth
        node.region = region
        node.features = features
        node.is_leaf = false
        node
    end
    end

    struct Tree{S}
        root   :: NodeMeta{S}
        labels :: Vector{Int}
    end

    # find an optimal split that satisfy the given constraints
    # (max_depth, min_samples_split, min_purity_increase)
    function _split!(
            X                   :: AbstractMatrix{S}, # the feature array
            Y                   :: AbstractVector{Float64}, # the label array
            W                   :: AbstractVector{U},
            node                :: NodeMeta{S}, # the node to split
            max_features        :: Int, # number of features to consider
            max_depth           :: Int, # the maximum depth of the resultant tree
            min_samples_leaf    :: Int, # the minimum number of samples each leaf needs to have
            min_samples_split   :: Int, # the minimum number of samples in needed for a split
            min_purity_increase :: Float64, # minimum purity needed for a split
            indX                :: AbstractVector{Int}, # an array of sample indices,
                                                # we split using samples in indX[node.region]
            # the two arrays below are given for optimization purposes
            Xf                  :: AbstractVector{S},
            Yf                  :: AbstractVector{Float64},
            Wf                  :: AbstractVector{U},
            rng                 :: Random.AbstractRNG) where {S, U}

        region = node.region
        n_samples = length(region)
        r_start = region.start - 1

        @inbounds @simd for i in 1:n_samples
            Yf[i] = Y[indX[i + r_start]]
            Wf[i] = W[indX[i + r_start]]
        end

        tssq = zero(U)
        tsum = zero(U)
        wsum = zero(U)
        @inbounds @simd for i in 1:n_samples
            tssq += Wf[i]*Yf[i]*Yf[i]
            tsum += Wf[i]*Yf[i]
            wsum += Wf[i]
        end

        node.label =  tsum / wsum
        if (min_samples_leaf * 2 >  n_samples
         || min_samples_split    >  n_samples
         || max_depth            <= node.depth
          # equivalent to old_purity > -1e-7
         || tsum * node.label    > -1e-7 * wsum + tssq)
            node.is_leaf = true
            return
        end

        features = node.features
        n_features = length(features)
        best_purity = typemin(U)
        best_feature = -1
        threshold_lo = X[1]
        threshold_hi = X[1]

        indf = 1
        # the number of new constants found during this split
        n_constant = 0
        # true if every feature is constant
        unsplittable = true
        # the number of non constant features we will see if
        # only sample n_features used features
        # is a hypergeometric random variable
        total_features = size(X, 2)

        # this is the total number of features that we expect to not
        # be one of the known constant features. since we know exactly
        # what the non constant features are, we can sample at 'non_constants_used'
        # non constant features instead of going through every feature randomly.
        non_constants_used = util.hypergeometric(n_features, total_features-n_features, max_features, rng)
        @inbounds while (unsplittable || indf <= non_constants_used) && indf <= n_features
            feature = let
                indr = rand(rng, indf:n_features)
                features[indf], features[indr] = features[indr], features[indf]
                features[indf]
            end

            rssq = tssq
            lssq = zero(U)
            rsum = tsum
            lsum = zero(U)

            @simd for i in 1:n_samples
                Xf[i] = X[indX[i + r_start], feature]
            end

            # sort Yf and indX by Xf
            util.q_bi_sort!(Xf, indX, 1, n_samples, r_start)
            @simd for i in 1:n_samples
                Yf[i] = Y[indX[i + r_start]]
                Wf[i] = W[indX[i + r_start]]
            end
            nl, nr = zero(U), wsum
            # lo and hi are the indices of
            # the least upper bound and the greatest lower bound
            # of the left and right nodes respectively
            hi = 0
            last_f = Xf[1]
            is_constant = true
            while hi < n_samples
                lo = hi + 1
                curr_f = Xf[lo]
                hi = (lo < n_samples && curr_f == Xf[lo+1]
                    ? searchsortedlast(Xf, curr_f, lo, n_samples, Base.Order.Forward)
                    : lo)

                (lo != 1) && (is_constant = false)
                # honor min_samples_leaf
                if lo-1 >= min_samples_leaf && n_samples - (lo-1) >= min_samples_leaf
                    unsplittable = false
                    purity = (rsum * rsum / nr) + (lsum * lsum / nl)
                    if purity > best_purity && !isapprox(purity, best_purity)
                        # will take average at the end, if possible
                        threshold_lo = last_f
                        threshold_hi = curr_f
                        best_purity  = purity
                        best_feature = feature
                    end
                end

                # update, lssq, rssq, lsum, rsum in the direction
                # that would require the smaller number of iterations
                if (hi << 1) < n_samples + lo # i.e., hi - lo < n_samples - hi
                    @simd for i in lo:hi
                        nr   -= Wf[i]
                        rsum -= Wf[i]*Yf[i]
                        rssq -= Wf[i]*Yf[i]*Yf[i]
                    end
                else
                    nr = rsum = rssq = zero(U)
                    @simd for i in (hi+1):n_samples
                        nr   += Wf[i]
                        rsum += Wf[i]*Yf[i]
                        rssq += Wf[i]*Yf[i]*Yf[i]
                    end
                end
                lsum = tsum - rsum
                lssq = tssq - rssq
                nl   = wsum - nr

                last_f = curr_f
            end

            # keep track of constant features to be used later.
            if is_constant
                n_constant += 1
                features[indf], features[n_constant] = features[n_constant], features[indf]
            end

            indf += 1
        end

        # no splits honor min_samples_leaf
        @inbounds if (unsplittable
                || best_purity - tsum * node.label < min_purity_increase * wsum)
            node.is_leaf = true
            return
        else
            # new_purity - old_purity < stop.min_purity_increase
            @simd for i in 1:n_samples
                Xf[i] = X[indX[i + r_start], best_feature]
            end

            try
                node.threshold = (threshold_lo + threshold_hi) / 2.0
            catch
                node.threshold = threshold_hi
            end
            # split the samples into two parts: ones that are greater than
            # the threshold and ones that are less than or equal to the threshold
            #                                 ---------------------
            # (so we partition at threshold_lo instead of node.threshold)
            node.split_at = util.partition!(indX, Xf, threshold_lo, region)
            node.feature = best_feature
            node.features = features[(n_constant+1):n_features]
        end

    end

    @inline function fork!(node :: NodeMeta{S}) where S
        ind = node.split_at
        region = node.region
        features = node.features
        # no need to copy because we will copy at the end
        node.l = NodeMeta{S}(features, region[    1:ind], node.depth + 1)
        node.r = NodeMeta{S}(features, region[ind+1:end], node.depth + 1)
    end

    function _fit(
            X                     :: AbstractMatrix{S},
            Y                     :: AbstractVector{Float64},
            W                     :: AbstractVector{U},
            max_features          :: Int,
            max_depth             :: Int,
            min_samples_leaf      :: Int,
            min_samples_split     :: Int,
            min_purity_increase   :: Float64,
            rng=Random.GLOBAL_RNG :: Random.AbstractRNG) where {S, U}

        n_samples, n_features = size(X)

        Yf  = Array{Float64}(undef, n_samples)
        Xf  = Array{S}(undef, n_samples)
        Wf  = Array{U}(undef, n_samples)

        indX = collect(1:n_samples)
        root = NodeMeta{S}(collect(1:n_features), 1:n_samples, 0)
        stack = NodeMeta{S}[root]

        @inbounds while length(stack) > 0
            node = pop!(stack)
            _split!(
                X, Y, W,
                node,
                max_features,
                max_depth,
                min_samples_leaf,
                min_samples_split,
                min_purity_increase,
                indX,
                Xf, Yf, Wf,
                rng)
            if !node.is_leaf
                fork!(node)
                push!(stack, node.r)
                push!(stack, node.l)
            end
        end
        return (root, indX)
    end

    function fit(;
            X                     :: AbstractMatrix{S},
            Y                     :: AbstractVector{Float64},
            W                     :: Union{Nothing, AbstractVector{U}},
            max_features          :: Int,
            max_depth             :: Int,
            min_samples_leaf      :: Int,
            min_samples_split     :: Int,
            min_purity_increase   :: Float64,
            rng=Random.GLOBAL_RNG :: Random.AbstractRNG) where {S, U}

        n_samples, n_features = size(X)
        if W == nothing
            W = fill(1.0, n_samples)
        end

        util.check_input(
            X,
            Y,
            W,
            max_features,
            max_depth,
            min_samples_leaf,
            min_samples_split,
            min_purity_increase)

        root, indX = _fit(
            X,
            Y,
            W,
            max_features,
            max_depth,
            min_samples_leaf,
            min_samples_split,
            min_purity_increase,
            rng)

        return Tree{S}(root, indX)

    end
end


text:
module JuCLI

export CLI, ActionContext, Flag, run

struct Flag{T}
    name::String
    usage::String
    optional::Bool
    default::Union{T,Nothing}
end

# ActionContext is passed to the action function.
# Note: The action function can modify this for
# sub-commands. For example, calling `popfirst!`
# on `ActionContext.args` to pick the next action.
struct ActionContext
    args::Vector{String}
    flags::Vector{Flag}

    ActionContext(
        args::Vector{String}=Vector{String}(),
        flags::Vector{Flag}=Vector{Flag}()
    ) = new(args, flags)
end

struct CLI
    name::String
    usage::String
    flags::Vector{Flag}
    action::Function
    autohelp::Bool

    CLI(
        name="cli",
        usage="",
        flags=Vector{Flag}(),
        action=() -> nothing,
        autohelp=true
    ) = new(name, usage, flags, action, autohelp)
end

# run the CLI with the supplied arguments.
# This is typically expected to be called with
# the `ARGS` data.
function run(cli::CLI, args::String[])
    ctx = parseargs(args)
end

# Create a new CLI context from the supplied arguments.
function parseargs(args::String[])::ActionContext

end

function formathelp(cli::CLI)::String
    """
    NAME:
        $(cli.name) - $(cli.usage)
    """
end

function displayhelp(cli::CLI)
    println(formathelp(cli))
end

# Check if the user
function checkforhelp(args::String[])::Bool

end

end # module


text:
using FluxExtensions, LinearAlgebra, MLDataPattern

FluxExtensions.log_normal(x, μ, σ2::T, d::Int) where {T<:Real}  = (d > 0) ? - sum((@. ((x - μ)^2) ./ σ2), dims=1) / 2 .- d * (log.(σ2) + log(2π)) / 2 : 0

"""
	LinearAlgebra.det(f, x::Vector)
	determinant of a jacobian of `f` at point `x`
"""
function LinearAlgebra.det(f, x::Vector)
	Σ = Flux.data(Flux.Tracker.jacobian(f, x))
	S = svd(Σ)
	2*sum(log.(S.S .+ 1f-6))
end

jacodeco(model, x::Vector, z::Vector) = sum(log_normal(z)) + sum(log_pxexpectedz(model, x, z)) - det(model.g, z)
function jacodeco(m, x, z)
	@assert nobs(x) == nobs(z)
	[jacodeco(m, getobs(x, i), getobs(z, i)) for i in 1:nobs(x)]
end
jacodeco(m, x) = jacodeco(m, x, m.g(x))

function correctedjacodeco(model, x::Vector, z::Vector)
	Σ = Flux.data(Flux.Tracker.jacobian(model.g, z))
	S = svd(Σ)
	J = inv(S)
	logd = 2*sum(log.(S.S .+ 1f-6))
	if (det(I + J * transpose(J)) < 0)
		# println("determinant = $(det(I + J * transpose(J))) was negative, returning jacodeco = -Inf")
		# println(J)
		# println(J * transpose(J))
		return -Inf32
	end
	logpz = log_normal(z, zeros(Float32, size(z)...), I + J * J')[1] # TODO I changed the transpose(J) to J' - might be an issue
	logpx = FluxExtensions.log_normal(x, x_from_z(model, z).data, 1, length(x) - length(z))[1]
	logpz + logpx - logd
end

function correctedjacodeco(m, x, z)
	@assert nobs(x) == nobs(z)
	[correctedjacodeco(m, getobs(x, i), getobs(z, i)) for i in 1:nobs(x)]
end

correctedjacodeco(m, x) = correctedjacodeco(m, x, manifoldz(model, x))

"""
	manifoldz(model, x, steps = 100, rand_init = false)
	find the point in the latent space `z` such that it minimizes the
	reconstruction error. This means that it finds the point on manifold
	defined by `model.g(z)` and x
"""
function manifoldz(model, x, steps = 100, z = Flux.data(zparams(model,x)[1]))
	z = param(z)
	ps = Flux.Tracker.Params([z])
	# println("Size z: $(size(z)), size x: $(size(x))")
	li = Flux.data(mean(log_pxexpectedz(model, x, z)))
	Flux.train!((i) -> -mean(log_pxexpectedz(model, x, z)), ps, 1:steps, ADAM())
	le = Flux.data(mean(log_pxexpectedz(model, x, z)))
	@info "initial = $(li) final = $(le)"
	Flux.data(z)
end

max_stars_repo_name: honzabim/FewShotAnomalyDetection

max_issues_repo_name: honzabim/FewShotAnomalyDetection

max_forks_repo_name: honzabim/FewShotAnomalyDetection


text:
using Gtk
"""
TODO:
- find @guarded (via finding its Package )

"""
@testset "misc" begin

    unhandled = convert(Cint, false)

    foo1 = @guarded (x, y) -> x + y     #guarded (private? )   #ERROR: @guarded not defined (which package isa @guarded ?)
    bar1 = @guarded (x, y) -> x + y + k
    @guarded foo2(x, y) = x + y
    @guarded bar2(x, y) = x + y + k
    @guarded function foo3(x, y)
        x + y
    end
    @guarded function bar3(x, y)
        x + y + k
    end
    @guarded unhandled function bar4(x, y)
        x + y + k
    end

    @test foo1(3, 5) == 8
    @test @test_logs (:warn, "Error in @guarded callback") bar1(3, 5) === nothing #FP  #Debugging: UncommentMe
    @test @test_logs (:warn, "Error in @guarded callback") bar1(3, 5) === nothing # TruePositive compare against nothing with triple ===
    @test foo2(3, 5) == 8
    @test @test_logs (:warn, "Error in @guarded callback") bar2(3, 5) === nothing #FP
    #@test @test_logs (:warn, "Error in @guarded callback") bar2(3, 5) == nothing # TruePositive      ....
    @test foo3(3, 5) == 8
    @test @test_logs (:warn, "Error in @guarded callback") bar3(3, 5) === nothing #FP
    #@test @test_logs (:warn, "Error in @guarded callback") bar3(3, 5) == nothing #TP    #Debugging: UncommentMe    ....
    @test @test_logs (:warn, "Error in @guarded callback") bar4(3, 5) == unhandled

    # Do-block syntax
    c = GtkCanvas()
    win = GtkWindow(c)
    showall(win)
    @guarded draw(c) do widget
        error("oops")
    end
    @test !isempty(c.mouse.ids)  # check storage of signal-handler ids (see GtkReactive)
    destroy(win)

    @test isa(Gtk.GdkEventKey(), Gtk.GdkEventKey)

    @testset "Eventloop control" begin
        before = Gtk.auto_idle[]

        Gtk.enable_eventloop(true)
        @test Gtk.is_eventloop_running()

        Gtk.auto_idle[] = true
        Gtk.pause_eventloop() do
            @test !Gtk.is_eventloop_running()
        end
        @test Gtk.is_eventloop_running()

        Gtk.auto_idle[] = false
        Gtk.pause_eventloop() do
            @test Gtk.is_eventloop_running()
        end
        @test Gtk.is_eventloop_running()

        Gtk.pause_eventloop(force = true) do
            @test !Gtk.is_eventloop_running()
        end
        @test Gtk.is_eventloop_running()

        Gtk.auto_idle[] = before
    end

end

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
function backtest(strat, data; mode = "train")
    记忆仓位 = @staticvar zeros(Float32, 100000)
    F, N, T = size(data.特征)
    T < 1 && return 0f0
    @unpack sim, 最大持仓, 持仓天数, 最多交易次数, 禁止平今,
        夜盘最早开仓时间, 夜盘最晚开仓时间, 夜盘最晚平仓时间,
        早盘最早开仓时间, 早盘最晚开仓时间,早盘最晚平仓时间 = strat
    最大持仓 = iszero(最大持仓) ? ncodes(data) : clamp(最大持仓, 1, ncodes(data))
    虚拟信号, 综合评分 = simulate(sim, data)
    持仓天数′ = size(虚拟信号, 1) ÷ N
    if 持仓天数′ == 1
        虚拟信号 = lockphase(虚拟信号, 持仓天数)
        综合评分 = lockphase(综合评分, 持仓天数)
    else
        持仓天数 = 持仓天数′
    end
    data = repeat(data, 持仓天数)
    @unpack 涨幅, 时间戳, 代码, 最新价, 买1价, 卖1价, 手续费率, 涨停, 跌停, 交易池 = data
    map!(fillnan, 时间戳, 时间戳)
    if size(记忆仓位, 1) != size(虚拟信号, 1)
        fill!(resize!(记忆仓位, size(虚拟信号, 1)), 0f0)
    end
    转移函数 = transition(sim)
    # TODO: merge select
    if select(sim)
        是否为ST = repeat(getfeat(data, r"is_st", errors = "ignore"), 持仓天数)
        实际仓位 = select_stocks(时间戳, 代码, 买1价, 卖1价, 涨停, 跌停, 交易池, 是否为ST, 虚拟信号, 综合评分,
                    记忆仓位, 最大持仓, 持仓天数, 最多交易次数, 转移函数, 禁止平今, 夜盘最早开仓时间, 夜盘最晚开仓时间,
                    夜盘最晚平仓时间, 早盘最早开仓时间, 早盘最晚开仓时间, 早盘最晚平仓时间)
    else
        实际仓位 = constraint(时间戳, 代码, 买1价, 卖1价, 涨停, 跌停, 虚拟信号, 记忆仓位, 转移函数, 最多交易次数, 禁止平今,
                    夜盘最早开仓时间, 夜盘最晚开仓时间, 夜盘最晚平仓时间, 早盘最早开仓时间, 早盘最晚开仓时间, 早盘最晚平仓时间)
    end
    if size(实际仓位, 2) == 1 && !any(isnan, 实际仓位)
        记忆仓位 .= 实际仓位[:, end]
    end
    if T > 1 && parseenv("CLOSE_LAST_POS", true)
        实际仓位[:, end] .= 0
    end
    pnl = if mode == "train"
        summarize_train(涨幅, 时间戳, 代码, 最新价, 买1价, 卖1价, 手续费率, 交易池, 记忆仓位, 实际仓位, 综合评分, 最大持仓)
    else
        summarize_test(涨幅, 时间戳, 代码, 最新价, 买1价, 卖1价, 手续费率, 交易池, 记忆仓位, 实际仓位, 综合评分, 虚拟信号, 最大持仓 * 持仓天数)
    end
    return pnl
end

function constraint(时间戳, 代码, 买1价, 卖1价, 涨停, 跌停, 虚拟信号, 记忆仓位, 转移函数, 最多交易次数, 禁止平今, 夜盘最早开仓时间, 夜盘最晚开仓时间, 夜盘最晚平仓时间, 早盘最早开仓时间, 早盘最晚开仓时间, 早盘最晚平仓时间)
    实际仓位 = zero(虚拟信号)
    N, T = size(实际仓位)
    交易次数 = zeros(Float32, N)
    @inbounds for t in 1:T, n in 1:N
        小时 = to_trade_hour(时间戳[n, t])
        之前仓位 = t > 1 ? 实际仓位[n, t - 1] : 记忆仓位[n]
        是否平仓 = 夜盘最晚平仓时间 <= 小时 <= 8 || 早盘最晚平仓时间 <= 小时 <= 16
        当前仓位 = 是否平仓 ? 0f0 : 转移函数(之前仓位, 虚拟信号[n, t])
        if (涨停[n, t] == 1 || isnan(卖1价[n, t])) && 当前仓位 > 之前仓位 ||
            (跌停[n, t] == 1 || isnan(买1价[n, t])) && 当前仓位 < 之前仓位 ||
            abs(当前仓位) > abs(之前仓位) && (交易次数[n] >= 最多交易次数 ||
                !(夜盘最早开仓时间 <= 小时 <= 夜盘最晚开仓时间) ||
                !(早盘最早开仓时间 <= 小时 <= 早盘最晚开仓时间)) ||
            abs(当前仓位) < abs(之前仓位) && 交易次数[n] >= 1 && 禁止平今 && !是否平仓
            当前仓位 = 之前仓位
        end
        时间差 = 时间戳[n, t] - 时间戳[n, max(1, t - 1)]
        if 小时 > 8 && 时间差 > 3600 * 5
            交易次数[n] = 0
        else
            交易次数[n] += abs(当前仓位 - 之前仓位)
        end
        if 代码[n, t] != 代码[n, min(t + 1, end)] # || 时间戳[n, t] == 0
            当前仓位 = 0f0
            交易次数[n] = 0
        end
        实际仓位[n, t] = 当前仓位
    end
    return 实际仓位
end

function select_stocks(时间戳, 代码, 买1价, 卖1价, 涨停, 跌停, 交易池, 是否为ST, 虚拟信号, 综合评分,
                记忆仓位, 最大持仓, 持仓天数, 最多交易次数, 转移函数, 禁止平今, 夜盘最早开仓时间,
                夜盘最晚开仓时间, 夜盘最晚平仓时间, 早盘最早开仓时间, 早盘最晚开仓时间, 早盘最晚平仓时间)
    实际仓位 = zero(虚拟信号)
    N, T = size(虚拟信号)
    交易次数 = zeros(Int, N)
    选股池 = zeros(Bool, N)
    for t in 1:T
        fill!(选股池, 0)
        选股数 = 最大持仓 * 持仓天数
        for n in 1:N
            小时 = to_trade_hour(时间戳[n, t])
            之前仓位 = t > 1 ? 实际仓位[n, t - 1] : 记忆仓位[n]
            是否平仓 = 夜盘最晚平仓时间 <= 小时 <= 8 || 早盘最晚平仓时间 <= 小时 <= 16
            当前仓位 = 是否平仓 ? 0f0 : 转移函数(之前仓位, 虚拟信号[n, t])
            if (涨停[n, t] == 1 || isnan(卖1价[n, t])) && 当前仓位 > 之前仓位 ||
                (跌停[n, t] == 1 || isnan(买1价[n, t])) && 当前仓位 < 之前仓位 ||
                abs(当前仓位) > abs(之前仓位) && (交易次数[n] >= 最多交易次数 ||
                    !(夜盘最早开仓时间 <= 小时 <= 夜盘最晚开仓时间) ||
                    !(早盘最早开仓时间 <= 小时 <= 早盘最晚开仓时间)) ||
                abs(当前仓位) < abs(之前仓位) && 交易次数[n] >= 1 && 禁止平今 && !是否平仓
                当前仓位 = 之前仓位
            end
            时间差 = 时间戳[n, t] - 时间戳[n, max(1, t - 1)]
            if 小时 > 8 && 时间差 > 3600 * 5
                交易次数[n] = 0
            else
                交易次数[n] += 当前仓位 != 之前仓位
            end
            if 代码[n, t] != 代码[n, min(t + 1, end)] || 时间戳[n, t] == 0
                当前仓位 = 0f0
                交易次数[n] = 0
            end
            if isnan(虚拟信号[n, t]) || 当前仓位 > 之前仓位
                选股池[n] = 1
                实际仓位[n, t] = 0
            else
                实际仓位[n, t] = 当前仓位
            end
            if 代码[n, t] != 代码[n, min(t + 1, end)] || 时间戳[n, t] == 0
                实际仓位[n, t] = 0
                选股池[n] = 0
                交易次数[n] = 0
            elseif 交易池[n, t] == 0 || 是否为ST[n, t] == 1
                实际仓位[n, t] = 0
                选股池[n] = 0
            elseif 涨停[n, t] == 1 || 跌停[n, t] == 1
                实际仓位[n, t] = 之前仓位
                选股池[n] = 0
            end
            选股数 -= 实际仓位[n, t] == 1
        end
        选股数 = min(选股数, 最大持仓)
        选股列表 = findall(选股池)
        if length(选股列表) <= 选股数
            实际仓位[选股列表, t] .= 1
        else
            选股列表索引 = partialsortperm(综合评分[选股列表, t], 1:选股数, rev = true)
            实际仓位[选股列表[选股列表索引], t] .= 1
        end
    end
    return 实际仓位
end

function summarize_core(涨幅, 时间戳, 代码, 最新价, 买1价, 卖1价, 手续费率, 交易池, 记忆仓位, 实际仓位, 综合评分, 最大持仓)
    记忆收益率 = @staticvar Dict{UInt64, Array{Float32}}()
    N, T = size(实际仓位)
    天数 = 1
    for n in 1:size(时间戳, 1)
        时间戳一行 = filter(!iszero, 时间戳[n, :])
        if !isempty(时间戳一行)
            天数 = sortednunique(unix2date, 时间戳一行)
            break
        end
    end
    复利 = get(ENV, "USE_COMPLEX", "0") == "1"
    收益率 = zero(实际仓位)
    之前仓位 = copy(记忆仓位)
    之前盈亏 = zeros(Float32, N)
    @inbounds for t in 1:T
        当前平均涨幅 = @views mean(涨幅[:, t])
        for n in 1:N
            当前仓位 = 实际仓位[n, t]
            仓位变化 = 当前仓位 - 之前仓位[n]
            买滑点 = fillnan((卖1价[n, t] - 最新价[n, t]) ⧶ 最新价[n, t])
            卖滑点 = fillnan((最新价[n, t] - 买1价[n, t]) ⧶ 最新价[n, t])
            交易成本 = fillnan(手续费率[n, t]) * abs(仓位变化)
            滑点 = ifelse(仓位变化 > 0, 买滑点, -卖滑点) * 仓位变化
            收益率[n, t] = 之前盈亏[n] - 交易成本 - 滑点
            之前仓位[n] = 当前仓位
        end
        t < T && for n in 1:N
            当前仓位 = 实际仓位[n, t]
            未来涨幅 = fillnan(涨幅[n, t + 1])
            之前盈亏[n] = 当前仓位 * 未来涨幅
        end
    end
    倍数 = 年化收益率 = 1f0
    资金曲线 = ones(Float32, T)
    if 复利
        @inbounds for t in 1:size(实际仓位, 2)
            Δ = @views sum(收益率[:, t]) ⧶ 最大持仓
            资金曲线[t] = 倍数 = (1 + Δ)  * 倍数
        end
        年化收益率 = 倍数^(240f0 / 天数) - 1
    else
        @inbounds for t in 1:T
            Δ = @views sum(收益率[:, t]) ⧶ 最大持仓
            资金曲线[t] = 倍数 = 倍数 + Δ
        end
        年化收益率 = 240f0 * (倍数 - 1f0) / 天数
    end
    最大回撤, 最大回撤期 = drawdown(资金曲线)
    夏普率 = sharperatio(资金曲线, 240T / 天数)
    length(记忆收益率) > 10 && empty!(记忆收益率)
    记忆收益率[objectid(时间戳)] = 收益率
    return 收益率, 资金曲线, 年化收益率, 最大回撤, 夏普率
end

function summarize_train(涨幅, 时间戳, 代码, 最新价, 买1价, 卖1价, 手续费率, 交易池, 记忆仓位, 实际仓位, 综合评分, 最大持仓)
    收益率, 资金曲线, 年化收益率, 最大回撤, 夏普率 =
        summarize_core(涨幅, 时间戳, 代码, 最新价, 买1价, 卖1价, 手续费率, 交易池, 记忆仓位, 实际仓位, 综合评分, 最大持仓)
    评分模式 = parseenv("SCORE", "R")
    if 评分模式 == "R" # Return
        分数 = 年化收益率
    elseif 评分模式 == "SHARPE"
        分数 = 夏普率
    elseif 评分模式 == "RoMaD" # Return Over Maximum Drawdown
        分数 = 年化收益率 ⧶ abs(最大回撤)
    elseif occursin("LAR", 评分模式)  # Loss Amplified Return
        λ = Meta.parse(split(评分模式, '_')[2])
        分数 = sum(x -> amploss(x, λ), 年化收益率)
    elseif occursin("IR", 评分模式)  # Information Rank
        nonan(x) = ifelse(isnan(x), zero(x), x)
        分数 = mean(nonan(corspearman(综合评分[:, t], 涨幅[:, t])) for t in 1:size(实际仓位, 2))
    elseif occursin("WR", 评分模式)  # Average Winning Ratio
        分数 = mean(mean(涨幅[实际仓位[:, t] .> 0, t]) >
                    mean(涨幅[交易池[:, t] .> 0, t])
                    for t in 1:(size(实际仓位, 2) - 1))
    elseif occursin("ACC", 评分模式) # Accuracy
        分数 = mean(mean(涨幅[实际仓位[:, t] .> 0, t] .>
                        mean(涨幅[交易池[:, t] .> 0, t]))
                    for t in 1:(size(实际仓位, 2) - 1))
    end
    @assert !isnan(分数)
    return 分数
end

function summarize_test(涨幅, 时间戳, 代码, 最新价, 买1价, 卖1价, 手续费率, 交易池, 记忆仓位, 实际仓位, 综合评分, 虚拟信号, 最大持仓)
    收益率, 资金曲线, 年化收益率, 最大回撤, 夏普率 = summarize_core(涨幅, 时间戳, 代码, 最新价, 买1价, 卖1价, 手续费率, 交易池, 记忆仓位, 实际仓位, 综合评分, 最大持仓)
    目录名 = mkpath(产生目录(时间戳, 代码, 年化收益率, 最大回撤, 夏普率))
    @indir 目录名 begin
        输出每日持股明细(代码, 时间戳, 实际仓位)
        输出资金和仓位曲线(时间戳, 实际仓位, 资金曲线)
        输出个股交易记录(时间戳, 代码, 最新价, 买1价, 卖1价, 手续费率, 收益率, 实际仓位)
        新目录名 = 输出资金曲线(时间戳, 代码, 实际仓位, 收益率, 最大持仓)
        报告后处理()
        parseenv("OUTPUT_POS", false) && 输出仓位评分信号(代码, 时间戳, 实际仓位, 综合评分, 虚拟信号)
        parseenv("OUTPUT_MINPOS", false) && 输出分钟仓位(代码, 时间戳, 最新价, 实际仓位, 收益率)
    end
    return 年化收益率
end


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

using Markdown
using InteractiveUtils

# ╔═╡ bca1c92c-8561-11eb-0420-6f0123bacb7f
begin
	# dependencies
	using Plots
	using LightGraphs, SimpleWeightedGraphs, LinearAlgebra, Random
	import StatsBase:sample
end

# ╔═╡ 774fa124-8560-11eb-2317-53b0d106b67a
md"""
# Local Search
"""

# ╔═╡ 2b5e6876-8561-11eb-0685-dfe1cca39c45
md"""
## Simulated annealing
Simulated annealing is similar to hill climbing, but it allows downhill moves. There are several aspects to consider.
* What is considered a neighbor? And how to chose them?
* How to evaluate the energy (i.e. the fitness)?
* How to let the temperature evolve. Several options exist here:
    - exponential decrease: $T = T_0(1-\alpha)^k$, where $\alpha$ is the ‘cooling rate’ and $k$ is the iteration
    - fast decrease: $T = T_0/k$
    - boltzmann decrease: $T = T_0/log(k)$
"""

# ╔═╡ 4e190a2e-8561-11eb-3759-795e3bb4d8d2
md"""
To illustrate the principle, let's try to find the maximum of
\\[ f(x) = -\dfrac{1}{4}x^4 +  \dfrac{7}{3} x^3 - 7x^2 + 8x \\]
in the interval $[0, 5]$. You can observe that there is a global maximum in $x=4$, however, there is a local maximum in $x=1$.
"""

# ╔═╡ b49f0528-8561-11eb-36a7-fd42574edc81
begin
	f = x-> -1/4*x^4 + 7/3 *x^3 - 7*x^2 + 8*x
	x = range(0,5, length=100);
	plot(x,f.(x),label="",size=(300,200))
end

# ╔═╡ d5019d8a-8561-11eb-2a43-a9ce616532f5
begin
	"""
	Evaluate the fitness of a specific point
	"""
	function fitness(x::Float64,f::Function)
		return f(x)
	end

	"""
	Determine successors for a given location (with a fixed step size)
	You could include automatic step size modification if you want.
	"""
	function successors(x::Float64; step=0.1)
		return x .+ step*[-1; 1]
	end

	"""
	mapping of time to temperature. You might make this a wrapper function for
	different cooling schemes. Currently only has exponential cooling.
	"""
	function scheduletemp(t::Int64;T0::Float64=273.0, α::Float64=0.005)
		return T0*(1-α)^t
	end

	"""
		simulated_annealing(args; kwargs)

	Simulated annealing implementation for a 1D problem.

	Since this is stochastic (sort of a random walk), it will not always reach the optimal value, therefore, one typically runs several iterations and the best one is retained.

	The probabiliy of accepting a value that lead to a decrease in energy is given by exp(ΔE/T), which is between 0 (temp = 0) and one (ΔE for any T ≠ 0).

	In the literature you might also find  1/(1 + exp(ΔE/T)) ∈ [1/2, 1] as a treshold.

	arguments:
	- x0: starting point
	keywords:
	- tmax: maximum iterations
	- T0: initial temperature
	- α: cooling rate
	- fitness: fitness function
	- optimfun: the function we want to optimize
	"""
	function simulated_annealing(x0;
				tmax=1000, T0::Float64=273.0, α::Float64=0.005,
				fitness::Function=fitness, optimfun::Function=f)
		current = x0
		next = x0
		T = T0
		best = x0
		for t = 1:tmax
			T = scheduletemp(t,T0=T0,α=α)
			if T == 0
				return current
			else
				next = rand(successors(current))
				ΔE = fitness(next, optimfun) - fitness(current, optimfun)
				if ΔE > 0
					current = next
				else
					if rand() < exp(ΔE/T)
						current = next
					end
				end
				if fitness(current, optimfun) > fitness(best, optimfun)
					best = current
				end
			end
		end
		return best
	end

	"""
		SAsearch

	Wrapper function for multiple SA iterations. Returns the best n found value for N runs of simulated annealing.
	"""
	function SAsearch(x0, n::Int64=1, N::Int64=10;kwargs...)
		return sort!([simulated_annealing(x0; kwargs...) for _ in 1:N],
					 by=x->kwargs[:fitness](x,kwargs[:optimfun]),
					 rev=true)[1:n]
	end
end

# ╔═╡ 9c2e2032-8568-11eb-103f-fdc57c11ab1f
md"""
### Solving the TSP problem with simulated annealing
Use simulated annealing to solve the travelling saleman problem again. You will need to think in particalar about:
* an appropriate inital solution (e.g. a tour by using a nearest neigbor approach or purely random)
* generating appropriate neigbors (e.g. by flipping a subpart of a tour)
"""

# ╔═╡ c84f6af2-8568-11eb-06c4-53f7664294ac
begin
	"""
		plotgraph(g::SimpleWeightedGraph, d::Dict,draw_MST::Bool=true,wf::Real=5)

	Visualize the network with its edges and their weights

	parameters:
		- g : SimpleWeightedGraph
		- d : Dict
			dictionary that links a node to a coordinate pair
		- draw_MST : Bool
			if true, also draws the MST on the figure
		- wf : Real
			multiplication factor for the edge weigth to obtain a nicer visualisation
	"""
	function plotgraph(g::SimpleWeightedGraph, d::Dict;wf::Real=5,draw_MST::Bool=true)
		p = plot(legend=:bottomright)
		for e in edges(g)
			plot!([d[e.src][1], d[e.dst][1]],[d[e.src][2], d[e.dst][2]],
				markershape=:circle, linewidth=e.weight*wf, color=:black, markeralpha = 0.5,
				linealpha=0.2, label="", markercolor=:black, markersize=5)
		end
		plot!([],[],markershape=:circle, color=:black,
				linealpha=0.2, label="Full network", markercolor=:black, markersize=5, markeralpha = 0.5)
		if draw_MST
			MST = prim_mst(g, weights(g))
			for e in MST
				plot!([d[e.src][1], d[e.dst][1]],[d[e.src][2], d[e.dst][2]],
						markershape=:circle, linewidth=weights(g)[e.src,e.dst]*wf, color=:red,
						linealpha=0.8, label="", linestyle=:solid, markeralpha = 0)

			end
			plot!([],[],linealpha=0.8, label="MST", color=:red, linestyle=:solid)
		end
		xlims!(0,1)
		ylims!(0,1)
		plot!(legend=:best)
		return p
	end

	"""
		plotsol(fp::Array{Int64,1}, nodedict::Dict)

	Function that plots the final TSP solution
	"""
	function plotsol(fp::Array{Int64,1}, nodedict::Dict)
		p = plot()
		for i in 1:length(fp)-1
			x = [nodedict[fp[i]][1]; nodedict[fp[i+1]][1]]
			y = [nodedict[fp[i]][2]; nodedict[fp[i+1]][2]]
			plot!(x,y,color=:blue,label="", marker=:circle)
		end
		# close loop
		x = [nodedict[fp[1]][1]; nodedict[fp[end]][1]]
		y = [nodedict[fp[1]][2]; nodedict[fp[end]][2]]
		plot!(x,y,color=:red,label="", marker=:circle)
		xlims!(0,1)
		ylims!(0,1)
		return p
	end

	"plot functions"
end

# ╔═╡ 46e468e0-8569-11eb-0a30-d11e814ccf7e
begin
	"""
		problemgenerator(n::Int)

	generate n points in the unit square
	"""
	function problemgenerator(n::Int)
		return [(rand(), rand()) for _ in 1:n]
	end

	"""
		distancematrix(x::Array{Tuple{Float64,Float64},1}, p::Real=1)

	Generate distancematrix between all points.

	x is an array [n x 1] and p is the L_p-norm you want to use for the distancematrix
	"""
	function distancematrix(x::Array{Tuple{Float64,Float64},1}, p::Real=1)
		n = length(x)
		A = zeros(n,n)
		for (i,j) in collect(Base.Iterators.product(1:n, 1:n))
			A[i,j] = norm(x[i] .- x[j], p)
		end
		return A
	end
end

# ╔═╡ 7c5282fa-8569-11eb-0dcf-abc7b3b795ef
begin
	N = problemgenerator(15) # nodes
	nodedict = Dict( i => N[i] for i in 1:length(N)) # nodes dict
	A = distancematrix(N,2) # distance matrix using manhattan norm
	G = SimpleWeightedGraph(A) # graph
	p = plotgraph(G,nodedict;draw_MST=false) # layout of the current graph
end

# ╔═╡ eb27eada-8569-11eb-0739-7717087a5240
begin
	"""
		finit(A::Array{Float64,2})

	Generate an initial distribution based and a random order
	"""
	function finit(A::Array{Float64,2})
		return shuffle(1:size(A,1))
	end

	finit(A)
end

# ╔═╡ 1692c654-856a-11eb-1b02-77b0ab027c29
begin
	"""
		fitness(T::Array{Int64,1},A::Array{Float64,2})

	Given a tour, compute the tour distance based on the distance matrix
	"""
	function fitness(T::Array{Int64,1},A::Array{Float64,2})
		return -sum([A[T[i],T[i+1]] for i in 1:length(T)-1]) - A[T[1],T[end]]
	end

	fitness(finit(A),A)
end

# ╔═╡ 3c58cc1a-856a-11eb-17bd-a31c244f6229
begin
	"""
	Determine successors for a given tour (simple random permutation of partial tour)
	"""
	function successors(T::Array{Int64,1}; nsuc::Int64=5)
		res = Array{Array{Int64,1},1}()
		for _ in 1:nsuc
			i,j = sample(1:length(T),2,replace=false,ordered=true)
			cand = copy(T)
			cand[i:j] = shuffle(T[i:j])
			push!(res, cand)
		end

		return res
	end
	successors(finit(A))
end

# ╔═╡ 7ab96168-856c-11eb-166a-b7531078ac90
begin
	"""
		simulated_annealing(args; kwargs...)

	Simulated annealing implementation for TSP

	Since this is stochastic (sort of a random walk), it will not always reach the optimal value, therefore, one typically runs several iterations and the best one is retained.

	The probabiliy of accepting a value that lead to a decrease in energy is given by exp(ΔE/T). In the literature you might also find  1/(1 + exp(ΔE/T)) ∈ [1/2, 1] as a treshold.
	"""
	function simulated_annealing(A::Array{Float64,2}, init::Array{Int64,1};
								tmax=1000, T0::Float64=273.0, α::Float64=0.005,
								fitness::Function=fitness,nsuc::Int64=10, kwargs...)
		current = init
		next = init
		T = T0
		best = init
		for t = 1:tmax
			T = scheduletemp(t,T0=T0)
			if T == 0
				return current
			else
				next = rand(successors(current, nsuc=nsuc))
				ΔE = fitness(next, A) - fitness(current, A)
				if ΔE > 0
					current = next
				else
					if rand() < exp(ΔE/T)
						current = next
					end
				end
				if fitness(current, A) > fitness(best, A)
					best = current
				end
			end
		end
		return best
	end

	"""
	wrapper function for multiple SA iterations

	Returns the best n found value for N runs of simulated annealing.
	"""
	function SAsearch(A::Array{Float64,2}, n::Int64=1, N::Int64=10;kwargs...)
		sort!([simulated_annealing(A, finit(A); kwargs...) for _ in 1:N],
				by=x->kwargs[:fitness](x,A),rev=false)[1:n]
	end

	simulated_annealing(A, finit(A))
end

# ╔═╡ 934aa9e4-856e-11eb-1924-fd1df8fd3fe2
md"""
## Genetic algorithms
a genetic algorithm (GA) is inspired by the process of natural selection. In order to implement this you need the following items:
* an initial population (e.g. random)
* a fitness function that allows you to rank the population members
* a selection method (e.g. roulette based proportions, ‘elite selection’)
* a cross-over method (cutting; 50/50 per gene etc.)
* mutation rate and method
* \\(\dots\\)
"""

# ╔═╡ 8753fd9a-856f-11eb-2d91-13657e302031
md"""
### Toy example
Let's try write a message by using a genetic algorithm. E.g "172 Pol evolved!" starting from a random set of characters (a-zA-Z0-9. !)
"""

# ╔═╡ 48c05e06-8570-11eb-0367-1f3db0aa20d0
begin
	chars = vcat(65:90,97:122, 48:57,[Int(' ');Int('.'); Int('!')])
	const genes = [Char(i) for i = chars]
	const pop_size = 100
	const goal = collect("172 Pol evolved!");
end

# ╔═╡ 6faf987e-8570-11eb-1861-7de0f661bd63
begin
	function fitness(p::Array{Char,1}; goal=goal)
		sum(p.==goal)
	end

	testcase = rand(genes,length(goal))
	testcase, fitness(testcase)
end

# ╔═╡ 728b94da-8566-11eb-1b02-e94c64424845
begin
	options = Dict(	:tmax=>1000, :T0=>273.0, :α=>0.005,
               		:fitness=>fitness, :optimfun=>f)
	x0 = 0.0;
	simulated_annealing(x0; options...)
end

# ╔═╡ a72d7c3c-8569-11eb-36a3-cb68ab2622cf
let
	n = 10;
	res = SAsearch(x0,n;options...)
end

# ╔═╡ 0ce6432c-8567-11eb-1465-432d66051360
begin
	n = 10;
	res = SAsearch(x0,n;options...)
	plot(x,f.(x),label="")
	scatter!([x0],[f.(x0)],label="start",marker=(:circle, :green),size=(300,200), legend=:bottom)
	scatter!([res],[f.(res)],label="optimal values",marker=(:circle, :red),size=(300,200), legend=:bottom)
end

# ╔═╡ ef64ca86-856c-11eb-0cd0-178565e85576
let
	options = Dict(:tmax=>1000, :T0=>273.0, :α=>0.005, :fitness=>fitness)
	@show res = SAsearch(A; options...)[1];
	init = finit(A)
	println("Cost of a random trip: $(-fitness(init,A))")
	println("Cost of a optimised trip: $(-fitness(res,A))")
	plotsol(res, nodedict)
	title!("TSP solution using simulated annealing")
end

# ╔═╡ 57726a64-856d-11eb-2b37-877202c47ede
let
	# Larger example
	n = 50 # number of nodes
	N = problemgenerator(n) # nodes
	options = Dict(:tmax=>1000, :T0=>273.0, :α=>0.005, :fitness=>fitness)
	nodedict = Dict( i => N[i] for i in 1:length(N)) # nodes dict
	A = distancematrix(N,2) # distance matrix using manhattan norm
	init = finit(A)

	@time res = SAsearch(A; options...)[1]
	println("Cost of a random trip: $(-fitness(init,A))")
	println("Cost of a optimised trip: $(-fitness(res,A))")
	plotsol(res, nodedict)
end

# ╔═╡ d70478c8-856d-11eb-3f27-0b565955b6ed
let
	# trajectory cost distribution
	res = Array{Float64,1}()
	for _ in 1:50
		push!(res, -fitness(SAsearch(A;options...)[1], A))
	end
	histogram(res, normalize=:pdf, alpha=0.5,label="PDF optimised tour cost")
	scatter!([-fitness(finit(A),A)],[0],label="initial non-optimised cost")
	xlabel!("Tour cost")
	ylabel!("PDF")
	title!("Histogram for tour cost ($(size(A,1)) cities)")
end

# ╔═╡ 9cc2fa40-8570-11eb-2eba-c15433dcb90e
begin
	"""
		mate(p1,p2; pmut, genes)

	mating function to obtain children. includes both crossover and mutation.
	"""
	function mate(p1::Array{Char,1}, p2::Array{Char,1}; pmut=0.1, genes=genes)
		child = typeof(p1)() # empty child
		# walk over genome sequence
		for i in 1:length(p1)
			p = rand()
			if p < (1-pmut)/2
				# parent 1 wins
				push!(child, p1[i])
			elseif p < 1-pmut
				# parent 2 wins
				push!(child, p2[i])
			else
				# mutation
				push!(child, rand(genes))
			end
		end
		return child
	end

	# demo
	p1 = ['a';'b';'c']
	p2 = ['d';'e';'f']
	mate(p1,p2, pmut=1)
end

# ╔═╡ 47153f8e-8572-11eb-1310-c932bd478721
md"""
### Some ideas for crossover:
Below you can find some algorithms that have proven to work quite well for this type of application. Both of these are described and an example and an implementation is provided. You are also invited to come up with a method yourself if you think you have a good idea.

#### The partially mapped crossover
After choosing two random cut points on parents to build offspring, the portion between cut points, one parent’s string is mapped onto the other parent’s string and the remaining information is exchanged.

E.g. for parents $P_1$ and $P_2$ with cutting positions $c_1$ and $c_2$:
* ``P_1 = [3, 4, 8 | 2, 7, 1 | 6, 5 ]``
* ``P_` = [4, 2, 5 | 1, 6, 8 | 3, 7 ]``
We get the following mappings:
* ``m_{1 \rightarrow 2}: \{2 \rightarrow 1, 7 \rightarrow 6, 1 \rightarrow 8\}``
* ``m_{2 \rightarrow 1}: \{1 \rightarrow 2, 6 \rightarrow 7, 8 \rightarrow 1\}``
which leads to the following initial offspring:
* ``O_1 = [0, 0, 0 | 1, 6, 8 | 0, 0 ] ``
* ``O_2 = [0, 0, 0 | 2, 7, 1 | 0, 0 ] ``
we then add the non-conflicting values (i.e. the values not already present the offspring due do the crossover) from the original parents:
* ``O_1 = [3, 4, 0 | 1, 6, 8 | 0, 5 ]``
* ``O_2 = [4, 0, 5 | 2, 7, 1 | 3, 0 ] ``
For the remaining values, we (recursively) make use of the mappings using $m_{2 \rightarrow 1}$ for $P_1$ and $m_{1 \rightarrow 2}$ for $P_2$. This finally leads to the following offspring:
* ``O_1 = [3, 4, 2 | 1, 6, 8 | 7, 5 ] ``
* ``O_2 = [4, 8, 5 | 2, 7, 1 | 3, 6 ] ``

#### The cycle crossover
The Cycle Crossover operator identifies a number of so-called cycles between two parent chromosomes. To form Child 1, cycle 1 is copied from parent 1, cycle 2 from parent 2, cycle 3 again from parent 1, and so on. An example might help to better understand the principle of a cycle. E.g.
* ``P_1 = [1,2,3,4,5,6,7,8] ``
* ``P_2 = [8,5,2,1,3,6,4,7] ``
For the first cycle, we start at index 1 in $P_1$. Following the links from $P_1$ to $P_2$ we find the following mappings:
* ``1 \rightarrow 8`` with index 8
* ``8 \rightarrow 7`` with index 7
* ``7 \rightarrow 4`` with index 4
* ``4 \rightarrow 1`` with index 1
after this last mapping, we are back at the starting point, so our first cycle is complete and is $\{1,4,7,8 \}$.

The second cycle will start at the first index that is not yet present in a cycle. In this case the index is equal to 2. Following the links from $P_1$ to $P_2$ we find the following mappings
* ``2 \rightarrow 5`` with index 2
* ``5 \rightarrow 3`` with index 5
* ``3 \rightarrow 2`` with index 3

So our second cycle is equal to $\{2,3,5\}$.

The third an final cycle is equal to the only remaining part, i.e. $\{6 \}$

We have now determined all cycles. The children aren obtained as follows:
* ``O_1[1,4,7,8] =  P_1[1,4,7,8] ``
* ``O_1[2,3,5] = P_2[2,3,5] ``
* ``O_1[6]  = P_1[6] ``

This leads to $O_1 = [1,5,2,4,3,6,7,8]$. In a similar fashion, but inverting the role of $P_1$ and $P_2$ one finds $O_2 = [8,2,3,1,5,6,4,7]$.

A downside of this type of crossover is that it can sometimes lead to the offspring being identical to the parents. You can verify this for $[3,4,8,2,7,1,6,5]$ & $[4,2,5,1,6,8,3,7]$.

#### Your own method
Come up with something that seems adequate...

"""

# ╔═╡ 1f5d255c-8576-11eb-1463-89f94258168d
begin
	"""
		mate(p1,p2;method=partialcrossover, kwargs...)

	Function that handles the crossover for two different tours. The method keyword argument can be any of the functions you define as long as 1. the function exists 2. the function works in a similar way as the other ones (i.e. return two children from two parents.)

	Additional options for your own crossover implementation can be passed via the kwargs. You should document the options for each algorithm sperately.

	"""
	function mate(p1::Array{Int64,1},p2::Array{Int64,1}; method::Function=partialcrossover, kwargs...)
			return method(p1,p2; kwargs...)
	end

	"""
		partialcrossover(p1, p2; kwargs...)

	Implementation of the partial crossover function for TSP routing problems.
	"""
	function partialcrossover(p1::Array{Int64,1},p2::Array{Int64,1}; kwargs...)
		"""
			deconflicter!(O::Array{Int64,1}, mapping::Dict)

		Function that deals with the conflicts after the initial assignments.
		"""
		function deconflicter(O::Array{Int64,1},p::Array{Int64,1}, mapping::Dict)
			for i in 1:length(O)
				if O[i] == 0
					cand = p[i]
					while cand in O
						cand = mapping[cand]
					end
					O[i] = cand
				end
			end
			return O
		end

		# choose cut points
		c1, c2 = sort(sample(1:length(p1),2, replace=false))
		# initiate offspring
		O1, O2 = [zeros(Int64,length(p1)) for _ in 1:2]
		# determine mapping
		map12 = Dict(p1[i]=> p2[i] for i in c1:c2)
		map21 = Dict(value => key for (key, value) in map12)
		# set offspring initial values
		O1[c1:c2] = p2[c1:c2]
		O2[c1:c2] = p1[c1:c2]
		# fill non-conflicting positions
		for i in setdiff(1:length(p1),c1:c2)
			O1[i] = p1[i] in O1 ? 0 : p1[i]
			O2[i] = p2[i] in O2 ? 0 : p2[i]
		end
		# fill conflicting positions using the mapping
		O1, O2 = deconflicter(O1, p1, map21), deconflicter(O2, p2, map12)

		return O1, O2
	end

	"""
		cyclecrossover(p1, p2)

	Implementation of the cycle crossover function for TSP routing problems.
	"""
	function cyclecrossover(p1::Array{Int64,1},p2::Array{Int64,1})
		# initiate children
		o1=zeros(Int64,length(p1))
		o2=zeros(Int64,length(p2))

		# Pt. 1 - find all cycles

		# make hashtable: value => index
		h1 = Dict{Int64, Int64}(p1[i]=>i for i in 1:length(p1))
		h2 = Dict{Int64, Int64}(p2[i]=>i for i in 1:length(p2))
		# make hashtable to track usage in a cycle: index => bool
		u1 = Dict{Int64, Bool}(i=>false for i in 1:length(p1))
		# loop over parent creating cycles on the go
		# a cycle is stored in a Set because this is hashable. This is (much) faster than an array for ∉ check
		# the values stores in the cycles are the indices to be used
		cycles = Array{Set{Int64},1}()
		for i in 1:length(p1)
			if !u1[h1[p1[i]]] # if not used yet, start a new cycle
				cycle = Set{Int64}()
				next = h1[p1[i]] # add initial values
				while next ∉ cycle
					push!(cycle, next)  # update cycle
					u1[next] = true     # update tracking
					next = h1[p2[next]] # next location
				end
				push!(cycles, cycle)
			else
				continue
			end
		end

		# Pt. 2 - make children from cycles
		for i in 1:length(cycles)
			inds = collect(cycles[i])
			if isodd(i)
				o1[inds] = p1[inds]
				o2[inds] = p2[inds]
			else
				o1[inds] = p2[inds]
				o2[inds] = p1[inds]
			end
		end
		return o1, o2
	end

	"""
		mycrossover(p1, p2; myoption::Float64=0.5, kwargs...)

	Your own method that takes predefined options and also accepts others.
	"""
	function mycrossover(p1::Array{Int64,1},p2::Array{Int64,1}; myoption::Float64=0.5, kwargs...)
		error("Function <mycrossover> not defined yet... Get creative!")
	end
end

# ╔═╡ fd95d464-8570-11eb-1a38-879d8079d127
begin
	function textga(;goal::Array{Char,1}=goal, popsize::Int64=pop_size,
                	 ngen::Int64=500, next::Float64=0.4, pmut=0.1, kwargs...)
		# initiate population
		population = [rand(genes, length(goal)) for _ in 1:popsize]
		gen = 0
		for _ in 1:ngen
			# sort by fitness
			sort!(population; by=x->fitness(x), rev=true)
			if population[1] == goal
				break
			end
			gen += 1
			# select adequate parents
			goodparents = population[1: round(Int,popsize*next)]
			# make children
			new_generation = typeof(population)()
			for _ in 1:popsize
				push!(new_generation, mate(rand(goodparents), rand(goodparents),
											pmut=pmut ))
			end
			population = new_generation
		end
		@info "stopped after $(gen) generations"

		return (gen, prod(population[1]), fitness(population[1])/fitness(goal))
	end

	textga(ngen=500, popsize=100, pmut=0.1)
end

# ╔═╡ 7c4ef658-8576-11eb-3f6d-d3587a1af7e8
begin
	let
		# small demo matching the explanation of the algorithm above.
		# Matches the detailed example if you force c1, c2 = 4, 6 in the partialcrossover algorithm
		p1 = [3,4,8,2,7,1,6,5]
		p2 = [4,2,5,1,6,8,3,7]
		@time mate(p1,p2)
	end

	let
		# Longer route does not necessarily imply greater time for offspring generation
		# with the partial crossover algorithm (why could that be?)
		@info "100 items"
		N = 100
		p1 = collect(1:N)
		p2 = sample(p1,N,replace=false)
		@time mate(p1,p2);
	end

	let
		# demo from the method description
		p1 = [1,2,3,4,5,6,7,8]
		p2 = [8,5,2,1,3,6,4,7]
		@time res = mate(p1,p2, method=cyclecrossover)
		@assert res ==  ([1,5,2,4,3,6,7,8],  [8,2,3,1,5,6,4,7])
	end

	let
		# demo for children equal to parents
		p1 = [3,4,8,2,7,1,6,5]
		p2 = [4,2,5,1,6,8,3,7]
		@time res = mate(p1,p2, method=cyclecrossover)
		@assert res == (p1 ,p2)
	end

	let
		# Longer routes scales nicely for cyclecrossover as well (≈10x longer route and only ≈2x slower)
		N = 100
		p1 = collect(1:N)
		p2 = sample(p1,N,replace=false)
		@time res = mate(p1,p2, method=cyclecrossover);
	end
end

# ╔═╡ 8b77cf1c-8571-11eb-2ea7-531abac106cc
md"""
### Genetic algorithm - TSP
Now that you are familiar with a genetic algorithm, let's try to solve the TSP problem again, but this time, using a genetic algorithm.

We should think about:
* representation
* an initial configuration
* a fitness function
* a selection method (e.g. roulette based proportions, ‘elite selection’ etc.)
* cross-over method: (what is the difference with the text application?)
* mutation rate and method (is this required? why (not)?)
* when do you stop making new generations?
* \\( \dots \\)

Tasks:
* Come up with your own suggestions for cross-over/mutation
* Implement the genetic algorithm to solve the TSP
* Compare the performance of the tree cross-over algorithms to deal with eachother. Give that there is some randomness involved, a statistical approached might be more appropriate. Since the initial population will most likely be chosen at random, you should impose the random seed to make sure you have the same starting point for the comparisons.
* Compare the scalability with what you know from previous practical sessions (tree search etc.) What is better and what is worse?
"""

# ╔═╡ 712cb76c-8576-11eb-37c3-87beed1c3fa1
begin
	function TSP_ga(d::Array{Float64,2};
					popsize::Int64=30, ngen::Int64=500, next::Float64=0.4,
					method::Function=partialcrossover, kwargs...)
		# initiate population
		n = size(d,1)
		population = [sample(1:n,n,replace=false) for _ in 1:popsize]
		gen = 0
		# go over the generations
		for _ in 1:ngen
			# sort by fitness
			sort!(population; by=x->fitness(x,d), rev=true)
			#@show "current tour cost: $(fitness(population[1],d))"
			# include a stopcondition?
			gen += 1
			# select adequate parents
			goodparents = population[1: round(Int,popsize*next)]
			# make children
			new_generation = typeof(population)()
			for _ in 1:2:popsize
				p1,p2 = rand(goodparents,2)
				offspring = mate(p1,p2;method=method, kwargs...)
				append!(new_generation, offspring)
			end
			population = new_generation
		end

		return (gen, population[1], fitness(population[1],d))
	end
end

# ╔═╡ fa780328-8576-11eb-3a16-b7109ffbad63
let
	n = 50 # number of nodes
	N = problemgenerator(n) # nodes
	nodedict = Dict( i => N[i] for i in 1:length(N)) # nodes dict
	A = distancematrix(N,1) # distance matrix using manhattan norm
	@info fitness(finit(A), A)
	G = SimpleWeightedGraph(A) # graph
	@time begin
    	D = Dict()
    	methods = [cyclecrossover, partialcrossover]
		for method in methods
			res = []
			for _ in 1:100
				(_, tour, cost) = TSP_ga(A, ngen=20, popsize=n,method=method)
				push!(res,(tour,cost))
			end
			D[method] = [x[2] for x in res]
		end
	end

	# illustration
	histogram()
	for (key, val) in D
		histogram!(-val, label=string(key), normalize=:pdf, alpha=0.5)
	end
	scatter!([-fitness(finit(A),A)],[0],label="initial non-optimised cost")
	xlabel!("Tour cost")
	ylabel!("PDF")
	title!("Histogram for tour cost ($(size(A,1)) cities)")
end

# ╔═╡ 1a79eb5a-8577-11eb-168e-8fa9230f26fb

# ╔═╡ Cell order:
# ╟─774fa124-8560-11eb-2317-53b0d106b67a
# ╠═bca1c92c-8561-11eb-0420-6f0123bacb7f
# ╠═2b5e6876-8561-11eb-0685-dfe1cca39c45
# ╠═4e190a2e-8561-11eb-3759-795e3bb4d8d2
# ╠═b49f0528-8561-11eb-36a7-fd42574edc81
# ╠═d5019d8a-8561-11eb-2a43-a9ce616532f5
# ╠═728b94da-8566-11eb-1b02-e94c64424845
# ╠═a72d7c3c-8569-11eb-36a3-cb68ab2622cf
# ╠═0ce6432c-8567-11eb-1465-432d66051360
# ╟─9c2e2032-8568-11eb-103f-fdc57c11ab1f
# ╟─c84f6af2-8568-11eb-06c4-53f7664294ac
# ╠═46e468e0-8569-11eb-0a30-d11e814ccf7e
# ╠═7c5282fa-8569-11eb-0dcf-abc7b3b795ef
# ╠═eb27eada-8569-11eb-0739-7717087a5240
# ╠═1692c654-856a-11eb-1b02-77b0ab027c29
# ╠═3c58cc1a-856a-11eb-17bd-a31c244f6229
# ╠═7ab96168-856c-11eb-166a-b7531078ac90
# ╠═ef64ca86-856c-11eb-0cd0-178565e85576
# ╠═57726a64-856d-11eb-2b37-877202c47ede
# ╠═d70478c8-856d-11eb-3f27-0b565955b6ed
# ╟─934aa9e4-856e-11eb-1924-fd1df8fd3fe2
# ╟─8753fd9a-856f-11eb-2d91-13657e302031
# ╠═48c05e06-8570-11eb-0367-1f3db0aa20d0
# ╠═6faf987e-8570-11eb-1861-7de0f661bd63
# ╠═9cc2fa40-8570-11eb-2eba-c15433dcb90e
# ╠═fd95d464-8570-11eb-1a38-879d8079d127
# ╟─47153f8e-8572-11eb-1310-c932bd478721
# ╠═1f5d255c-8576-11eb-1463-89f94258168d
# ╠═7c4ef658-8576-11eb-3f6d-d3587a1af7e8
# ╟─8b77cf1c-8571-11eb-2ea7-531abac106cc
# ╠═712cb76c-8576-11eb-37c3-87beed1c3fa1
# ╠═fa780328-8576-11eb-3a16-b7109ffbad63
# ╠═1a79eb5a-8577-11eb-168e-8fa9230f26fb

max_stars_repo_path: PS/PS05 - Local Search.jl

max_issues_repo_path: PS/PS05 - Local Search.jl

max_forks_repo_path: PS/PS05 - Local Search.jl