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

using Markdown
using InteractiveUtils

# ╔═╡ 73112674-61f9-4915-9b09-216c70df0e74
begin
	using MortalityTables

	vbt = MortalityTables.table("2001 VBT Residual Standard Select and Ultimate - Male Nonsmoker, ANB") #or any other table
end

# ╔═╡ ce999823-acab-47bb-bb15-babb8018a06a
using DataFrames

# ╔═╡ a5611f56-cca3-4a06-842e-715a89d401e4
begin
	# add a table of contents to the page
	using PlutoUI
	TableOfContents()
end

# ╔═╡ 059c84bc-5d87-4781-b0f2-15563bfae09d
md"# Using MortaltiyTables.jl with DataFrames

## MortalityTables.jl standard representation

MortalityTables.jl stores the rates in a very efficient manner as a collection of vectors indexed by attained age.
"

# ╔═╡ 6e38cf1f-0c8e-47b9-baa1-e9525f0a1faa
md"First, we include the package, and then we'll pick a table, where all of the `mort.soa.org` tables are mirrored into your MortalityTables.jl installation."

# ╔═╡ c9770fc4-e460-46d9-ba6b-2ab0fc17830e
md" To see how the data is represented, we can look at the the select data for a 55 year old and see the attained age and mortality rates:"

# ╔═╡ 6bce5873-d064-46ea-9969-7809051300db
vbt.select[55]

# ╔═╡ 5f0e15fc-1107-11ec-2a58-3351c91877fc
md" This is very efficient and convienent for modeling, but a lot of times you want the data matched up with policy data in a DataFrame."

# ╔═╡ f0b7e725-b433-458a-a5bf-efd5dbc990ce
md"## Getting data into a dataframe"

# ╔═╡ b5a00124-c994-4b79-99dc-6b966bec88cd
md"### Generate sample data"

# ╔═╡ 2d497d32-b088-4027-ad5b-0c4597db3922
sample_size = 10_000

# ╔═╡ dbb95d40-6517-491b-b6a9-abccbd2a3a78
sample_data = let
	# generate fake data
	df = DataFrame(
		"sex" => rand(["Male","Female"],sample_size),
		"smoke" => rand(["Smoker","Nonsmoker"],sample_size),
		"issue_age" => rand(25:65,sample_size),
		)

	# a random offset of issue age is the current attained age
	df.attained_age = df.issue_age .+ rand(1:10,sample_size)
	df
end

# ╔═╡ e0acb2d4-668e-40aa-ac70-de4527c7a891
md"### Define the table set you want to use"

# ╔═╡ 423f35bd-6b79-464d-a8f6-a2405e5ba3ce
md"
There are a lot of different possible combinations of parameters that you might want to use, such as rates that vary by sex, risk class, table set (VBT/CSO/etc), smoking status, relative risk, ALB/ANB, etc.

It's easy to define the parameters applicable to your assumption set. Here, we'll use a dictionary to define the relationship:
"

# ╔═╡ c351704b-4bf8-4fc0-bda2-c2b1970c2a11
rate_map = Dict(
	"Male" => Dict(
		"Smoker" => MortalityTables.table("2001 VBT Residual Standard Select and Ultimate - Male Smoker, ANB"),
		"Nonsmoker" => MortalityTables.table("2001 VBT Residual Standard Select and Ultimate - Male Nonsmoker, ANB"),
		),

	"Female" => Dict(
		"Smoker" => MortalityTables.table("2001 VBT Residual Standard Select and Ultimate - Female Smoker, ANB"),
		"Nonsmoker" => MortalityTables.table("2001 VBT Residual Standard Select and Ultimate - Female Nonsmoker, ANB"),
		)
	)

# ╔═╡ 4d88a451-f350-4a72-ae5b-8698b0527511
md" and then we'll define a function to look up the relevant rate. Note how the function matches the levels we defined for the assumption set dictionary above."

# ╔═╡ 45f4378e-9299-429d-b361-8f99b2be6efd
function rate_lookup(assumption_map,sex,smoke,issue_age,attained_age)
	# pick the relevant table
	table = assumption_map[sex][smoke]

	# check if the select rate exists, otherwise look to the ulitmate table
	if issue_age in eachindex(table.select)
		table.select[issue_age][attained_age]
	else
		table.ultimate[attained_age]
	end
end

# ╔═╡ 6d5b5548-7cbc-435e-8058-ebd529f839f1
md"### Lining up with dataframe

By mapping each row's data to the lookup function, we get a vector of rates for our data:"

# ╔═╡ 372694fa-60a8-4305-a922-f020a4deb896
rates = map(eachrow(sample_data)) do row
	rate_lookup(rate_map,row.sex,row.smoke,row.issue_age,row.attained_age)
end

# ╔═╡ a8cfd44f-5e40-4624-b98e-00de4a623af2
md"And finally, we can just add this to the dataframe:"

# ╔═╡ 7fe85397-51a5-4d82-86ee-f6881dcf8808
sample_data.expectation = rates;

# ╔═╡ 47da227d-72ab-4eae-bb27-3a61d842bb03
sample_data

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
MortalityTables = "4780e19d-04b9-53dc-86c2-9e9aa59b5a12"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
DataFrames = "~1.2.2"
MortalityTables = "~2.1.4"
PlutoUI = "~0.7.9"
"""

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

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

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

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

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

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

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

[[DataAPI]]
git-tree-sha1 = "bec2532f8adb82005476c141ec23e921fc20971b"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.8.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"

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

[[EzXML]]
deps = ["Printf", "XML2_jll"]
git-tree-sha1 = "0fa3b52a04a4e210aeb1626def9c90df3ae65268"
uuid = "8f5d6c58-4d21-5cfd-889c-e3ad7ee6a615"
version = "1.1.0"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[MortalityTables]]
deps = ["OffsetArrays", "Parsers", "Pkg", "QuadGK", "Requires", "UnPack", "XMLDict"]
git-tree-sha1 = "cae9281196e40e24cb5885b8463a6f955ae0f3eb"
uuid = "4780e19d-04b9-53dc-86c2-9e9aa59b5a12"
version = "2.1.4"

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

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

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

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

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

[[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 = "a193d6ad9c45ada72c14b731a318bedd3c2f00cf"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.3.0"

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

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

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

[[QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "12fbe86da16df6679be7521dfb39fbc861e1dc7b"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.1"

[[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 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

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

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

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

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

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

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

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

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

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

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

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

[[XMLDict]]
deps = ["EzXML", "IterTools", "OrderedCollections"]
git-tree-sha1 = "d9a3faf078210e477b291c79117676fca54da9dd"
uuid = "228000da-037f-5747-90a9-8195ccbf91a5"
version = "0.4.1"

[[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:
# ╟─059c84bc-5d87-4781-b0f2-15563bfae09d
# ╟─6e38cf1f-0c8e-47b9-baa1-e9525f0a1faa
# ╠═73112674-61f9-4915-9b09-216c70df0e74
# ╟─c9770fc4-e460-46d9-ba6b-2ab0fc17830e
# ╠═6bce5873-d064-46ea-9969-7809051300db
# ╟─5f0e15fc-1107-11ec-2a58-3351c91877fc
# ╠═f0b7e725-b433-458a-a5bf-efd5dbc990ce
# ╟─b5a00124-c994-4b79-99dc-6b966bec88cd
# ╠═ce999823-acab-47bb-bb15-babb8018a06a
# ╠═2d497d32-b088-4027-ad5b-0c4597db3922
# ╠═dbb95d40-6517-491b-b6a9-abccbd2a3a78
# ╟─e0acb2d4-668e-40aa-ac70-de4527c7a891
# ╟─423f35bd-6b79-464d-a8f6-a2405e5ba3ce
# ╠═c351704b-4bf8-4fc0-bda2-c2b1970c2a11
# ╟─4d88a451-f350-4a72-ae5b-8698b0527511
# ╠═45f4378e-9299-429d-b361-8f99b2be6efd
# ╟─6d5b5548-7cbc-435e-8058-ebd529f839f1
# ╠═372694fa-60a8-4305-a922-f020a4deb896
# ╟─a8cfd44f-5e40-4624-b98e-00de4a623af2
# ╠═7fe85397-51a5-4d82-86ee-f6881dcf8808
# ╠═47da227d-72ab-4eae-bb27-3a61d842bb03
# ╟─a5611f56-cca3-4a06-842e-715a89d401e4
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: JeannotJeannot/Learn

max_issues_repo_name: JeannotJeannot/Learn

max_forks_repo_name: JeannotJeannot/Learn


text:
module Polyopt

export solve_mosek, monomials
export MomentProb, momentprob, momentprob_chordal, momentprob_chordalembedding
export BSOSProb, bsosprob_chordal

import Base.^, Base.start, Base.next, Base.done, Base.sub, Base.length

include("polynomial.jl")
include("cliques.jl")

immutable MomentProb{T<:Number}
    order :: Int
    basis :: Vector{Poly{Int}}
    obj   :: AbstractArray{T}
    mom   :: Array{Any,1}
    eq    :: Array{Any,1}
end

immutable BSOSProb{T<:Number}
    degree :: Int
    order  :: Int
    obj    :: AbstractArray{T}
    Al     :: Array{Any,1}
    El     :: Array{Any,1}
    As     :: Array{Any,1}
    lb     :: Array{Any,1}
end

include("solver_mosek.jl")
include("latex.jl")

immutable MonomialPowers
    n :: Int
    d :: Int
end

monomialpowers(n, d) = MonomialPowers(n, d)
start(m::MonomialPowers) = zeros(Int, m.n)
done(m::MonomialPowers, state::Vector{Int}) = (state[1] == m.d+1)
length(m::MonomialPowers) = binomial(m.n + m.d, m.d)

function next(m::MonomialPowers, powers::Vector{Int})

    # find index of element to increment
    k = m.n
    while (k>1)
        if sum(view(powers,1:k)) < m.d
            break
        end
        k -= 1
    end

    state = copy(powers)
    state[k+1:end] = 0
    state[k] += 1

    powers, state
end

function monomials{T<:Number}(deg::Int, vars::Array{Poly{T},1})
    n = length(vars)
    I = [ findfirst(vi.alpha) for vi=vars ]

    v = Array(Poly{Int}, binomial(n+deg, deg))
    alpha = zeros(Int, 1, vars[1].n)
    for (i,a)=enumerate(monomialpowers(n, deg))
        alpha[I] = a
        v[i] = Poly(vars[1].syms, copy(alpha), [1])
        alpha[I] = 0
    end
    v
end

function moment(order::Int, syms::Symbols)
    v = monomials(order, variables(syms))
    v*v'
end

function moment(order::Int, syms::Symbols, I::Array{Int})
    v = monomials(order, variables(syms)[I])
    v*v'
end

function moment{T<:Number}(order::Int, syms::Symbols, p::Poly{T})
    v = monomials(order - (p.deg+1) >> 1, variables(syms))
    p*v*v'
end

function moment{T<:Number}(order::Int, syms::Symbols, p::Poly{T}, I::Array{Int})
    v = monomials(order - (p.deg+1) >> 1, variables(syms)[I])
    p*v*v'
end

function basis_index(a::Vector{Int}, degree::Int)
    n = length(a)
    idx = 1

    @inbounds begin
    for j=1:n-1
        if a[j] > 0
            #dj = degree - sum(a[1:j-1])
            dj = degree - sum(view(a,1:j-1))
            for i=0:a[j]-1
                nj = n-j
                di = dj - i
                idx += binomial(nj+di, di)
            end
        end
    end

    end
    idx + a[n]
end

basis_index(p::Poly, degree::Int) = basis_index(vec(p.alpha), degree)

linear_index(p::Poly, d::Int) = Int[ basis_index(vec(p.alpha[i,:]), d) for i=1:p.m ]

vectorize{T<:Number}(p::Polyopt.Poly{T}, d::Int) = sparsevec(linear_index(p, d), p.c, binomial(p.n+d,d))

function vectorize{T<:Number}(A::AbstractArray{Poly{T}}, d::Int)
    dims = length(size(A))
    if dims==1
        m, n = length(A), 1
    elseif dims==2
        m, n = size(A)
    else
        error("dimension must be <= 2")
    end
    subi = Int[]
    subj = Int[]
    val  = T[]
    for j=1:n
        for i=1:m
            I = linear_index(A[i,j], d)
            v = A[i,j].c
            push!(subi, (i+(j-1)*m)*ones(Int, length(v))...)
            push!(subj, I...)
            push!(val, v...)
        end
    end
    sparse(subi, subj, val, m*n, binomial(A[1].n+d,d))
end

function momentprob{S,T,U}(order::Int, obj::Poly{S}, pineq::Array{Poly{T},1}, peq::Array{Poly{U},1})

    v = monomials(2*order, variables(obj.syms))
    obj.deg <= 2*order || error("obj has degree higher than 2*order")

    p = vectorize(obj, 2*order)
    mom = Array(Any, length(pineq)+1)
    mom[1] = vectorize(moment(order, obj.syms), 2*order)
    for k=1:length(pineq)
       pineq[k].deg <= 2*order || error("pineq[$(k)] has degree higher than 2*order")
       mom[k+1] = vectorize(moment(order, pineq[k].syms, pineq[k]), 2*order)
    end

    momeq = Array(Any, length(peq))
    for k=1:length(peq)
        peq[k].deg <= 2*order || error("peq[$(k)] has degree higher than 2*order")
        momeq[k] = vectorize(moment(order, peq[k].syms, peq[k]), 2*order)
    end

    MomentProb(order, v, p, mom, momeq)
end

momentprob{S}(order::Int, obj::Poly{S}) =
    momentprob(order, obj, Poly{Int}[], Poly{Int}[])

momentprob{S,T}(order::Int, obj::Poly{S}, pineq::Array{Poly{T},1}) =
    momentprob(order, obj, pineq, Poly{Int}[])

function momentprob_chordal{S,T,U}(order::Int, cliques::Array{Array{Int,1},1}, obj::Poly{S},
                                   pineq::Array{Poly{T},1}, peq::Array{Poly{U},1})

    v = monomials(2*order, variables(obj.syms))
    obj.deg <= 2*order || error("obj has degree higher than 2*order")

    p = vectorize(obj, 2*order)
    mom = []
    for k=1:length(cliques)
        push!(mom,vectorize(moment(order, obj.syms, cliques[k]), 2*order))
    end

    for k=1:length(pineq)
        pineq[k].deg <= 2*order || error("pineq[$(k)] has degree higher than 2*order")
        for j=clique_index(cliques, find(sum(pineq[k].alpha,1)), obj.n)
            push!(mom,vectorize(moment(order, pineq[k].syms, pineq[k], cliques[j]), 2*order))
        end
    end

    momeq = []
    for k=1:length(peq)
        peq[k].deg <= 2*order || error("peq[$(k)] has degree higher than 2*order")
        for j=clique_index(cliques, find(sum(peq[k].alpha,1)), obj.n)
            push!(momeq, vectorize(moment(order, peq[k].syms, peq[k], cliques[j]), 2*order))
        end
    end

    MomentProb(order, v,  p, mom, momeq)
end

momentprob_chordal{S,T}(order::Int, cliques::Array{Array{Int,1},1}, obj::Poly{S}, pineq::Array{Poly{T},1}) =
    momentprob_chordal(order, cliques, obj, pineq, Poly{Int}[])

momentprob_chordal{S}(order::Int, cliques::Array{Array{Int,1},1}, obj::Poly{S}) =
    momentprob_chordal(order, cliques, obj, Poly{Int}[], Poly{Int}[])

function correlative_sparsity{S,T}(obj::Poly{S}, p::Array{Poly{T},1})
    A = eye(Int,obj.n,obj.n)

    for k=1:obj.m
        I = find(obj.alpha[k,:])
        A[I,I] = 1
    end

    for j=1:length(p)
        I = find(sum(p[j].alpha,1))
        A[I,I] = 1
    end
    sparse(A)
end

function momentprob_chordalembedding{S,T,U}(order::Int, obj::Poly{S}, pineq::Array{Poly{T},1}, peq::Array{Poly{U},1})
    C = correlative_sparsity(obj, [pineq; peq])
    cliques = chordal_embedding(C)
    momentprob_chordal(order, cliques, obj, pineq, peq)
end

momentprob_chordalembedding{S,T}(order::Int, obj::Poly{S}, pineq::Array{Poly{T},1}) =
    momentprob_chordalembedding(order, obj, pineq, Poly{Int}[])

momentprob_chordalembedding{S}(order::Int, obj::Poly{S}) =
    momentprob_chordalembedding(order, obj, Poly{Int}[], Poly{Int}[])

function symbol_restrict{T<:Number}(p::Poly{T}, syms::Symbols, I::Array{Int,1})

    if p.n == 0
        p
    else
        mask = ones(Int, p.n)
        mask[I] = 0
        J = find(p.alpha * mask .== 0)
        Poly{T}(syms, p.alpha[J,I], p.c[J])
    end
end

function bsosprob_chordal{S,T,V}(degree::Int, order::Int, cliques::Array{Array{Int,1}},
                                    obj::Poly{S}, pineq::Array{Poly{T},1}, peq::Array{Poly{V},1})

    n = obj.n

    allp = [pineq; peq]
    dmax = max(obj.deg, 2*order, degree*maximum([gi.deg for gi in allp]))
    m = binomial(n+dmax,dmax)

    x = variables(obj.syms)
    f = vectorize(obj, dmax)

    Ji = Vector{Int}[]
    for j=1:length(cliques)
        push!(Ji, [])
    end

    for (i, gi) in enumerate(pineq)
        for j=clique_index(cliques, find(sum(gi.alpha,1)), n)
            append!(Ji[j], i)
        end
    end

    Je = Vector{Int}[]
    for j=1:length(cliques)
        push!(Je, [])
    end

    for (i, hi) in enumerate(peq)
        for j=clique_index(cliques, find(sum(hi.alpha,1)), n)
            append!(Je[j], i)
        end
    end

    As, Al, El, lb = [], [], [], []
    for (j,c) in enumerate(cliques)
        symc = Symbols(obj.syms.names[c])
        xc = Poly{Int}[symbol_restrict(x[i], symc, c) for i in c ]

        u = monomials(order, xc)
        M = u*u'
        push!(As, vectorize(M, dmax))

        pj = [ symbol_restrict(pineq[i], symc, c) for i=Ji[j] ]

        p = vcat( pj,
                  [ (1-pji) for pji in pj ],
                  [ symbol_restrict(peq[i], symc, c) for i=Je[j] ])

        ab_d = monomialpowers(length(p), degree)
        #len_ab_d = length(ab_d)
        len_ab_d = length(ab_d) - (j > 1)

        #lbj = zeros(length(ab_d))
        lbj = zeros(len_ab_d)
        ai, aj, av = Int[], Int[], Float64[]
        lgta = binomial(obj.n+dmax,dmax)

        for (i, ab) in enumerate(ab_d)

            if j > 1
                if i==1
                    continue
                else
                    i = i - 1
                end
            end

            # for efficiency reasons we unroll   h_ab = prod(p .^ ab)
            h_ab = p[1]^ab[1]
            for l=2:length(ab)
                if ab[l] == 1
                    h_ab *= p[l]
                elseif ab[l] == 2
                    h_ab *= p[l]*p[l]
                elseif ab[l] > 2
                    h_ab *= p[l]^(ab[l])
                end
            end
            #println("h_$(ab)=$(h_ab)")
            ak = SparseVector(lgta, linear_index(h_ab, dmax), h_ab.c)

            append!(ai, i*ones(ak.nzind))
            append!(aj, ak.nzind)
            append!(av, ak.nzval)

            lbj[i] = (sum(ab[2*length(pj)+1:end]) > 0 ? -Inf : 0.0)
        end
        #push!(Al, sparse(ai, aj, av, length(ab_d), binomial(length(c)+dmax,dmax)))
        push!(Al, sparse(ai, aj, av, len_ab_d, binomial(length(c)+dmax,dmax)))
        push!(lb, lbj)

        k = Array(Int, binomial(length(c)+dmax, dmax))
        alpha = zeros(Int, n)
        for (i,a)=enumerate(monomialpowers(length(c), dmax))
            alpha[c] = a
            k[i] = basis_index(alpha, dmax)
            alpha[c] = 0
        end

        push!(El, SparseMatrixCSC{Int,Int}(m,length(k),collect(1:length(k)+1),k,ones(Int,length(k))))

    end

    BSOSProb(degree, order, f, Al, El, As, lb)
end

bsosprob_chordal{S,T}(degree::Int, order::Int, cliques::Array{Array{Int,1}}, obj::Poly{S}, pineq::Array{Poly{T},1}) =
    bsosprob_chordal(degree, order, cliques, obj, pineq, Poly{Int}[])

end

ext: jl


text:
export DefaultArray

# To solve this problem, you'll have to improve the `struct` definition and modify/add constructors
struct DefaultArray{N} <: AbstractArray{Any,N}
    parentarray::AbstractArray
    defaultvalue
end
DefaultArray(parentarray, defaultvalue) = DefaultArray{ndims(parentarray)}(parentarray, defaultvalue)

@inline Base.getindex(a::DefaultArray{N}, i::Vararg{Int,N}) where N = checkbounds(Bool, a, i...) ? a.parentarray[i...] : a.defaultvalue
Base.size(a::DefaultArray) = size(a.parentarray)


text:
function ggplot(dp, directory)
  DFres = dp.DF
  DF = dp.data.DF
  @rput DF
  @rput DFres
  @rput directory
  R"""
  library(ggplot2)
  library(cowplot)
  g <- ggplot(DFres, aes(x = x, y = mean)) +
  geom_histogram(data = DF, aes(x = VAF, y = ..density..), bins = 100,   fill = "gray68") +
  geom_line(col = "lightblue4", size = 1.0, alpha = 1.0) +
  geom_ribbon(aes(ymin = lq, ymax = uq), size = 0.2, col = "lightblue4", fill = "lightblue3", alpha = 0.7) +
  xlab("VAF") +
  ylab("Density")

  save_plot(filename = directory, plot = g)

  g

  """

end


text:
module GeneralizedGauss

using BasisFunctions, LinearAlgebra
import BasisFunctions: moment

export compute_moments,
    compute_gauss_rule, compute_gauss_rules,
    oscillatory_basis_exp, oscillatory_basis_cos, log_basis,
    jacobi_matrix, newton

import Base:
    eltype,
    length,
    size

include("rootfinding.jl")

include("sets.jl")

include("quadrule.jl")

include("representations.jl")

include("gengauss.jl")

end # module


text:
module SparseGaussianProcesses

using Requires

# include("transform.jl")
include("hyperprior.jl")
include("kernel.jl")
include("prior_basis.jl")
include("inducing.jl")
include("gp.jl")
include("loss.jl")

include("utils.jl")
include("gpu.jl")

function __init__()
  @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("plot.jl")
end

end # module


text:
using Documenter, Elfel

makedocs(
	modules = [Elfel],
	doctest = false, clean = true,
	format = Documenter.HTML(prettyurls = false),
	authors = "Petr Krysl",
	sitename = "Elfel.jl",
	pages = Any[
			"Home" => "index.md",
			"Tutorials" => "tutorials/tutorials.md",
			"How to guide" => "guide/guide.md",
			"Reference" => "man/reference.md",
			"Concepts" => "concepts/concepts.md"
		],
	)

deploydocs(
    repo = "github.com/PetrKryslUCSD/Elfel.jl.git",
)


text:
function comb_sort(seq::Array{Int64})::Array{Int64}
    gap = length(seq)
    swap = true
    while gap > 1 || swap
        gap = max(1, Int64(floor((gap / 1.25))))
        swap = false
        for i in (0:(length(seq)-gap)-1)
            if seq[i+1] > seq[(i+gap)+1]
                seq[i+1], seq[(i+gap)+1] = (seq[(i+gap)+1], seq[i+1])
                swap = true
            end
        end
    end
    return seq
end

function main()
    unsorted = [14, 11, 19, 5, 16, 10, 19, 12, 5, 12]
    expected = [5, 5, 10, 11, 12, 12, 14, 16, 19, 19]
    @assert(comb_sort(unsorted) == expected)
    println(join(["OK"], " "))
end

main()


text:
using availability

println(repeat("-", 50))
for i in (Int8, Int16, Int32, Int64,Int128)
    println("$i range = $(typemin(i)), $(typemax(i))")
end
println(repeat("-", 50))

for i in 1:1000
    a = Char(i)
    length(a) > 0 ? print("$(a) : ") : print("NOT_FOUND");break
end

a = """hello !!
I am Sumit Kawade.
I am learning "Julia language".
I find it easier for my next growth.
"""

b = a[1:10]
c = a[12:50]
d = replace(a, Pair(" ", "             "))

e = a[end - 50: end]

first(a, 10)

firstindex(a)
lastindex(a)

f = collect(a)

sizeof(f)

@. length(f)

g = "f: $f"

h = sort([1, 3 , 5, 2, 6, 2, 3 , 7])

i = sort(["abc", "bcd", "cde", "abd"])

j = findall("I", a)

k = join([1, 2, 3, 4, 5], " \n", "\n$(repeat("a", 10))")

codeunit(a, 1)

length(a)

nextind(a, 1, 4)

l = b"$a"

m = (4 ,2, 4)
m1 = m[2]
length(m)
sizeof(m)

sort(collect(m))
for i in m; println(i); end

n = (x=1, y=2, z=3, a=1, b=2)

typeof(n)

keys(n)

for a in keys(n)
    println("$a = $(n[a])")
end

function getBeforeAfter(in::NamedTuple, indx::Int)
    if length(in) < indx
        return nothing
    end
    in[indx-1], in[indx], length(in) < indx ? in[indx+1] : nothing

end

getBeforeAfter(n, 6)

o(x, y, z...) = (x, y, z)

o1 = o(1, 2, 3, 4, 5)

for i in o1
    println(i)
end

function withOptionalValues(x=1, y=2, z=3)
    println("x=$x, y=$y, z=$z")
end

withOptionalValues()

p = map(x -> x^2, 1:10)

p1 = map(1:100) do x
     x^2, (x+1)^2
 end

open("README.md", "r") do io
        println(readuntil(io, "232221"))
end

q = map(first ∘ uppercase, split("sumit is learning Julia language"))

square(x) = x.^2

1:10000 |> sum |> square

 r = [1, 2, 3, 4, 5]

r1 = square.(square.(r))
r2 = @. square(square(r))
isequal(r1, r2)
r1 == r2

try
    throw(AvailabilityException(1, "random"))
catch e
    println("Exception occurred: ", e)
finally
    println("exiting try/catch block")
end

# -------------------------------------------------
using JuliaDB

supplyTemp = Supply("item00001", "each", "node001", "onhand", 100)

t = table([1, 2, 3, 4, 5], [2, 3, 4, 5, 6]; names = [:x, :y])
# t = table((itemId=supplyTemp.itemId, uom=supplyTemp.uom, nodeId=supplyTemp.nodeId, type=supplyTemp.type, quantity=supplyTemp.quantity); pkey=[:itemId, :uom, :nodeId, :type])

save(t, "temp.dat")

filter(r -> r.x == 2, t)

push!(rows(t), (x=20, y=40))

# ------------------------------------------------
# Tasks - async / parallel/ multi-threaded programming

varsize = 1_000_000
function supplyProducer(c::Channel)
    for i in 1:varsize
        s = Supply(string("item_", rand(1:1000)), "each", string("node_", rand(1:100)), "onhand", rand(1:1000))
        put!(c, s)
    end
end

consumer = Channel(supplyProducer)

supplyArray = []

@time for i in 1:varsize
    u = take!(consumer)
    # println(u)
    push!(supplyArray, u)
end

u1 = table([], [], [], [], []; names = (:itemId, :uom, :nodeId, :type, :quantity))

@time for supply in supplyArray
    push!(rows(u1), (itemId = supply.itemId, uom = supply.uom, nodeId =
        supply.nodeId, type=supply.type, quantity=supply.quantity))
    # push!(([1], [2], [3], [4], [5]))
    # push!(rows(u1), (1, 2, 3, 4, 5))
end

save(u1, "temp.dat")
u1 = load("temp.dat")

u2 = filter(r -> r.quantity > 1 && r.itemId == "item_453", u1)

u3 = filter(r -> isequal(r.itemId, "item_453"), u1)
u31 = select(u3, :quantity) |> sum

u4 = select(filter(r -> isequal(r.itemId, "item_453"), u1), :quantity) |> sum

function myItemFilter(tab::IndexedTable)
    filter(r -> isequal(r.itemId, "item_453"), tab)
end

function mySelectQty(tab::IndexedTable)
    select(tab, :quantity)
end

u5 = mySelectQty(myItemFilter(u3)) |> sum

u5 = (mySelectQty ∘ myItemFilter)(u1) |> sum

u6 = ((t -> select(t, :quantity)) ∘ (t -> filter(r -> isequal(r.itemId, "item_453"), t)))(u1) |> sum

function squareAndPrint(arg)
    out = square(arg)
    println("square($arg) = $out")
end

u4 = Task(() -> squareAndPrint(6))
istaskstarted(u4)
istaskdone(u4)
istaskfailed(u4)
a = schedule(u4, 3)

# b = yieldto(u4) TODO it is not returning control back to the main thread. Also, function is not executing. process gets stuck
# yieldto(current_task())

println(b)

# ----------------------------------------
# types, parameteric types, abstract types, etc

u7 = ((applySupplyAdjusters) ∘ (t -> select(t, :quantity)) ∘ (t -> filter(r -> isequal(r.itemId, "item_453"), t)))(u1) |> sum
@time u8 = ((applySupplyAdjusters) ∘ (t -> select(t, :quantity)) ∘ (t -> filter(r -> isequal(r.itemId, "item_453"), t)))(u1) |> sum

# ----------------------------------------
# interfaces

struct Square
    num::Int
end

v1 = Square(10)

Base.iterate(in::Square, i=1) = i > in.num ? nothing : (i*i, i+1)

@benchmark for val in Square(10)
    # println("val=$(val)")
    val * val
end

125 in v1

using Statistics
mean(v1)
sum(v1)

mean(Square(2))
sum(Square(2))

Base.eltype(::Type{Square}) = Int
Base.length(s::Square) = s.num

collect(v1)

function Base.getindex(S::Square, i::Int)
    1 <=i <=S.num || throw(BoundsError(S, i))

    return i * i
end

Square(1000)[15]

Base.firstindex(S::Square) = 1
Base.lastindex(S::Square) = length(S::Square)

firstindex(Square(10000))
lastindex(Square(10010))

Square(10010)[end]

Base.getindex(S::Square, i) = [S[i1] for i1 in i]

Square(100)[(1, 5, 4, 3, 9)]

Square(100)[[1, 5, 4, 3, 9]]

max_stars_repo_name: kawadesumit/availability

max_issues_repo_name: kawadesumit/availability

max_forks_repo_name: kawadesumit/availability


text:
using JLD2
using Infiltrator
using LightGraphs
using Random
using JSON
using BenchmarkTools
using MultiAgentPathFinding

using AerialGroundPathFinding

const BASE_SEED = 6789
const MANHATTAN_NODEFILE = "../data/manhattan_node_attribs.json"
const MANHATTAN_WEIGHTFILE = "../data/manhattan_sparse_wts.jld2"

# Only caring about distance for now
const ALPHA_WEIGHT_DISTANCE = 1.0

const NCARS = parse(Int64, ARGS[1])
const NDRONES = parse(Int64, ARGS[2])
const OUTFILEDIR = ARGS[3]
const CAPACITY = parse(Int64, ARGS[4])
const TRIALS = parse(Int64, ARGS[5])

## Stage 1 paths are computed directly
## Stage 2 paths are drones over transit
function main()

    @load MANHATTAN_WEIGHTFILE manhattan_sparse_wts

    graph = LightGraphs.SimpleDiGraph(manhattan_sparse_wts)
    location_list = get_location2D_list(MANHATTAN_NODEFILE)

    for t = 1:TRIALS+1

        rng = MersenneTwister(BASE_SEED+t)

        ground_task_list = get_tasks_with_valid_path(graph, manhattan_sparse_wts, location_list, NCARS, rng)
        aerial_task_list = get_tasks_with_valid_path(graph, manhattan_sparse_wts, location_list, NDRONES, rng)

        env = CoordinatedMAPFEnv(ground_task_list=ground_task_list, aerial_task_list=aerial_task_list, road_graph=graph, road_graph_wts=manhattan_sparse_wts, location_list=location_list, alpha_weight_distance=ALPHA_WEIGHT_DISTANCE, car_capacity=CAPACITY)

        time_comps = Float64[]

        t1 = @elapsed dir_drone_paths = compute_independent_paths(env, env.aerial_task_list)
        t2 = @elapsed dir_car_paths = compute_independent_paths(env, env.ground_task_list)
        push!(time_comps, t1+t2)

        # Just set env ground and aerial paths to direct and update
        env.ground_paths = dir_car_paths
        env.aerial_paths = dir_drone_paths
        t3 = @elapsed set_ground_transit_graph!(env, [ddp.total_dist for ddp in dir_drone_paths])
        push!(time_comps, t3)

        initial_amapf_states = [AerialMAPFState(idx=1, ground_transit_idx=sg[1]) for (_, sg) in enumerate(env.gtg_drone_start_goal_idxs)]
        t4 = @elapsed aerial_pp_ordering = get_increasing_dist_ordering(dir_drone_paths)
        push!(time_comps, t4)

        t5 = @elapsed aerial_pp_soln = plan_prioritized_paths!(env, aerial_pp_ordering, initial_amapf_states, GroundTransitConstraint, PlanResult{AerialMAPFState,GroundTransitAction,Float64})
        push!(time_comps, t5)

        t6 = @elapsed update_aerial_ground_paths_cbs_solution!(env, aerial_pp_soln, dir_drone_paths)
        push!(time_comps, t6)

        total_dist = compute_total_cost(env, env.aerial_paths, env.ground_paths)

        res_dict = Dict{String,Float64}("compute_time" => sum(time_comps), "total_path_cost" => total_dist)

        if t > 1
            fname = string(OUTFILEDIR, "_Cars_", NCARS, "_Drones_", NDRONES, "_Capacity_", CAPACITY, "_Trial_", (t-1), ".json")
            println(fname)
            open(fname, "w") do f
                JSON.print(f, res_dict, 2)
            end
        end
    end
end

main()


text:
abstract type AbstractParticle{T,Dim} end
abstract type Shape{T<:AbstractFloat,Dim} end

struct Particle{T<:AbstractFloat,Dim,P<:PhysicalProperties,S<:Shape} <: AbstractParticle{T,Dim}
    medium::P
    shape::S
    # Enforce that the Dims and Types are all the same
    function Particle{T,Dim,P,S}(medium::P,shape::S) where {T,Dim,FieldDim,P<:PhysicalProperties{T,Dim,FieldDim},S<:Shape{T,Dim}}
        new{T,Dim,P,S}(medium,shape)
    end
end

struct Specie{T<:Real,Dim,P<:PhysicalProperties,S<:Shape} <: AbstractParticle{T,Dim}
  medium::P
  shape::S
  num_density::T # number density
end

shape(p::Particle) = p.shape
shape(p::CapsuleParticle) = p.outer.shape


text:
using JuMP
using CPLEX
#using Gurobi
using LightGraphs
using DelimitedFiles
using LinearAlgebra
using Random
using TimerOutputs

include("../src/algorithm.jl")
include("../src/RCG.jl")

const ϵ = 10^-5
const TIME_LIM = 7200
#const Optimizer = "CPLEX"
const Optimizer = "CPLEX"
#const GUROBI_ENV = Gurobi.Env()
const to = TimerOutput()
"For non-linear objective function, different choices arise for the objective of the local search step.
Set the following to LINEAR in case of KP or SP"
#const local_search_mode = "PRODUCT"
const local_search_mode = "LINEAR"

#-----------------------------------------------------------------------------------

function create_model(TIME_LIM)
  if Optimizer=="CPLEX"
    M = Model(CPLEX.Optimizer)
    MOI.set(M, MOI.RawParameter("CPX_PARAM_THREADS"), 1)
    MOI.set(M, MOI.RawParameter("CPX_PARAM_SCRIND"), 0)
    MOI.set(M, MOI.RawParameter("CPX_PARAM_MIPDISPLAY"), 1)
    MOI.set(M, MOI.RawParameter("CPX_PARAM_TILIM"), max(0,TIME_LIM))
  elseif Optimizer=="Gurobi"
    M = Model(optimizer_with_attributes(() -> Gurobi.Optimizer(GUROBI_ENV)))
    set_optimizer_attribute(M, "OutputFlag", 0)
    set_optimizer_attribute(M, "TimeLimit", max(0,TIME_LIM))
    set_optimizer_attribute(M, "Threads", 1)
    #M = Model(with_optimizer(Gurobi.Optimizer,GUROBI_ENV,OutputFlag= 0,TimeLimit=max(0,TIME_LIM), Threads=1))
  end
  return M
end

max_stars_repo_name: mjposs/min-max-min

max_issues_repo_name: mjposs/min-max-min

max_forks_repo_name: mjposs/min-max-min


text:
function value(obj::MultivariateDifferentiableFunction, arg0::Vector{DerivativeStructure})
    return jcall(obj, "value", DerivativeStructure, (Vector{DerivativeStructure},), arg0)
end

function value(obj::MultivariateFunction, arg0::Vector{jdouble})
    return jcall(obj, "value", jdouble, (Vector{jdouble},), arg0)
end


text:
module Dispersal
# Use the README as the module docs
@doc read(joinpath(dirname(@__DIR__), "README.md"), String) Dispersal

using ConstructionBase,
      Colors,
      Dates,
      DimensionalData,
      Distributed,
      Distributions,
      DocStringExtensions,
      FieldDefaults,
      FieldDocTables,
      FieldMetadata,
      Flatten,
      LinearAlgebra,
      LossFunctions,
      Mixers,
      PoissonRandom,
      Reexport,
      Setfield,
      Statistics

@reexport using DynamicGrids

using LossFunctions: ZeroOneLoss
using DimensionalData: Time

import Base: getindex, setindex!, lastindex, size, length, push!

import DynamicGrids: applyrule, applyrule!, applyrule, applyrule!,
       neighbors, sumneighbors, neighborhood, setneighbor!, mapsetneighbor!,
       radius, gridsize, mask, overflow, cellsize, ruleset, inbounds, extent,
       currenttime, currenttimestep, timestep, tspan,
       buffer, WritableGridData, aux, unwrap

import ConstructionBase: constructorof

import FieldMetadata: @default, @description, @bounds, @flattenable,
                      default, description, bounds, flattenable

export AbstractDispersalKernel, DispersalKernel

export KernelFormulation, ExponentialKernel, GeometricKernel, GaussianKernel

export DistanceMethod, CentroidToCentroid, CentroidToArea, AreaToArea, AreaToCentroid

export AbstractInwardsDispersal, InwardsBinaryDispersal, InwardsPopulationDispersal,
       SwitchedInwardsPopulationDispersal

export AbstractOutwardsDispersal, OutwardsBinaryDispersal, OutwardsPopulationDispersal

export AlleeExtinction, JumpDispersal, HumanDispersal

export BatchGroups, HeirarchicalGroups

export ExactExponentialGrowth, ExactLogisticGrowth

export DiscreteGrowth, DiscreteGrowthMap, DiscreteGrowth2, DiscreteGrowth3

export MaskGrowthMap, ExactExponentialGrowthMap, ExactLogisticGrowthMap, ExactLogisticGrowthMap2, ExactLogisticGrowthMap3

export SurvLogLogistic, SurvLogLogisticMap, SurvLogLogisticMap2

export GrowthSurvLogLogisticMap3

export MaskSurvMap

export SelectionGradientSurv, SelectionGradientSurvMap

export SelectionGradient1locusSurv, SelectionGradient1locusSurvMap, LandeVariable

export DeltaAlleleFrequencySurv, DeltaAlleleFrequencySurvMap, DeltaAlleleFrequencySurvMap_noFastmath

export MatingPopulation, MatingDispersal

export LayerCopy

export Parametriser, AbstractObjective, SimpleObjective, RegionObjective, RegionOutput,
       ColorRegionFit, Accuracy

export ThreadedReplicates, DistributedReplicates, SingleCoreReplicates

export targets, predictions

const FIELDDOCTABLE = FieldDocTable((:Description, :Default, :Bounds),
                                    (description, default, bounds);
                                    truncation=(100,40,100))

# Documentation templates
@template TYPES =
    """
    $(TYPEDEF)
    $(DOCSTRING)
    """

include("rules/mixins.jl")
include("downsampling.jl")
include("layers.jl")
include("rules/kernel/common.jl")
include("rules/kernel/inwards.jl")
include("rules/kernel/outwards.jl")
include("rules/growth.jl")
include("rules/human.jl")
include("rules/jump.jl")
include("rules/allee.jl")
include("rules/survival.jl")
include("rules/alleleFrequency.jl")
include("rules/selectionGradient.jl")
include("rules/matingPopulation.jl")
include("optimisation/optimisation.jl")
include("optimisation/objectives.jl")
include("optimisation/output.jl")

end # module


text:
# TextREST.jl
#
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

module TextREST

using Text
using Stage, Ollam, DataStructures
using HttpServer, JSON
using JLD, HDF5

export text_rest_server

# function to get resource name from HTTP Request
function get_resource_name(req::Request)
    return strip(split(split(req.resource, "?")[1], "/")[2])
end

# function to return TEXT REST Server object
function text_rest_server(fextractor::Function, model::LinearModel{String})
    http = HttpHandler() do req::Request, res::Response
        result      = Dict()
        resource    = get_resource_name(req)
        result["resource"]  = resource
        if ismatch(r"^$", resource)
            result["resource"]  = "index"
        elseif ismatch(r"^lid$", resource)
            lid(req, result, fextractor, model)
        else
            result["resource"]  = 404
        end
        return Response(JSON.json(result))
    end
    http.events["error"]    = (client, err) -> println(err)
    http.events["listen"]   = (saddr) -> println("Running on http://$saddr (Press CTRL+C to quit)")

    server  = Server(http)

    return server
end

include("lid.jl")

end

max_forks_repo_name: trevorlewis/TEXT-Language-REST


text:
# Regularize

"""
  Regularize(f::ProximableFunction, rho::Real, a::AbstractArray)

Given function `f`, returns `g(x) = f(x) + (rho/2)||x-a||²`.
"""

immutable Regularize{T <: ProximableFunction, S <: Real, A <: Union{Real, AbstractArray}} <: ProximableFunction
  f::T
  rho::S
  a::A
  function Regularize(f::T, rho::S, a::A)
    if rho <= 0.0
      error("parameter ρ must be positive")
    else
      new(f, rho, a)
    end
  end
end

Regularize{T <: ProximableFunction, S <: Real, A <: AbstractArray}(f::T, rho::S, a::A) = Regularize{T, S, A}(f, rho, a)

Regularize{T <: ProximableFunction, S <: Real}(f::T, rho::S=1.0, a::S=0.0) = Regularize{T, S, S}(f, rho, a)

function (g::Regularize){T <: RealOrComplex}(x::AbstractArray{T})
	return g.f(x) + g.rho/2*vecnorm(x-g.a)^2
end

function prox!{T <: RealOrComplex}(g::Regularize, x::AbstractArray{T}, y::AbstractArray{T}, gamma::Real=1.0)
  gr= g.rho*gamma
  gr2= 1/(1+gr)
  v = prox!(g.f, gr2.*(x+gr.*g.a), y, gr2*gamma)
  return v + g.rho/2*vecnorm(y-g.a)^2
end

function prox_naive{T <: RealOrComplex}(g::Regularize, x::AbstractArray{T}, gamma::Real=1.0)
  y, v = prox_naive(g.f, x./(1+gamma*g.rho)+g.a./(1/(gamma*g.rho)+1), gamma/(1+gamma*g.rho))
  return y, v + g.rho/2*vecnorm(y-g.a)^2
end

fun_name(f::Regularize) = string("Regularized ", fun_name(f.f))
fun_dom(f::Regularize) = fun_dom(f.f)
fun_expr(f::Regularize) = string(fun_expr(f.f),"+(ρ/2)||x-a||²")
fun_params(f::Regularize) = "ρ = $(f.rho), λ = $(f.f.lambda), a = $( typeof(f.a)<:Real ? f.a :typeof(f.a) )"


text:
struct SimulationResultsReference
    ref::Dict
    results_folder::String
    chronologies::Dict
    base_powers::Dict
end

function SimulationResultsReference(sim::Simulation; kwargs...)
    date_run = convert(String, last(split(dirname(sim.internal.raw_dir), "/")))
    ref = make_references(sim, date_run; kwargs...)
    chronologies = Dict()
    base_powers = Dict()
    for (stage_number, stage_name) in sim.sequence.order
        stage = get_stage(sim, stage_name)
        interval = get_stage_interval(sim, stage_name)
        resolution = PSY.get_forecasts_resolution(get_sys(stage))
        chronologies["stage-$stage_name"] = convert(Int, (interval / resolution))
        base_powers[stage_name] = PSY.get_basepower(sim.stages[stage_name].sys)
    end
    return SimulationResultsReference(
        ref,
        sim.internal.results_dir,
        chronologies,
        base_powers,
    )
end

# internal function for differentiating variables from duals in file names
function _concat_dual(duals::Vector{Symbol})
    dual = []
    for d in duals
        dual = vcat(dual, "dual_" * String(d))
    end
    return dual
end

# internal function for differentiating variables from parameters in file names
function _concat_param(params::Vector{Symbol})
    param = []
    for p in params
        param = vcat(param, "parameter_" * String(p))
    end
    return param
end

"""
    make_references(sim::Simulation, date_run::String; kwargs...)

Creates a dictionary of variables with a dictionary of stages
that contains dataframes of date/step/and desired file path
so that the results can be parsed sequentially by variable
and stage type.

**Note:** make_references can only be run after run_sim_model
or else, the folder structure will not yet be populated with results

# Arguments
- `sim::Simulation = sim`: simulation object created by Simulation()
- `date_run::String = "2019-10-03T09-18-00"``: the name of the file created
that contains the specific simulation run of the date run and "-test"

# Example
```julia
sim = Simulation("Test", 7, stages, "/Users/yourusername/Desktop/"; system_to_file = false)
execute!(sim::Simulation; kwargs...)
references = make_references(sim, "2019-10-03T09-18-00-test")
```
"""
function make_references(sim::Simulation, date_run::String)
    sim.internal.date_ref[1] = sim.initial_time
    sim.internal.date_ref[2] = sim.initial_time
    references = Dict()
    for (stage_number, stage_name) in sim.sequence.order
        stage = sim.stages[stage_name]
        references["stage-$stage_name"] = make_result_reference(stage, sim)
    end
    return references
end

function make_result_reference(
    stage::Stage{T},
    sim::Simulation,
) where {T <: PowerSimulationsOperationsProblem}
    stage_number = get_number(stage)
    stage_name = get_stage_name(sim, stage)
    stage_container = get_psi_container(stage)
    variables = Dict{Symbol, Any}()
    interval = get_stage_interval(sim, stage_name)
    variable_names = (collect(keys(stage_container.variables)))
    if !is_milp(get_psi_container(stage))
        constraint_duals = get_constraint_duals(stage_container.settings)
        constraint_duals_names = Symbol.(_concat_dual(constraint_duals))
        variable_names = vcat(variable_names, constraint_duals_names)
    end
    params = collect(keys(get_parameters_value(get_psi_container(stage))))
    param_keys = Symbol.(_concat_param(params))
    variable_names = vcat(variable_names, param_keys)
    for name in variable_names
        variables[name] = DataFrames.DataFrame(
            Date = Dates.DateTime[],
            Step = String[],
            File_Path = String[],
        )
    end
    for s in 1:(sim.steps)
        for run in 1:(stage.internal.executions)
            sim.internal.current_time = sim.internal.date_ref[stage_number]
            for name in variable_names
                full_path = joinpath(
                    sim.internal.raw_dir,
                    "step-$(s)-stage-$(stage_name)",
                    replace_chars("$(sim.internal.current_time)", ":", "-"),
                    "$(name).feather",
                )
                if isfile(full_path)
                    date_df = DataFrames.DataFrame(
                        Date = sim.internal.current_time,
                        Step = "step-$(s)",
                        File_Path = full_path,
                    )
                    variables[name] = vcat(variables[name], date_df)
                else
                    @error "$full_path, does not contain any simulation result raw data"
                end
            end
            sim.internal.run_count[s][stage_number] += 1
            sim.internal.date_ref[stage_number] =
                sim.internal.date_ref[stage_number] + interval
        end
    end
    return variables
end

make_result_reference(stage::Stage{T}, sim::Simulation) where {T} = nothing

struct SimulationResults <: IS.Results
    base_power::Float64
    variable_values::Dict{Symbol, DataFrames.DataFrame}
    total_cost::Dict
    optimizer_log::Dict
    time_stamp::DataFrames.DataFrame
    dual_values::Dict{Symbol, Any}
    results_folder::Union{Nothing, String}
    parameter_values::Dict{Symbol, DataFrames.DataFrame}
end

IS.get_base_power(result::SimulationResults) = result.base_power
IS.get_variables(result::SimulationResults) = result.variable_values
IS.get_total_cost(result::SimulationResults) = result.total_cost
IS.get_optimizer_log(results::SimulationResults) = results.optimizer_log
IS.get_time_stamp(result::SimulationResults) = result.time_stamp
get_duals(result::SimulationResults) = result.dual_values
IS.get_parameters(result::SimulationResults) = result.parameter_values

function deserialize_sim_output(file_path::String)
    path = joinpath(file_path, "output_references")
    list = setdiff(
        collect(readdir(path)),
        ["results_folder.json", "chronologies.json", "base_power.json"],
    )
    ref = Dict()
    for stage in list
        ref[stage] = Dict{Symbol, Any}()
        for variable in collect(readdir(joinpath(path, stage)))
            var = splitext(variable)[1]
            ref[stage][Symbol(var)] = Feather.read(joinpath(path, stage, variable))
            ref[stage][Symbol(var)][!, :Date] =
                convert(Array{Dates.DateTime}, ref[stage][Symbol(var)][!, :Date])
        end
    end
    results_folder = read_json(joinpath(path, "results_folder.json"))
    chronologies = Dict{Any, Any}(read_json(joinpath(path, "chronologies.json")))
    base_powers = Dict{Any, Any}(read_json(joinpath(path, "base_power.json")))
    sim_output = SimulationResultsReference(ref, results_folder, chronologies, base_powers)
    return sim_output
end

# internal function to parse through the reference dictionary and grab the file paths
function _read_references(
    results::Dict,
    list::Array,
    stage::String,
    step::Array,
    references::Dict,
    time_length::Int,
)

    for name in list
        date_df = references[stage][name]
        step_df = DataFrames.DataFrame(
            Date = Dates.DateTime[],
            Step = String[],
            File_Path = String[],
        )
        for n in 1:length(step)
            step_df = vcat(step_df, date_df[date_df.Step .== step[n], :])
        end
        results[name] = DataFrames.DataFrame()
        for (ix, time) in enumerate(step_df.Date)
            file_path = step_df[ix, :File_Path]
            var = Feather.read("$file_path")
            results[name] = vcat(results[name], var[1:time_length, :])
        end
    end
    return results
end
# internal function to parse through the reference dictionary and grab the file paths
function _read_references(
    results::Dict,
    list::Array,
    stage::String,
    references::Dict,
    time_length::Int,
)
    for name in list
        date_df = references[stage][name]
        results[name] = DataFrames.DataFrame()
        for (ix, time) in enumerate(date_df.Date)
            file_path = date_df[ix, :File_Path]
            var = Feather.read(file_path)
            var_length = min(time_length, size(var, 1))
            results[name] = vcat(results[name], var[1:var_length, :])
        end
    end
    return results
end
# internal function to remove the overlapping results and only use the most recent
function _read_time(file_path::String, time_length::Number)
    time_file_path = joinpath(dirname(file_path), "time_stamp.feather")
    temp_time_stamp = Feather.read("$time_file_path")
    time_stamp = temp_time_stamp[(1:time_length), :]
    time_stamp = convert.(Dates.DateTime, time_stamp)
    return time_stamp
end

"""
    load_simulation_results(stage, step, variable, SimulationResultsReference)

This function goes through the reference table of file paths and
aggregates the results over time into a struct of type OperationsProblemResults
for the desired step range and variables

# Arguments
- `SimulationResultsReference::SimulationResultsReference`: the container for the reference dictionary created in execute!
- `stage_number::Int = 1``: The stage of the results getting parsed: 1 or 2
- `step::Array{String} = ["step-1", "step-2", "step-3"]`: the steps of the results getting parsed
- `variable::Array{Symbol} = [:P__ThermalStandard, :P__RenewableDispatch]`: the variables to be parsed

# Example
```julia
stage = "stage-1"
step = ["step-1", "step-2", "step-3"] # has to match the date range
variable = [:P__ThermalStandard, :P__RenewableDispatch]
results = load_simulation_results(stage,step, variable, SimulationResultsReference)
```
# Accepted Key Words
"""

function load_simulation_results(
    path::String,
    stage_name::String,
    step::Array,
    variable::Array;
    kwargs...,
)
    sim_results = deserialize_sim_output(path)
    load_simulation_results(sim_results, stage_name, step, variable; kwargs...)
end

function load_simulation_results(
    sim_output::SimulationResultsReference,
    stage_name::String,
    step::Array,
    variable::Array;
    kwargs...,
)
    results_folder = sim_output.results_folder
    stage = "stage-$stage_name"
    references = sim_output.ref
    base_power = sim_output.base_powers[stage_name]
    variables = Dict{Symbol, DataFrames.DataFrame}()
    duals = Dict{Symbol, Any}()
    params = Dict{Symbol, DataFrames.DataFrame}()
    time_stamp = DataFrames.DataFrame(Range = Dates.DateTime[])
    time_length = sim_output.chronologies[stage]
    dual = _find_duals(collect(keys(references[stage])))
    param = _find_params(collect(keys(references[stage])))
    variable = setdiff(variable, vcat(param, dual))
    for l in 1:length(variable)
        date_df = references[stage][variable[l]]
        step_df = DataFrames.DataFrame(
            Date = Dates.DateTime[],
            Step = String[],
            File_Path = String[],
        )
        for n in 1:length(step)
            step_df = vcat(step_df, date_df[date_df.Step .== step[n], :])
        end
        variables[(variable[l])] = DataFrames.DataFrame()
        for (ix, time) in enumerate(step_df.Date)
            file_path = step_df[ix, :File_Path]
            var = Feather.read("$file_path")
            variables[(variable[l])] = vcat(variables[(variable[l])], var[1:time_length, :])
            if l == 1
                time_stamp = vcat(time_stamp, _read_time(file_path, time_length))
                check_file_integrity(dirname(file_path))
            end
        end
    end
    time_stamp[!, :Range] = convert(Array{Dates.DateTime}, time_stamp[!, :Range])
    file_path = dirname(references[stage][variable[1]][1, :File_Path])
    optimizer = read_json(joinpath(file_path, "optimizer_log.json"))
    obj_value = Dict{Symbol, Any}(:OBJECTIVE_FUNCTION => optimizer["obj_value"])
    duals = _read_references(duals, dual, stage, step, references, time_length)
    param_values = _read_references(params, param, stage, step, references, time_length)
    return SimulationResults(
        base_power,
        variables,
        obj_value,
        optimizer,
        time_stamp,
        duals,
        results_folder,
        param_values,
    )
end

"""
    load_simulation_results(file_path, stage)

This function goes through the reference table of file paths and
aggregates the results over time into a struct of type OperationsProblemResults

# Arguments
- `file_path::String`: the file path to the dated folder with the raw results
- `stage_number::String`: The stage of the results getting parsed

# Example
```julia
execute!(simulation)
results = load_simulation_results("file_path", "stage_name")
```
# Accepted Key Words
"""
function load_simulation_results(path::String, stage_name::String, kwargs...)
    sim_results = deserialize_sim_output(path)
    load_simulation_results(sim_results, stage_name; kwargs...)
end
"""
    load_simulation_results(sim_output, stage)

This function goes through the reference table of file paths and
aggregates the results over time into a struct of type OperationsProblemResults

# Arguments
- `sim_output::SimulationResultsReference`: the container for the reference dictionary created in execute!
- `stage_number::String`: The stage of the results getting parsed

# Example
```julia
sim_output = execute!(simulation)
results = load_simulation_results(sim_output, "stage_name")
```
# Accepted Key Words
"""
function load_simulation_results(
    sim_output::SimulationResultsReference,
    stage_name::String;
    kwargs...,
)
    results_folder = sim_output.results_folder
    stage = "stage-$stage_name"
    references = sim_output.ref
    base_power = sim_output.base_powers[stage_name]
    variables = Dict{Symbol, DataFrames.DataFrame}()
    duals = Dict{Symbol, Any}()
    params = Dict{Symbol, DataFrames.DataFrame}()
    variable = (collect(keys(references[stage])))
    dual = _find_duals(variable)
    param = _find_params(collect(keys(references[stage])))
    variable = setdiff(variable, vcat(param, dual))
    time_stamp = DataFrames.DataFrame(Range = Dates.DateTime[])
    time_length = sim_output.chronologies[stage]

    for l in 1:length(variable)
        date_df = references[stage][variable[l]]
        variables[(variable[l])] = DataFrames.DataFrame()
        for (ix, time) in enumerate(date_df.Date)
            file_path = date_df[ix, :File_Path]
            var = Feather.read(file_path)
            variables[(variable[l])] = vcat(variables[(variable[l])], var[1:time_length, :])
            if l == 1
                time_stamp = vcat(time_stamp, _read_time(file_path, time_length))
                check_file_integrity(dirname(file_path))
            end
        end
    end
    time_stamp[!, :Range] = convert(Array{Dates.DateTime}, time_stamp[!, :Range])
    file_path = dirname(references[stage][variable[1]][1, :File_Path])
    optimizer = read_json(joinpath(file_path, "optimizer_log.json"))
    obj_value = Dict{Symbol, Any}(:OBJECTIVE_FUNCTION => optimizer["obj_value"])
    param_values = _read_references(params, param, stage, references, time_length)
    duals = _read_references(duals, dual, stage, references, time_length)
    return SimulationResults(
        base_power,
        variables,
        obj_value,
        optimizer,
        time_stamp,
        duals,
        results_folder,
        param_values,
    )
end
"""
    check_file_integrity(path::String)

Checks the hash value for each file made with the file is written with the new hash_value to verify the file hasn't been tampered with since written

# Arguments
- `path::String`: this is the folder path that contains the results and the check.sha256 file
"""
function check_file_integrity(path::String)
    matched = true
    for file_info in read_file_hashes(path)
        filename = file_info["filename"]
        expected_hash = file_info["hash"]
        actual_hash = compute_sha256(filename)
        if expected_hash != actual_hash
            @error "hash mismatch for file" filename expected_hash actual_hash
            matched = false
        end
    end

    if !matched
        throw(IS.HashMismatchError("The hash value in the written files does not match the read files, results may have been tampered."))
    end
end

function get_variable_names(sim::Simulation, stage::Any)
    return get_variable_names(sim.stages[stage].internal.psi_container)
end

function get_reference(
    sim_results::SimulationResultsReference,
    stage::String,
    step::Int,
    variable::Symbol,
)
    file_paths = sim_results.ref["stage-$stage"][variable]
    return filter(file_paths -> file_paths.Step == "step-$step", file_paths)[:, :File_Path]
end

get_psi_container(sim::Simulation, stage::Any) = sim.stages[stage].internal.psi_container

"""
    write_results(results::SimulationResults)

Exports Simulations Results to the path where they come from in the results folder

# Arguments
- `results::SimulationResults`: results from the simulation
- `save_path::String`: folder path where the files will be written
- `results_folder`: name of the folder where the results will be written

# Accepted Key Words
- `file_type = CSV`: only CSV and featherfile are accepted
"""

function IS.write_results(res::SimulationResults; kwargs...)
    folder_path = res.results_folder
    if !isdir(folder_path)
        throw(IS.ConflictingInputsError("Specified path is not valid. Set up results folder."))
    end
    write_data(res.variable_values, res.time_stamp, folder_path; kwargs...)
    write_optimizer_log(res.optimizer_log, folder_path)
    write_data(res.time_stamp, folder_path, "time_stamp"; kwargs...)
    write_data(res.base_power, folder_path)
    write_data(res.dual_values, folder_path; kwargs...)
    write_data(res.parameter_values, folder_path; kwargs...)
    files = collect(readdir(folder_path))
    compute_file_hash(folder_path, files)
    @info("Files written to $folder_path folder.")
    return
end

function serialize_sim_output(sim_results::SimulationResultsReference)
    file_path = mkdir(joinpath(dirname(sim_results.results_folder), "output_references"))
    for (k, stage) in sim_results.ref
        try
            for (i, v) in stage
                path = joinpath(file_path, "$k")
                !isdir(path) && mkdir(path)
                # TODO: Remove this line. There shouldn't be empties coming here.
                !isempty(v) && Feather.write(joinpath(path, "$i.feather"), v)
            end
        catch
            @warn("Results Reference not compatible with serialization")
        end
    end
    JSON.write(
        joinpath(file_path, "results_folder.json"),
        JSON.json(sim_results.results_folder),
    )
    JSON.write(
        joinpath(file_path, "chronologies.json"),
        JSON.json(sim_results.chronologies),
    )
    JSON.write(joinpath(file_path, "base_power.json"), JSON.json(sim_results.base_powers))
end

# writes the results to CSV files in a folder path, but they can't be read back
function write_to_CSV(res::SimulationResults; kwargs...)
    folder_path = res.results_folder
    if !isdir(folder_path)
        throw(IS.ConflictingInputsError("Specified path is not valid. Set up results folder."))
    end
    for (k, v) in IS.get_variables(res)
        if decode_symbol(k)[1] == "P"
            IS.get_variables(res)[k] = IS.get_base_power(res) .* v
        end
    end
    for (p, v) in IS.get_parameters(res)
        IS.get_parameters(res)[p] = IS.get_base_power(res) .* v
    end
    write_data(
        IS.get_variables(res),
        res.time_stamp,
        folder_path;
        file_type = CSV,
        kwargs...,
    )
    write_optimizer_log(IS.get_total_cost(res), folder_path)
    write_data(
        IS.get_time_stamp(res),
        folder_path,
        "time_stamp";
        file_type = CSV,
        kwargs...,
    )
    write_data(get_duals(res), folder_path; file_type = CSV, kwargs...)
    write_data(IS.get_parameters(res), folder_path; file_type = CSV, kwargs...)
    files = collect(readdir(folder_path))
    compute_file_hash(folder_path, files)
    @info("Files written to $folder_path folder.")
    return
end

"""
    get_result_variable(IS.results, Symbol, PSY.DataType)

Retrieve a specific variable dataframe from the results.

# Arguments
- `results::IS.Results`
- `name::Symbol`: The prefix for a type of variable or parameter
- `PSY.DataType`: The datatype of the variable from Power Systems

# Example
```julia
variable = get_result_variable(results, :ON, ThermalStandard)
```
"""

function get_result_variable(results::IS.Results, sym::Symbol, data_type::PSY.DataType)
    variable_name = encode_symbol(data_type, sym)
    if variable_name in keys(IS.get_variables(results))
        variable = IS.get_variables(results)[variable_name]
        return variable
    else
        @info "Variable $variable_name not found in results."
    end
end


text:
"""
    Vessel
        target::AbstractTargetingBehavior
        catchability::Catchability
        ξ::Real
        ϕ::Real

Describes the behavior of a fishing vessel, including its targeting strategy
(random or preferential), its catchability coefficient (which may vary across
space), and the parameters of the Tweedie distribution that govern its catch
distribution.
"""
struct Vessel{Tt, Tq, Tf}
    target::Tt
    catchability::Tq
    ξ::Tf
    ϕ::Tf
    function Vessel(target::Tt,
                    catchability::Tq,
                    ξ::Tf,
                    ϕ::Tf) where {Tt<:AbstractTargetingBehavior,
                                  Tq<:AbstractCatchability,
                                  Tf<:Real}
        new{Tt, Tq, Tf}(target, catchability, ξ, ϕ)
    end
end

## Vector of vessels needs to be abstract type for now; need to figure out small
## unions to type more concretely
struct Fleet{Tv, Te<:Integer}
    vessels::Vector{Tv}
    total_effort::Vector{Te}
    priority::Vector{Te}

    function Fleet(vessels::Vector{Tv},
                   total_effort::Vector{Te},
                   priority::Vector{Te}) where {Tv<:Vessel, Te<:Integer}
        length(vessels) == length(total_effort) == length(priority) ||
            throw(DimensionMismatch("Must have an effort and priority for each vessel"))
        priority = denserank(priority)
        new{Tv, Te}(vessels, total_effort, priority)
    end
end

function Fleet(vessels::Vector{Tv},
               total_effort::Vector{Te}) where {Tv<:Vessel, Te<:Integer}
    Fleet(vessels, total_effort, ones(eltype(total_effort), length(total_effort)))
end

getindex(F::Fleet, i) = F.vessels[i]
vessels(F::Fleet) = F.vessels
length(F::Fleet) = length(F.vessels)

"""
    Catch
        time::Ti
        loc_idx::Ti
        coordinates::Tuple{Tf, Tf}
        effort::Tf
        c::Tf

Holds catch and effort data from a given fleet's fishing effort. Does
not account for within-season depletion.
"""
struct Catch{Tf, Ti} <: Any
    time::Ti
    vessel_idx::Ti
    loc_idx::Ti
    coord1::Tf
    coord2::Tf
    effort::Tf
    catch_biomass::Tf

    function Catch(time::Ti, vessel_idx::Ti,
                   loc_idx::Ti, coord1::Tf, coord2::Tf,
                   effort::Tf, catch_biomass::Tf) where {Tf<:Real, Ti<:Integer}
        new{Tf, Ti}(time, vessel_idx, loc_idx, coord1, coord2, effort, catch_biomass)
    end
end
function Catch(time::Ti, vessel_idx::Ti,
               loc_idx::Ti, coordinates::Tuple{Tf, Tf},
               effort::Tf, catch_biomass::Tf) where {Tf<:Real, Ti<:Integer}
    Catch(time, vessel_idx, loc_idx, coordinates[1], coordinates[2], effort, catch_biomass)
end

"""
    fish!(P::PopState, V::Vessel, Ω::AbstractFisheryDomain, t::Integer)

Fish down a population state P. Returns a `Catch` object, mutates the
PopState in place. Limited to biomass available in a cell.
"""
function fish!(P::PopState,
               V::Vessel,
               Ω::AbstractFisheryDomain,
               t::Integer = 1,
               vessel_idx = 0)
    target_location = (l = target(V.target, Ω), t = t)
    μ = P.P[target_location.l] * V.catchability[target_location]
    if μ == 0
        catch_biomass = μ
    else
        catch_biomass = rand(CompoundPoissonGamma(μ, V.ξ, V.ϕ))
    end
    if catch_biomass > P.P[target_location.l]
        catch_biomass = P[target_location.l]
    end
    setindex!(P, P[target_location.l] - catch_biomass, target_location.l)
    C = Catch(t, vessel_idx, target_location.l, Ω.locs[target_location.l],
              convert(typeof(catch_biomass), 1), catch_biomass)
end

"""
    fish!(P::PopState, F::Fleet, Ω::AbstractFisheryDomain, t::Integer)

Fish down a population P using vessels in fleet F and their associated efforts.
Vessels harvest in random order and cells are sequentially depleted.
"""
function fish!(P::PopState{Tf},
               F::Fleet,
               Ω::AbstractFisheryDomain,
               t::Ti = 1) where {Tf<:Real, Ti<:Integer}
    effort_vec = order_effort(F)
    catch_record = Vector{Catch{Tf, Ti}}()
    for idx in effort_vec
        c = fish!(P, F[idx], Ω, t, idx)
        push!(catch_record, c)
    end
    # FIXME need a better way to deal with resetting stateful targeting each
    # year and dealing with changing targeting over time
    for vessel in vessels(F)
        reset!(vessel.target, P)
    end
    catch_record
end

"Return a vector with vessel indices ordered for fishing."
function order_effort(fleet::Fleet{Tv, Te}) where {Tv, Te}
    tot_eff = fleet.total_effort
    pri = fleet.priority
    vessel_idx = 1:length(fleet)
    eff_vvec = [fill(i, tot_eff[i]) for i in vessel_idx]

    eff = [Te[] for _ in unique(pri)]
    for (i, p) in enumerate(pri)
        append!(eff[p], eff_vvec[i])
    end
    shuffle!.(eff)
    reduce(vcat, eff)
end


text:
struct Game
    date :: Date
    opponent :: String
    team_spi :: Float64
    opp_spi :: Float64
    team_score :: Float64
    opp_score :: Float64
end

function load_games(games_df :: DataFrame)
    games = Dict{String, Vector{Game}}()
    for row ∈ eachrow(games_df[:, :team1])
        games[row[1, 1]] = Vector{Game}()
    end

    for row ∈ eachrow(games_df[:, [:date, :team1, :team2, :spi1, :spi2, :xg1, :xg2, :score1, :score2]])
        if !ismissing(row[:xg1])
            push!(games[row[:team1]], Game(row[:date], row[:team2], row[:spi1], row[:spi2], row[:xg1], row[:xg2]))
            push!(games[row[:team2]], Game(row[:date], row[:team1], row[:spi2], row[:spi1], row[:xg2], row[:xg1]))
        end
    end
    return games
end

discounted_score(g :: Game) = (g.team_score - g.opp_score) * (g.opp_spi / g.team_spi)

max_stars_repo_name: shalinkpatel/football-form

max_issues_repo_name: shalinkpatel/football-form

max_forks_repo_name: shalinkpatel/football-form


text:
using GitForge
const GF = GitForge

using HTTP: HTTP
using JSON2: JSON2
using Test: @test, @testset

function capture(f::Function)
    so = stdout
    r, w = redirect_stdout()
    out = @async read(r, String)
    result = f()
    redirect_stdout(so)
    close(w)
    return result, fetch(out)
end

struct TestPostProcessor <: GF.PostProcessor end
GF.postprocess(::TestPostProcessor, ::HTTP.Response, ::Type) = :foo

struct TestForge <: GF.Forge end
GF.base_url(::TestForge) = "https://httpbin.org"
GF.request_headers(::TestForge, ::Function) = ["Foo" => "Bar"]
GF.request_query(::TestForge, ::Function) = Dict("foo" => "bar")
GF.request_kwargs(::TestForge, ::Function) = Dict(:verbose => 1)
GF.postprocessor(::TestForge, ::Function) = TestPostProcessor()
GF.endpoint(::TestForge, ::typeof(get_user)) = GF.Endpoint(:GET, "/get")
GF.into(::TestForge, ::typeof(get_user)) = Symbol

@testset "GitForge.jl" begin
    f = TestForge()
    (val, resp), out = capture(() -> get_user(f))

    @testset "Basics" begin
        @test val === :foo
        @test resp isa HTTP.Response
    end

    @testset "Request options" begin
        body = JSON2.read(IOBuffer(resp.body))
        @test startswith(get(body, :url, ""), "https://httpbin.org")
        @test get(get(body, :headers, Dict()), :Foo, "") == "Bar"
        @test get(get(body, :args, Dict()), :foo, "") == "bar"
        @test occursin("/get", get(body, :url, ""))
        @test !isempty(out)
    end

    @testset "Per-call request options" begin
        (_, resp), out = capture() do
            get_user(
                f;
                query=Dict("a" => "b"),
                headers=["A" => "B"],
                request_opts=(; verbose=0),
            )
        end

        @test isempty(out)

        body = JSON2.read(IOBuffer(resp.body))
        @test haskey(body.headers, :Foo)
        @test get(body.headers, :A, "") == "B"
        @test haskey(get(body, :args, Dict()), :foo)
        @test get(get(body, :args, Dict()), :a, "") == "b"
    end
end


text:
using Documenter, BellScenario

DocMeta.setdocmeta!(BellScenario, :DocTestSetup, :(using BellScenario); recursive=true)

makedocs(;
    modules=[BellScenario],
    format=Documenter.HTML(
        assets = ["assets/custom.css"],
    ),
    pages=[
        "Home" => "index.md",
        "User Guide" => "user_guide.md",
        "BellScenario.jl" => [
            "Overview" => "BellScenario/overview.md",
            "Scenarios" => "BellScenario/scenarios.md",
            "Strategies" => "BellScenario/strategies.md",
            "Games" => "BellScenario/games.md",
            "Combinatorics" => "BellScenario/combinatorics.md",
        ],
        "LocalPolytope.jl" => [
            "Overview" => "LocalPolytope/overview.md",
            "Vertices" => "LocalPolytope/vertices.md",
            "Facets" => "LocalPolytope/facets.md",
            "Generators" => "LocalPolytope/generators.md",
            "Adjacency Decomposition" => "LocalPolytope/adjacency_decomposition.md",
            "Utilities" => "LocalPolytope/utils.md",
        ],
        "Nonlocality.jl" => [
            "Overview" => "Nonlocality/overview.md",
            "Optimize Measurements" => "Nonlocality/optimize_measurements.md",
        ],
        "Development Manual" => "development_manual.md",
    ],
    repo="https://github.com/ChitambarLab/BellScenario.jl/blob/{commit}{path}#L{line}",
    sitename="BellScenario.jl",
    authors="Brian Doolittle",
)

deploydocs(;
    repo="github.com/ChitambarLab/BellScenario.jl",
    devbranch = "main",
)


text:
# Uses method from Parviainen et al. 2020 to compute the sky-separation.

"""Compute the derivative components of the series expansion."""
function components(d::AbstractVector{T},h::T) where T<:AbstractFloat
    A = d[7]-d[1]; B = d[1]+d[7]; C = d[2]-d[6]
    D = d[2]+d[6]; E = d[5]-d[3]; F = d[3]+d[5]
    v = dot(SVector(A, C, E), COEFF.vd)/h
    a = dot(SVector(B, D, F, d[4]), COEFF.ad)/(h*h)
    j = dot(SVector(-A, -C, -E), COEFF.jd)/(h*h*h)
    s = dot(SVector(B, D, F , d[4]), COEFF.sd)/(h*h*h*h)
    return SVector(d[4],v,0.5*a,(1/6)*j,(1/24)*s)
end

"""Compute the impact parameter at t=tc using a series expansion about t0."""
function compute_impact_parameter(tc::T, t0::T, xc::AbstractVector{T}, yc::AbstractVector{T}) where T<:Real
    t = tc - t0
    ts = SVector(1.0, t, t*t, t*t*t, t*t*t*t)
    lx = dot(xc, ts)
    ly = dot(yc, ts)
    return sqrt(lx*lx + ly*ly)
end

function compute_impact_parameter(tc::T,t0::T,h::T,points::AbstractMatrix{T}) where T<:Real
    t = tc - t0
    ts = SVector(1.0,t,t*t,t*t*t,t*t*t*t)
    xc = components(points[:,1], h)
    yc = components(points[:,2], h)
    lx = dot(xc,ts)
    ly = dot(yc,ts)
    return sqrt(lx*lx + ly*ly)
end

function compute_impact_parameter!(tc::T, t0::T, xc, yc, dxc, dyc, grad) where T<:Real
    t = tc - t0
    ts = SVector(1.0, t, t*t, t*t*t, t*t*t*t)
    lx = dot(xc, ts)
    ly = dot(yc, ts)
    b = sqrt(lx*lx + ly*ly)

    dbdlx = lx/b; dbdly = ly/b;
    N = length(grad)
    for p in 1:N
        dlxdp = dot(dxc[p],ts)
        dlydp = dot(dyc[p],ts)
        grad[p] = dbdlx*dlxdp + dbdly*dlydp
    end

    return b
end


text:
s = SingleStock()
a = 100.00USD

d1 = Dates.today()
d2 = d1 + Dates.Day(120)

m = CoreModel(d1, a, 0.05, 0.0)

@test value(m, s) == a

@test value(m, Miletus.Zero()) == 0*a
@test value(m, Miletus.Receive(a)) == a

@test value(m, Pay(a)) == -a

@test value(m, Buy(s, a)) == 0USD
@test value(m, Sell(s, a)) == 0USD

@test value(m, Option(s)) == a

@test value(m, Forward(d2, s, 100.00USD)) ≈ 100.0USD * (1-exp(-0.05*120/365))

# assumes cost-of-carry
m = GeomBMModel(d1, a, 0.05, -0.1, 0.1)

@test value(m, Miletus.Zero()) == 0*a
@test value(m, Miletus.Receive(a)) == a

@test value(m, Pay(a)) == -a

@test value(m, Buy(s, a)) == 0USD
@test value(m, Sell(s, a)) == 0USD

@test value(m, Option(s)) == a

@test value(m, Scale(3, Option(s))) == 3*a

@test value(m, Forward(d2, s, 100.00USD)) ≈ 100.0USD * (exp(0.1*120/365)-exp(-0.05*120/365))

m = CoreForwardModel(d1,
                     [d1 + Dates.Day(30), d1 + Dates.Day(60), d1 + Dates.Day(120)],
                     [100.00USD, 110.00USD, 120.00USD],
                     0.05)

@test value(m, Forward(d2, s, 100.00USD)) ≈ (120.00USD - 100.0USD) * exp(-0.05*120/365)


text:
"""
    reduce_solution!(integrator::DDEIntegrator, tmax)

Remove time points of ODE solution of `integrator` up to time `tmax` that are not required
for calculation of DDE solution.
"""
@inline function reduce_solution!(integrator::DDEIntegrator, tmax)
    if integrator.saveiter < integrator.prev_idx
        # last time point of ODE solution that is known to be required for DDE solution
        @inbounds t_sol_prev = integrator.integrator.sol.t[integrator.saveiter]
        # next time point of ODE solution that might be required for DDE solution
        @inbounds t_sol = integrator.integrator.sol.t[integrator.saveiter + 1]

        @inbounds while t_sol < tmax && integrator.saveiter < integrator.prev_idx &&
            # do not handle time point if it is required for interpolation of tmax
            (integrator.saveiter + 2 > integrator.integrator.saveiter ||
             tmax >= integrator.integrator.sol.t[integrator.saveiter + 2])

            needed = false

            while !isempty(integrator.saveat) && top(integrator.saveat) <= t_sol
                # do not remove time point if it is upper bound of a required
                # interpolation interval
                if t_sol_prev < top(integrator.saveat)
                    needed = true
                end

                # remove passed time points
                pop!(integrator.saveat)
            end

            # do not remove time point if it is lower bound of a required interpolation
            # interval
            if !needed && integrator.saveiter + 1 < integrator.integrator.saveiter
                t_sol_next = integrator.integrator.sol.t[integrator.saveiter + 2]

                needed = !isempty(integrator.saveat) && top(integrator.saveat) < t_sol_next
            end

            if !needed
                # delete not required time point, function value, and interpolation data
                deleteat!(integrator.integrator.sol.t, integrator.saveiter + 1)
                deleteat!(integrator.integrator.sol.u, integrator.saveiter + 1)
                deleteat!(integrator.integrator.sol.k, integrator.saveiter + 1)

                # correct other arrays
                if typeof(integrator.integrator.alg) <: OrdinaryDiffEq.OrdinaryDiffEqCompositeAlgorithm
                    deleteat!(integrator.integrator.alg_choice, integrator.saveiter + 1)
                end

                # update counter of saved time points
                integrator.integrator.saveiter -= 1
                integrator.integrator.saveiter_dense -= 1

                # update indices of tprev and u(tprev), and t and u(t) in solution
                integrator.prev_idx -= 1
                integrator.prev2_idx -= 1

                # update time point of loop
                t_sol = integrator.integrator.sol.t[integrator.saveiter + 1]
            else
                # increase counter of time points that are required for DDE solution
                integrator.saveiter += 1

                # update time points of loop
                t_sol_prev = t_sol
                t_sol = integrator.integrator.sol.t[integrator.saveiter + 1]
            end
        end
    end
end

"""
    build_solution_array(integrator::DDEIntegrator)

Create a `DiffEqArray` of the time points and values that form the solution of `integrator`.
"""
function build_solution_array(integrator::DDEIntegrator)
    if integrator.opts.save_everystep && isempty(integrator.opts.saveat)
        # use solution of ODE integrator if no additional time points provided
        if integrator.opts.save_start
            t = integrator.sol.t
            if integrator.opts.save_idxs === nothing
                u = integrator.sol.u
            else
                u = [u[integrator.opts.save_idxs] for u in integrator.sol.u]
            end
        else # remove initial time point
            t = integrator.sol.t[2:end]
            if integrator.opts.save_idxs === nothing
                u = integrator.sol.u[2:end]
            else
                u = [u[integrator.opts.save_idxs] for u in @view(integrator.sol.u[2:end])]
            end
        end
    else
        # calculate number of additional time points of solution
        saveat_length = count(
            x -> integrator.tdir * x < integrator.tdir * integrator.sol.t[end],
            integrator.opts.saveat.valtree)

        # create vectors of time points and corresponding values which form final solution
        n = integrator.opts.save_everystep ? length(integrator.sol.t) : 2
        n += saveat_length
        integrator.opts.save_start || (n -= 1)
        t = Vector{typeof(integrator.t)}(undef,n)
        if integrator.opts.save_idxs === nothing
          u = Vector{typeof(integrator.u)}(undef,n)
        else
          u = Vector{typeof(integrator.u[integrator.opts.save_idxs])}(undef,n)
        end

        # output initial time point if desired
        write_idx = 1 # next index of solution to write to
        if integrator.opts.save_start
            t[1] = integrator.sol.t[1]
            if integrator.opts.save_idxs === nothing
                u[1] = integrator.sol.u[1]
            else
                u[1] = integrator.sol.u[1][integrator.opts.save_idxs]
            end
            write_idx = 2
        end

        # merge additional time points and time points of solution of ODE integrator
        # both data structures are already sorted
        if integrator.opts.save_everystep
            sol_idx = 2 # next index of solution of ODE integrator to read from
            sol_length = length(integrator.sol.t)

            @inbounds while saveat_length > 0 && sol_idx < sol_length
                if integrator.tdir * top(integrator.opts.saveat) <
                    integrator.tdir * integrator.sol.t[sol_idx]

                    # copy additional time points and calculate corresponding values
                    t[write_idx] = pop!(integrator.opts.saveat)
                    u[write_idx] = integrator.sol(t[write_idx];
                                                  idxs=integrator.opts.save_idxs)

                    saveat_length -= 1
                else
                    # copy time points and values of solution of ODE integrator
                    t[write_idx] = integrator.sol.t[sol_idx]
                    if integrator.opts.save_idxs === nothing
                        u[write_idx] = integrator.sol.u[sol_idx]
                    else
                        u[write_idx] = integrator.sol.u[sol_idx][integrator.opts.save_idxs]
                    end

                    sol_idx += 1
                end

                write_idx += 1
            end

            # copy remaining time points of solution of ODE integrator
            # except of final time point
            copyto!(t, write_idx, integrator.sol.t, sol_idx, sol_length - sol_idx)
            if integrator.opts.save_idxs === nothing
                copyto!(u, write_idx, integrator.sol.u, sol_idx, sol_length - sol_idx)
            else
                copyto!(u, write_idx,
                        (u[integrator.opts.save_indxs] for u in integrator.sol.u),
                        sol_idx, sol_length - sol_idx)
            end

            write_idx += sol_length - sol_idx
        end

        # copy remaining additional time points and calculate corresponding values
        @inbounds while saveat_length > 0
            t[write_idx] = pop!(integrator.opts.saveat)
            u[write_idx] = integrator.sol(t[write_idx]; idxs=integrator.opts.save_idxs)

            saveat_length -= 1
            write_idx += 1
        end

        # always output final time point
        t[end] = integrator.sol.t[end]
        if integrator.opts.save_idxs === nothing
            u[end] = integrator.sol.u[end]
        else
            u[end] = integrator.sol.u[end][integrator.opts.save_idxs]
        end
    end

    DiffEqArray(u, t)
end

"""
    build_solution_interpolation(integrator::DDEIntegrator, sol::DiffEqArray)

Create interpolation data to solution of `integrator`, which is formed by time points and
values in `sol`.
"""
function build_solution_interpolation(integrator::DDEIntegrator, sol::DiffEqArray)
    if integrator.opts.dense
        if integrator.opts.save_idxs === nothing
            integrator.sol.interp
        else # update interpolation data if only a subset of indices is returned
            if typeof(integrator.alg) <: OrdinaryDiffEq.OrdinaryDiffEqCompositeAlgorithm
                OrdinaryDiffEq.CompositeInterpolationData(
                    integrator.sol.interp.f, [u[integrator.opts.save_idxs] for u in
                                              integrator.sol.interp.timeseries],
                    integrator.sol.interp.ts, [[k[integrator.opts.save_idxs] for k in ks]
                                               for ks in integrator.sol.interp.ks],
                    integrator.sol.interp.alg_choice,
                    true, integrator.sol.interp.cache)
            else
                OrdinaryDiffEq.InterpolationData(
                    integrator.sol.interp.f, [u[integrator.opts.save_idxs] for u in
                                              integrator.sol.interp.timeseries],
                    integrator.sol.interp.ts, [[k[integrator.opts.save_idxs] for k in ks]
                                               for ks in integrator.sol.interp.ks],
                    true, integrator.sol.interp.cache)
            end
        end
    else # create not dense interpolation data if desired
        if typeof(integrator.alg) <: OrdinaryDiffEq.OrdinaryDiffEqCompositeAlgorithm
            OrdinaryDiffEq.CompositeInterpolationData(
                integrator.sol.interp.f, sol.u, sol.t, typeof(integrator.sol.k)(undef, 0),
                Int[], false, integrator.sol.interp.cache)
        else
            OrdinaryDiffEq.InterpolationData(
                integrator.sol.interp.f, sol.u, sol.t, typeof(integrator.sol.k)(undef, 0),
                false, integrator.sol.interp.cache)
        end
    end
end

"""
    advance_ode_integrator!(integrator::DDEIntegrator)

Advance ODE integrator of `integrator` to next time interval, values and complete
interpolation data of `integrator`.
"""
function advance_ode_integrator!(integrator::DDEIntegrator)
    # algorithm only works if current time of DDE integrator equals final time point
    # of solution
    integrator.t != integrator.sol.t[end] && error("cannot advance ODE integrator")

    # complete interpolation data of DDE integrator for time interval [t, t+dt]
    # and copy it to ODE integrator
    # has to be done before updates to ODE integrator, otherwise history function
    # is incorrect
    if typeof(integrator.cache) <: OrdinaryDiffEq.CompositeCache
        addsteps!(integrator.k, integrator.t, integrator.uprev,
                                     integrator.u, integrator.dt, integrator.f,
                                     integrator.p,
                                     integrator.cache.caches[integrator.cache.current],
                                     false, true, true)
    else
        addsteps!(integrator.k, integrator.t, integrator.uprev,
                                     integrator.u, integrator.dt, integrator.f,
                                     integrator.p,
                                     integrator.cache, false, true,
                                     true)
    end
    recursivecopy!(integrator.integrator.k, integrator.k)

    # move ODE integrator to interval [t, t+dt]
    integrator.integrator.t = integrator.t + integrator.dt
    integrator.integrator.tprev = integrator.t
    integrator.integrator.dt = integrator.dt
    if typeof(integrator.alg) <: OrdinaryDiffEq.OrdinaryDiffEqCompositeAlgorithm
      integrator.integrator.cache.current = integrator.cache.current
    end
    if isinplace(integrator.sol.prob)
        recursivecopy!(integrator.integrator.u, integrator.u)
    else
        integrator.integrator.u = integrator.u
    end

    # u(t) is not modified hence we do not have to copy it
    integrator.integrator.uprev = integrator.sol.u[end]

    # update prev_idx to index of t and u(t) in solution
    integrator.prev_idx = length(integrator.sol.t)
end

#=
Dealing with discontinuities

If we hit a discontinuity (this is checked in `apply_step!`), then we remove the
discontinuity, additional discontinuities at the current time point (if present), and
maybe also discontinuities and time stops coming shortly after the current time point
in `handle_discontinuities!`. The order of the discontinuity at the current time point is
defined as the lowest order of all these discontinuities.

If the problem is not neutral, we will only add additional discontinuities if
this order is less or equal to the order of the algorithm in
`add_next_discontinuities!`. If we add discontinuities, we add discontinuities
of the next order caused by constant lags (these we can calculate explicitly and
just add them to `d_discontinuities` and `tstops`) and we add the current
discontinuity to `tracked_discontinuities` which is the array of old
discontinuities that are checked by a `DiscontinuityCallback` (if existent).
=#

"""
    handle_discontinuities!(integrator::DDEIntegrator)

Handle discontinuities at the current time point of `integrator`.
"""
function handle_discontinuities!(integrator::DDEIntegrator)
    # remove all discontinuities at current time point and calculate minimal order
    # of these discontinuities
    d = pop!(integrator.opts.d_discontinuities)
    order = d.order
    while !isempty(integrator.opts.d_discontinuities) &&
        top(integrator.opts.d_discontinuities) == integrator.t

        d2 = pop!(integrator.opts.d_discontinuities)
        order = min(order, d2.order)
    end

    # remove all discontinuities close to the current time point as well and
    # calculate minimal order of these discontinuities
    # integrator.EEst has unitless type of integrator.t
    if typeof(integrator.EEst) <: AbstractFloat
        maxΔt = 10eps(integrator.t)

        while !isempty(integrator.opts.d_discontinuities) &&
            abs(top(integrator.opts.d_discontinuities).t - integrator.t) < maxΔt

            d2 = pop!(integrator.opts.d_discontinuities)
            order = min(order, d2.order)
        end

        # also remove all corresponding time stops
        while !isempty(integrator.opts.tstops) &&
            abs(top(integrator.opts.tstops) - integrator.t) < maxΔt

            pop!(integrator.opts.tstops)
        end
    end

    # add discontinuities of next order to integrator
    add_next_discontinuities!(integrator, order)
end

"""
    add_next_discontinuities!(integrator::DDEIntegrator, order[, t=integrator.t])

Add discontinuities of next order that are propagated from discontinuity of
order `order` at time `t` in `integrator`, but only if `order` is less or equal
than the order of the applied method or the problem is neutral.

Discontinuities caused by constant delays are immediately calculated, and
discontinuities caused by dependent delays are tracked by a callback.
"""
function add_next_discontinuities!(integrator, order, t=integrator.t)
    # obtain delays
    neutral = integrator.sol.prob.neutral
    constant_lags = integrator.sol.prob.constant_lags

    # only track discontinuities up to order of the applied method
    order > alg_maximum_order(integrator.alg) && !neutral && return nothing

    # discontinuities caused by constant lags
    if constant_lags !== nothing
        maxlag = abs(integrator.sol.prob.tspan[end] - t)
        for lag in constant_lags
            if abs(lag) < maxlag
                # calculate discontinuity and add it to heap of discontinuities and time stops
                d = Discontinuity(t + lag, order + 1)
                push!(integrator.opts.d_discontinuities, d)
                push!(integrator.opts.tstops, d.t)
            end
        end
    end

    # track propagated discontinuities with callback
    push!(integrator.tracked_discontinuities, Discontinuity(t, order))
end


text:
# function to check if the current character represents a mine or not
ismine(ch) = (ch == '*' ? 1 : 0)

function annotate(minefield)

    # ans will store our final output of the program
    ans = []

    # i represents row number
    for i in 1 : length(minefield)

        # ans_row string will be an element of our ans
        ans_row = ""

        # j represents column number
        for j in 1 : length(minefield[i])

            # check if it is a mine or not
            if minefield[i][j] != '*'

                # calculate the boundary indices of neighbourhood
                r1 = max(1, i - 1)
                r2 = min(i + 1, length(minefield))
                c1 = max(1, j - 1)
                c2 = min(j + 1, length(minefield[i]))

                # count will store the number of mines in the neighbourhood
                count = 0
                for x in r1:r2, y in c1:c2
                    count += ismine(minefield[x][y])
                end

                # append space char if no mine is found else append count
                if count == 0
                    ans_row *= ' '
                else
                    ans_row *= string(count)
                end

            # if current cell represents mine then append it
            else
                ans_row *= '*'
            end
        end

        # push the calculated row_string to our ans
        push!(ans, ans_row)
    end

    return ans
end


text:
using CTMCEnsemble
using Base.Test

@testset "example" begin
    preds = [([0.5, 0.5], [1, 2]), ([0.5, 0.5], [2, 3])]
    G = CTMCEnsemble.build(preds)
    @test CTMCEnsemble.stationdist(G) ≈ [1/3, 1/3, 1/3]

    @test [1/3, 1/3, 1/3] ≈ average(preds)
    @test [1/4, 1/2, 1/4] ≈ average(preds, multiplicity=false)
    @test [1/3, 1/3, 1/3] ≈ product(preds)
    @test [2/5, 1/5, 2/5] ≈ product(preds, multiplicity=false)
    @test [1/3, 1/3, 1/3] ≈ powermethod(preds)
    @test [1/3, 1/3, 1/3] ≈ svdmethod(preds)
    @test [1/3, 1/3, 1/3] ≈ ctmc(preds)
end

@testset "single" begin
    # single
    for i = 1:5
        d = 1000
        p = normalize!(rand(d), 1)
        preds = [(p, 1:d)]

        @test p ≈ average(preds)
        @test p ≈ average(preds, multiplicity=false)
        @test p ≈ product(preds)
        @test p ≈ product(preds, multiplicity=false)
        @test p ≈ powermethod(preds)
        @test p ≈ svdmethod(preds)
        @test p ≈ ctmc(preds)
    end
end

@testset "multiple" begin
    newv() = normalize!(rand(300), 1)

    for i = 1:5
        d = 600
        preds = [
            (newv(),     1:300),
            (newv(),   201:500),
            (newv(), ((401:700) - 1) .% 600 .+ 1),
        ]

        v_powermethod = powermethod(preds, maxiter=64)
        v_svdmethod   = svdmethod(preds)
        v_ctmc        = ctmc(preds)

        @test v_svdmethod ≈ v_ctmc
        @test_approx_eq_eps v_svdmethod v_powermethod 1/d
    end
end

@testset "love bufan" begin
    @test 1 == 1
end


text:
CONFIG = Dict(
    "USER_KEY" => "YOUR_USER_KEY",
    "API_TOKEN" => "YOUR_API_TOKEN",
    "TITLE" => "YOUR_TITLE",
    "PRIORITY" => 0
)