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

using Markdown
using InteractiveUtils

# ╔═╡ a57a7e52-f314-4bd5-8a23-644085702ed6
using PlutoUI; PlutoUI.TableOfContents(title="目录")

# ╔═╡ 35cda3c5-f6a9-4be5-86c5-e265a417f5fc
using EDF # EDF.read() 读入 bdf+ 数据

# ╔═╡ 4a805a3f-51d5-4271-8353-a6d10f1cda74
using DataFrames

# ╔═╡ 1949fef3-1c95-4c3b-85d7-3ee7c20908c3
md"""
## 载入包
"""

# ╔═╡ d57abd01-4081-48ce-b1eb-9515d011665b
import Pipe: @pipe

# ╔═╡ 91c41e82-f69f-4e34-b05d-b981a7c0009e
md"""
## 设定当前文件夹
"""

# ╔═╡ 517a6404-4448-40d8-b867-666a03335618
cd(@__DIR__) # 把当下 notebook 所在文件夹定义为当前文件夹

# ╔═╡ 6116fed2-1cef-443b-b0b8-81de40037769
md"""
## 分析过程
"""

# ╔═╡ 9785790c-e5bf-42e9-a687-b931ecb2eff6
md"""
## 定义函数
"""

# ╔═╡ bfe3b2ca-becb-4c28-8f9d-a99f75d8c4f0
md"""
- 定义函数 `base_line` 提取被试的 `baseline` 时间点
"""

# ╔═╡ c5341c98-0773-43ba-a9b7-f369f920cb40
md"""
- 定义函数 `retrieve_triggers()` 把函数 `read_bdfplus_event()` 和 `retrieve_connective_pause()` 组合起来
"""

# ╔═╡ 9ee7b8a6-14c1-447a-9311-19584a9ee814
md"""
- 定义函数 `retrieve_bdfplus_event` 以从 `evt.bdf` 中提取事件信息
"""

# ╔═╡ ed357abb-551e-4bc6-8d30-1b35bc2af390
function read_bdfplus_event(evtbdf)
	# 读入数据
	evtfile = EDF.read(evtbdf)

	# 提取事件序号: 0-1080, 共 1081 个
	index = [evntpair[1].onset_in_seconds for evntpair in evtfile.signals[2].records]

	# 提取事件触发时间， 单位是秒
	onset = [evntpair[2].onset_in_seconds for evntpair in evtfile.signals[2].records]

	# 提取事件编码
	annotation =
		[evntpair[2].annotations[1] for evntpair in evtfile.signals[2].records]

	# 汇总成一个数据框
	return DataFrame(Number = index, Latency = onset, Type = annotation)
end

# ╔═╡ f07abe75-7079-4f18-8f6e-7bb46aa90da7
base_line(bdf) = @pipe read_bdfplus_event(bdf) |>
	filter(:Type => ==("Annotation"), _)       |>
	transform(_, :Latency => ByRow(x -> x * 1000) => :Latency)

# ╔═╡ a423d6fa-2243-424b-ac2d-42278d01053d
sbj = base_line("08/evt.bdf")

# ╔═╡ 04838004-a30e-4bfa-81e4-f4b2f222568f
md"""
- 定义 `trigger` 信息中的编码常数
"""

# ╔═╡ 6f6f8f19-f0ea-4193-bdd8-5f4f2ee627ac
begin
	const Conditions      = ["211", "212", "221", "222"]
	const TrialNumber     = lpad.(1:180, 3, '0')
	const Noun2Onset      = "202"
	const ConnectiveOnset = "203"
	const Noun3Onset      = "204"
	const TrialOffset     = "255"
end;

# ╔═╡ fe2ab09d-68c6-4cef-9d03-a05d67856b76
PlutoUI.LocalResource("TriggerInformation.png")

# ╔═╡ e926ed83-c2b0-40d7-9f27-db34fa135a7d
md"""
- 定义函数 `retrieve_connective_pause` 以从读取的`事件文件`中提取的连接词和停顿信息。
"""

# ╔═╡ a141d9f2-b6ba-49b9-8456-6ea330638ae5
function retrieve_connective_pause(dt)
 	@pipe dt                                                                |>

	# 添加 `TrialOnset` 列, 表示该行是否处于试次开始, 注： `Ref` 用法参见 `?∈`
	transform(_, :Type => ByRow(x -> x .∈ Ref(Conditions) ) => :TrialOnset) |>

	# 添加 `TrialCount` 列顺序累加试次的个数, 接下来用于分组 `groupby`
	transform(_, :TrialOnset => (cumsum) => :TrialCount)                    |>

	# 去掉 'Annotation' 行
	filter!(:Type => !=("Annotation"), _)                                   |>

	# 把 `Type` 和 `Number` 列修正成长度相同, 即 3个字符和4个字符, 用 `0` 填充
	transform(_,
		:Type   => ByRow(x -> lpad(x, 3, '0'))              => :Type,
		:Number => ByRow(x -> lpad(string(Int(x)), 4, '0')) => :Number
	)                                                                       |>

	# 用 `TrialCount` 列把数据框分成 180 个子数据框
	groupby(_, :TrialCount)                                                 |>

	# 在每一个子数据中依据 `Type` 列添加 9 个新列
	#   `Trial`:      用 `Type` 列的第二个数据填充;
	#   `Connective`: 由 `Type` 列第一个数据的第二位数决定: 1 -> And; 2 -> Or
	#   `Pause`:      由 `Type` 列第一个数据的第三位数决定: 1 -> 000ms; 2 -> 200ms
	#   `Sentence_Onset` 等 triggers 出现时间
	transform(_,
		:Type => (x -> x[2]) => :Trial,
		:Type => (x -> SubString(x[1], 2, 2) == "1" ? "And" : "Or") => :Connective,
		:Type => (x -> SubString(x[1], 3, 3) == "1" ? "000ms" : "200ms") => :Pause,
		:Type => ByRow(∈(Conditions))       => :Sentence_Onset,
		:Type => ByRow(∈(TrialNumber))      => :Verb_Onset,
		:Type => ByRow(==(Noun2Onset))      => :Noun2_Onset,
		:Type => ByRow(==(ConnectiveOnset)) => :Connective_Onset,
		:Type => ByRow(==(Noun3Onset))      => :Noun3_Onset,
		:Type => ByRow(==(TrialOffset))     => :Sentence_Offset
	)                                                                     |>

	# 删掉 `TrialCount` 列
	select(_, Not(:TrialCount))                                           |>

	# 把秒转化成毫秒
	transform(_, :Latency => ByRow(x -> x * 1000) => :Latency)

end

# ╔═╡ a810fb79-edfa-465e-bbbd-564290d2790d
retrieve_triggers(bdf) = @pipe read_bdfplus_event(bdf) |> retrieve_connective_pause(_)

# ╔═╡ 4a6348de-c4e9-4d3c-bdee-7d34eaf9ce40
evts = retrieve_triggers("08/evt.bdf")

# ╔═╡ dc73980b-e532-4f2a-8bad-a41b8c58b0c7
describe(evts)

# ╔═╡ 6f3df1c5-71d1-4faf-9c6c-3c3130da40ff
md"""
## 数据结构说明
"""

# ╔═╡ bb8626e9-bb55-47ca-9372-30f27a58fe4d
evtfile = EDF.read("08/evt.bdf")

# ╔═╡ 54efce69-939d-479c-bcf9-60d1ec9aa776
propertynames(evtfile)        # 数据记录了三个字段: `io`, `header`, `signals`

# ╔═╡ 725d9991-7e5a-4d07-a963-d822ddb3ad58
evtfile.io                    # 原始文件

# ╔═╡ 5b1423fc-7f8f-40d7-b625-b83d4f16cd60
evtfile.header               # 文件信息， 包括被试信息、 采集日期等

# ╔═╡ c0ab84d4-841c-4029-83d5-36d73a48a3b2
signals = evtfile.signals    # 信号信息： 两个 channel, 实际有用的是 `AnnotationalSignal`

# ╔═╡ f731b484-b506-4441-aff4-44a5a0b25240
event_signal = signals[2]    # 提取第二 channel 的信息

# ╔═╡ 6fbf4365-8a84-40ac-9ab2-895fc6aa24d6
typeof(event_signal)        # 查看第二个 channel 的信息存储类型

# ╔═╡ a4fb8c4d-1f86-4715-b075-d17a7f29efda
propertynames(event_signal) # 类型为 `AnnotationsSignal` 的对象有两个字段， `:samples_per_second` - 采样率, `:records` - 事件信息

# ╔═╡ 1ea979b9-644d-4b91-af79-13d9206375ca
evtss = event_signal.records # 实际记录到的事件

# ╔═╡ 3ab3caff-bd9f-4a40-953b-9817e4aa5cdd
length(evtss)                # 共 1081 个事件

# ╔═╡ 6c34d2ee-0d21-4b78-b32f-bcf5cb15ae95
event_first_pair = evtss[2]  # 每个事件由两个`时间事件列表`组成, `TimestampedAnnotationList`

# ╔═╡ d14d6f11-ed68-4307-a746-32685422f88a
event_two = event_first_pair[2]

# ╔═╡ f98cfe5b-afbb-49c8-82ef-87965c395d3e
propertynames(event_two)       # 每个`时间事件列表`有三个字段： `onset_in_seconds` - 事件开始时间 (s), `duration_in_seconds` - 事件持续时间, `:annotations` - 事件文字信息

# ╔═╡ 2849dca1-5886-4c77-90e0-3b72b8a2901c
event_two.onset_in_seconds

# ╔═╡ 9a278aee-b22f-4b95-8f5d-595ed89ef953
event_two.duration_in_seconds # 事件持续时间默认为零

# ╔═╡ ad91ea6e-c626-4b97-9622-2986aaf99c6a
annot = event_two.annotations # 事件文字信息

# ╔═╡ 2dc037d5-972a-46c1-b9cb-fb2658741532
typeof(annot)                 # 事件文字信息是一个`数组`

# ╔═╡ 6aa6c49c-4657-4fb1-8cf1-7604ef7a8e19
annot[1]                      # 把文字信息从数组中提取出来

# ╔═╡ 7aa5c780-e0df-4273-b41c-299f6f450ae2
annot[]                       # 如果数组中只有一个值， 可以这样提取

# ╔═╡ 0d3f8d18-ec69-48cb-a30f-731673388059
md"""
## `EDF.File` 帮助文档
"""

# ╔═╡ d7bfd9b4-9fb6-4fa3-9725-5feb26c8ee78
@doc EDF.File

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
EDF = "ccffbfc1-f56e-50fb-a33b-53d1781b2825"
Pipe = "b98c9c47-44ae-5843-9183-064241ee97a0"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
DataFrames = "~1.2.2"
EDF = "~0.7.2"
Pipe = "~1.3.0"
PlutoUI = "~0.7.18"
"""

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

[[AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "0ec322186e078db08ea3e7da5b8b2885c099b393"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.0"

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

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

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

[[BitIntegers]]
deps = ["Random"]
git-tree-sha1 = "f50b5a99aa6ff9db7bf51255b5c21c8bc871ad54"
uuid = "c3b6d118-76ef-56ca-8cc7-ebb389d030a1"
version = "0.2.5"

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

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

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

[[EDF]]
deps = ["BitIntegers", "Dates", "Printf"]
git-tree-sha1 = "bac5d5d2e936115a4c838d3a93a26f892acd3c48"
uuid = "ccffbfc1-f56e-50fb-a33b-53d1781b2825"
version = "0.7.2"

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

[[Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[HypertextLiteral]]
git-tree-sha1 = "5efcf53d798efede8fee5b2c8b09284be359bf24"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.2"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

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

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

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

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

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

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

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

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

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

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

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

[[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 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

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

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

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

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

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

[[Pipe]]
git-tree-sha1 = "6842804e7867b115ca9de748a0cf6b364523c16d"
uuid = "b98c9c47-44ae-5843-9183-064241ee97a0"
version = "1.3.0"

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

[[PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "57312c7ecad39566319ccf5aa717a20788eb8c1f"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.18"

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

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

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

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

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

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

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

[[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 = "fed34d0e71b91734bf0a7e10eb1bb05296ddbcd0"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.0"

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

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

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

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

[[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:
# ╟─a57a7e52-f314-4bd5-8a23-644085702ed6
# ╟─1949fef3-1c95-4c3b-85d7-3ee7c20908c3
# ╠═35cda3c5-f6a9-4be5-86c5-e265a417f5fc
# ╠═4a805a3f-51d5-4271-8353-a6d10f1cda74
# ╠═d57abd01-4081-48ce-b1eb-9515d011665b
# ╟─91c41e82-f69f-4e34-b05d-b981a7c0009e
# ╠═517a6404-4448-40d8-b867-666a03335618
# ╟─6116fed2-1cef-443b-b0b8-81de40037769
# ╠═4a6348de-c4e9-4d3c-bdee-7d34eaf9ce40
# ╠═dc73980b-e532-4f2a-8bad-a41b8c58b0c7
# ╠═a423d6fa-2243-424b-ac2d-42278d01053d
# ╟─9785790c-e5bf-42e9-a687-b931ecb2eff6
# ╟─bfe3b2ca-becb-4c28-8f9d-a99f75d8c4f0
# ╠═f07abe75-7079-4f18-8f6e-7bb46aa90da7
# ╟─c5341c98-0773-43ba-a9b7-f369f920cb40
# ╠═a810fb79-edfa-465e-bbbd-564290d2790d
# ╟─9ee7b8a6-14c1-447a-9311-19584a9ee814
# ╠═ed357abb-551e-4bc6-8d30-1b35bc2af390
# ╟─04838004-a30e-4bfa-81e4-f4b2f222568f
# ╠═6f6f8f19-f0ea-4193-bdd8-5f4f2ee627ac
# ╟─fe2ab09d-68c6-4cef-9d03-a05d67856b76
# ╟─e926ed83-c2b0-40d7-9f27-db34fa135a7d
# ╠═a141d9f2-b6ba-49b9-8456-6ea330638ae5
# ╟─6f3df1c5-71d1-4faf-9c6c-3c3130da40ff
# ╠═bb8626e9-bb55-47ca-9372-30f27a58fe4d
# ╠═54efce69-939d-479c-bcf9-60d1ec9aa776
# ╠═725d9991-7e5a-4d07-a963-d822ddb3ad58
# ╠═5b1423fc-7f8f-40d7-b625-b83d4f16cd60
# ╠═c0ab84d4-841c-4029-83d5-36d73a48a3b2
# ╠═f731b484-b506-4441-aff4-44a5a0b25240
# ╠═6fbf4365-8a84-40ac-9ab2-895fc6aa24d6
# ╠═a4fb8c4d-1f86-4715-b075-d17a7f29efda
# ╠═1ea979b9-644d-4b91-af79-13d9206375ca
# ╠═3ab3caff-bd9f-4a40-953b-9817e4aa5cdd
# ╠═6c34d2ee-0d21-4b78-b32f-bcf5cb15ae95
# ╠═d14d6f11-ed68-4307-a746-32685422f88a
# ╠═f98cfe5b-afbb-49c8-82ef-87965c395d3e
# ╠═2849dca1-5886-4c77-90e0-3b72b8a2901c
# ╠═9a278aee-b22f-4b95-8f5d-595ed89ef953
# ╠═ad91ea6e-c626-4b97-9622-2986aaf99c6a
# ╠═2dc037d5-972a-46c1-b9cb-fb2658741532
# ╠═6aa6c49c-4657-4fb1-8cf1-7604ef7a8e19
# ╠═7aa5c780-e0df-4273-b41c-299f6f450ae2
# ╟─0d3f8d18-ec69-48cb-a30f-731673388059
# ╠═d7bfd9b4-9fb6-4fa3-9725-5feb26c8ee78
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
#    Bessel functions of the first kind of order zero and one
#                       besselj0, besselj1
#
#    Calculation of besselj0 is done in three branches using polynomial approximations
#
#    Branch 1: x <= pi/2
#              besselj0 is calculated using a 9 term, even minimax polynomial
#
#    Branch 2: pi/2 < x < 26.0
#              besselj0 is calculated by one of 16 different degree 13 minimax polynomials
#       Each polynomial is an expansion around either a root or extrema of the besselj0.
#       This ensures accuracy near the roots. Method taken from [2]
#
#   Branch 3: x >= 26.0
#              besselj0 = sqrt(2/(pi*x))*beta(x)*(cos(x - pi/4 - alpha(x))
#   See modified expansions given in [2]. Exact coefficients are used.
#
#   Calculation of besselj1 is done in a similar way as besselj0.
#   See [2] for details on similarities.
#
# [1] https://github.com/deepmind/torch-cephes
# [2] Harrison, John. "Fast and accurate Bessel function computation."
#     2009 19th IEEE Symposium on Computer Arithmetic. IEEE, 2009.
#

"""
    besselj0(x::T) where T <: Union{Float32, Float64}

Bessel function of the first kind of order zero, ``J_0(x)``.
"""
function besselj0(x::Float64)
    T = Float64
    x = abs(x)

    if x < 26.0
        x < pi/2 && return evalpoly(x * x, J0_POLY_PIO2(T))
        n = unsafe_trunc(Int, TWOOPI(T) * x)
        root = @inbounds J0_ROOTS(T)[n]
        r = x - root[1] - root[2]
        return evalpoly(r, @inbounds J0_POLYS(T)[n])
    else
        xinv = inv(x)
        iszero(xinv) && return zero(T)
        x2 = xinv * xinv

        if x < 120.0
            p1 = (one(T), -1/16, 53/512, -4447/8192, 3066403/524288, -896631415/8388608, 796754802993/268435456, -500528959023471/4294967296)
            q1 = (-1/8, 25/384, -1073/5120, 375733/229376, -55384775/2359296, 24713030909/46137344, -7780757249041/436207616)
            p = evalpoly(x2, p1)
            q = evalpoly(x2, q1)
        else
            p2 = (one(T), -1/16, 53/512, -4447/8192)
            q2 = (-1/8, 25/384, -1073/5120, 375733/229376)
            p = evalpoly(x2, p2)
            q = evalpoly(x2, q2)
        end

        a = SQ2OPI(T) * sqrt(xinv) * p
        xn = muladd(xinv, q, -PIO4(T))

        # the following computes b = cos(x + xn) more accurately
        # see src/misc.jl
        b = cos_sum(x, xn)
        return a * b
    end
end
function besselj0(x::Float32)
    T = Float32
    x = abs(x)

    if x <= 2.0f0
        z = x * x
        if x < 1.0f-3
            return 1.0f0 - 0.25f0 * z
        end
        DR1 = 5.78318596294678452118f0
        p = (z - DR1) * evalpoly(z, JP_j0(T))
        return p
    else
        q = inv(x)
        iszero(q) && return zero(T)
        w = sqrt(q)
        p = w * evalpoly(q, MO_j0(T))
        w = q * q
        xn = q * evalpoly(w, PH_j0(T)) - PIO4(Float32)
        p = p * cos(xn + x)
        return p
    end
end

"""
    besselj1(x::T) where T <: Union{Float32, Float64}

Bessel function of the first kind of order one, ``J_1(x)``.
"""
function besselj1(x::Float64)
    T = Float64
    s = sign(x)
    x = abs(x)

    if x <= 26.0
        x <= pi/2 && return x * evalpoly(x * x, J1_POLY_PIO2(T))
        n = unsafe_trunc(Int, TWOOPI(T) * x)
        root = @inbounds J1_ROOTS(T)[n]
        r = x - root[1] - root[2]
        return evalpoly(r, @inbounds J1_POLYS(T)[n]) * s
    else
        xinv = inv(x)
        iszero(xinv) && return zero(T)
        x2 = xinv * xinv
        if x < 120.0
            p1 = (one(T), 3/16, -99/512, 6597/8192, -4057965/524288, 1113686901/8388608, -951148335159/268435456, 581513783771781/4294967296)
            q1 = (3/8, -21/128, 1899/5120, -543483/229376, 8027901/262144, -30413055339/46137344, 9228545313147/436207616)
            p = evalpoly(x2, p1)
            q = evalpoly(x2, q1)
        else
            p2 = (one(T), 3/16, -99/512, 6597/8192)
            q2 = (3/8, -21/128, 1899/5120, -543483/229376)
            p = evalpoly(x2, p2)
            q = evalpoly(x2, q2)
        end

        a = SQ2OPI(T) * sqrt(xinv) * p
        xn = muladd(xinv, q, -3 * PIO4(T))

        # the following computes b = cos(x + xn) more accurately
        # see src/misc.jl
        b = cos_sum(x, xn)
        return a * b * s
    end
end
function besselj1(x::Float32)
    T = Float32
    s = sign(x)
    x = abs(x)

    if x <= 2.0f0
        z = x * x
        Z1 = 1.46819706421238932572f1
        p = (z - Z1) * x * evalpoly(z, JP32)
        return p * s
    else
        q = inv(x)
        iszero(q) && return zero(T)
        w = sqrt(q)
        p = w * evalpoly(q, MO132)
        w = q * q
        xn = q * evalpoly(w, PH132) - THPIO4(T)
        p = p * cos(xn + x)
        return p * s
    end
end

ext: jl


text:
""" Generic Julia interface for RDF graphs.
"""
module RDF
export Expression, Statement, Node, Edge, BaseURI, Prefix,
  Resource, ResourceURI, ResourceCURIE, Literal, Blank, Triple, Quad

using AutoHashEquals

# Data types
############

abstract type Expression end
abstract type Statement <: Expression end

@auto_hash_equals struct BaseURI <: Statement
  uri::String
end

@auto_hash_equals struct Prefix <: Statement
  name::String
  uri::String
end

abstract type Node <: Expression end
abstract type Resource <: Node end

@auto_hash_equals struct ResourceURI <: Resource
  uri::String
end

@auto_hash_equals struct ResourceCURIE <: Resource
  prefix::String
  name::String
end

@auto_hash_equals struct Literal <: Node
  value::Any
  language::String
  Literal(value::Any, language::String="") = new(value, language)
end

@auto_hash_equals struct Blank <: Node
  name::String
end

abstract type Edge <: Statement end

@auto_hash_equals struct Triple <: Edge
  subject::Union{Resource,Blank}
  predicate::Resource
  object::Node
end

@auto_hash_equals struct Quad <: Edge
  subject::Union{Resource,Blank}
  predicate::Resource
  object::Node
  graph::Resource
end

# Convenience constructors
Prefix(name::String) = Prefixes.prefix(name)
Resource(uri::String) = ResourceURI(uri)
Resource(prefix::String, name::String) = ResourceCURIE(prefix, name)

Edge(subj::Node, pred::Node, obj::Node) = Triple(subj, pred, obj)
Edge(subj::Node, pred::Node, obj::Node, graph::Node) = Quad(subj, pred, obj, graph)
Edge(subj::Node, pred::Node, obj::Node, graph::Nothing) = Triple(subj, pred, obj)

function Edge(subj::Node, pred::Node, obj::Node, graph::Union{<:Node,Nothing})
  isnothing(graph) ? Triple(subj, pred, obj) : Quad(subj, pred, obj, graph)
end

# Prefixes
##########

module Prefixes
export prefix, add_prefix!

using ..RDF: Prefix

const _prefixes = Dict{String,Prefix}()

function prefix(name::String)::Prefix
  _prefixes[name]
end

function add_prefix!(name::String, uri::String)
  if haskey(_prefixes, name)
    error("Prefix \"$name\" already defined")
  end
  _prefixes[name] = Prefix(name, uri)
end

add_prefix!("xsd", "http://www.w3.org/2001/XMLSchema#")
add_prefix!("rdf", "http://www.w3.org/1999/02/22-rdf-syntax-ns#")
add_prefix!("rdfs", "http://www.w3.org/2000/01/rdf-schema#")
add_prefix!("owl", "http://www.w3.org/2002/07/owl#")

add_prefix!("schema", "http://schema.org/")
add_prefix!("skos", "http://www.w3.org/2004/02/skos/core#")
add_prefix!("prov", "http://www.w3.org/ns/prov#")

end

end


text:
# BigBed Overlap
# ==============

struct OverlapIterator
    reader::Reader
    chromid::UInt32
    chromstart::UInt32
    chromend::UInt32
end

function Base.eltype(::Type{OverlapIterator})
    return Record
end

function Base.IteratorSize(::Type{OverlapIterator})
    return Base.SizeUnknown()
end

function GenomicFeatures.eachoverlap(reader::Reader, interval::GenomicFeatures.Interval)
    if haskey(reader.chroms, seqname(interval))
        id, _ = reader.chroms[seqname(interval)]
    else
        id = typemax(UInt32)
    end
    return OverlapIterator(reader, id, leftposition(interval) - 1, rightposition(interval))
end

mutable struct OverlapIteratorState
    state::BioCore.Ragel.State
    data::Vector{UInt8}
    done::Bool
    record::Record
    blocks::Vector{BBI.Block}
    current_block::Int
end

function Base.iterate(iter::OverlapIterator)
    data = Vector{UInt8}(undef, iter.reader.header.uncompress_buf_size)
    blocks = BBI.find_overlapping_blocks(iter.reader.index, iter.chromid, iter.chromstart, iter.chromend)
    if !isempty(blocks)
        seek(iter.reader.stream, blocks[1].offset)
    end
    state = OverlapIteratorState(
        BioCore.Ragel.State(
            data_machine.start_state,
            Libz.ZlibInflateInputStream(iter.reader.stream, reset_on_end=false)),
        data,
        isempty(blocks), Record(), blocks, isempty(blocks) ? 1 : 2)
    return iterate(iter, state)
end

function Base.iterate(iter::OverlapIterator, state::OverlapIteratorState)
    advance!(iter, state)
    if state.done
        return nothing
    end
    return copy(state.record), state
end

function advance!(iter::OverlapIterator, state::OverlapIteratorState)
    while true
        while state.current_block ≤ lastindex(state.blocks) && eof(state.state.stream)
            block = state.blocks[state.current_block]
            seek(iter.reader.stream, block.offset)
            size = BBI.uncompress!(state.data, read(iter.reader.stream, block.size))
            state.state = BioCore.Ragel.State(data_machine.start_state, BufferedStreams.BufferedInputStream(state.data[1:size]))
            state.current_block += 1
        end
        if state.done || (state.current_block > lastindex(state.blocks) && eof(state.state.stream))
            state.done = true
            return state
        end

        _read!(iter.reader, state.state, state.record)
        state.record.reader = iter.reader
        if overlaps(state.record, iter.chromid, iter.chromstart, iter.chromend)
            return state
        end
    end
end

function overlaps(record::Record, chromid::UInt32, chromstart::UInt32, chromend::UInt32)
    return record.chromid == chromid && !(record.chromend ≤ chromstart || record.chromstart ≥ chromend)
end


text:
using Images, TestImages
using Deconvolution
using FFTW
using Colors
using ZernikePolynomials
using FFTW
using NumberTheoreticTransforms

trim = 225:288 # need to limit size due to O(N^4) complexity
image_float = channelview(testimage("cameraman"))[trim.-110, trim.+40]
image_int = map(x -> x.:i, image_float) .|> Int64

blur_float = evaluateZernike(LinRange(-41,41,512), [12, 4, 0], [1.0, -1.0, 2.0], index=:OSA)
blur_float ./= (sum(blur_float))
blur_float = blur_float[trim, trim]
blur_int = blur_float .|> Normed{UInt8, 8} .|> x -> x.:i .|> Int64

blur_int = fftshift(blur_int)

@assert sum(blur_int) == 256 # to not change brightness

blurred_image_float = ifft(fft(image_float) .* fft(fftshift(blur_float))) |> real

(g, q, n) = (4, 274177, 64)
X = ntt(image_int, g, q)
@assert intt(X, g, q) == image_int
H = ntt( blur_int, g, q)
@assert intt(H, g, q) == blur_int
Y = mod.((X .* H), q)
y = intt(Y, g, q)
blurred_image_int = y / 2^16

using Statistics
@show maximum(abs.(blurred_image_float .- blurred_image_int))
@show mean(abs.(blurred_image_float .- blurred_image_int))
@show var(abs.(blurred_image_float .- blurred_image_int))

save("doc/src/ntt/original.jpg", Images.clamp01nan.(image_float))
save("doc/src/ntt/blurred_fft.jpg", Images.clamp01nan.(Gray.(Normed{UInt8, 8}.(blurred_image_float))))
save("doc/src/ntt/blurred_ntt.jpg", Images.clamp01nan.(Gray.(Normed{UInt8, 8}.(blurred_image_int))))


text:
using DiffEqBase, ParameterizedFunctions, DiffEqOperators, Random, LinearAlgebra
using Markdown

Random.seed!(100)

#ODE Example Problems
export prob_ode_linear, prob_ode_bigfloatlinear, prob_ode_2Dlinear,
       prob_ode_large2Dlinear, prob_ode_bigfloat2Dlinear, prob_ode_rigidbody,
       prob_ode_2Dlinear_notinplace, prob_ode_vanderpol, prob_ode_vanderpol_stiff,
       prob_ode_lorenz, prob_ode_rober, prob_ode_threebody, prob_ode_mm_linear, prob_ode_pleiades,
       prob_ode_brusselator_1d, prob_ode_brusselator_2d, prob_ode_orego,
       prob_ode_hires, prob_ode_pollution, prob_ode_filament,
       SolverDiffEq, filament_prob

include("ode_linear_prob.jl")
include("ode_simple_nonlinear_prob.jl")
include("brusselator_prob.jl")
include("pollution_prob.jl")
include("filament_prob.jl")


text:
# TODO: add a type for kinetic energy

struct Hamiltonian{M<:AbstractMetric, Tlogπ, T∂logπ∂θ}
    metric::M
    logπ::Tlogπ
    ∂logπ∂θ::T∂logπ∂θ
end

# Create a `Hamiltonian` with a new `M⁻¹`
(h::Hamiltonian)(M⁻¹) = Hamiltonian(h.metric(M⁻¹), h.logπ, h.∂logπ∂θ)

∂H∂θ(h::Hamiltonian, θ::AbstractVector) = -h.∂logπ∂θ(θ)

∂H∂r(h::Hamiltonian{<:UnitEuclideanMetric}, r::AbstractVector) = copy(r)
∂H∂r(h::Hamiltonian{<:DiagEuclideanMetric}, r::AbstractVector) = h.metric.M⁻¹ .* r
∂H∂r(h::Hamiltonian{<:DenseEuclideanMetric}, r::AbstractVector) = h.metric.M⁻¹ * r

function hamiltonian_energy(h::Hamiltonian, θ::AbstractVector, r::AbstractVector)
    K = kinetic_energy(h, r, θ)
    if isnan(K)
        K = Inf
        @warn "Kinetic energy is `NaN` and is set to `Inf`."
    end
    V = potential_energy(h, θ)
    if isnan(V)
        V = Inf
        @warn "Potential energy is `NaN` and is set to `Inf`."
    end
    return K + V
end

potential_energy(h::Hamiltonian, θ::AbstractVector) = -h.logπ(θ)

# Kinetic energy
# NOTE: the general form of K depends on both θ and r
kinetic_energy(h::Hamiltonian{<:UnitEuclideanMetric}, r, θ) = sum(abs2, r) / 2
function kinetic_energy(h::Hamiltonian{<:DiagEuclideanMetric}, r, θ)
    return sum(abs2(r[i]) * h.metric.M⁻¹[i] for i in 1:length(r)) / 2
end
function kinetic_energy(h::Hamiltonian{<:DenseEuclideanMetric}, r, θ)
    mul!(h.metric._temp, h.metric.M⁻¹, r)
    return dot(r, h.metric._temp) / 2
end

# Momentum sampler
function rand_momentum(rng::AbstractRNG, h::Hamiltonian{<:UnitEuclideanMetric})
    return randn(rng, h.metric.dim)
end
function rand_momentum(rng::AbstractRNG, h::Hamiltonian{<:DiagEuclideanMetric})
    r = randn(rng, h.metric.dim)
    r ./= h.metric.sqrtM⁻¹
    return r
end
function rand_momentum(rng::AbstractRNG, h::Hamiltonian{<:DenseEuclideanMetric})
    r = randn(rng, h.metric.dim)
    ldiv!(h.metric.cholM⁻¹, r)
    return r
end

rand_momentum(h::Hamiltonian) = rand_momentum(GLOBAL_RNG, h)


text:
# Base Thrift type system
struct TSTOP end
const TVOID     = Nothing
const TBOOL     = Bool
const TBYTE     = UInt8     # TBYTE is actually I8 (signed)
const TDOUBLE   = Float64
const TI16      = Int16
const TI32      = Int32
const TI64      = Int64
const TBINARY   = Vector{UInt8}
const TUTF8     = String
const TSTRING   = Union{TUTF8, TBINARY}
const TSTRUCT   = Any
const TMAP      = Dict
const TSET      = Set
const TLIST     = Array

abstract type TMsg end

struct _enum_TTypes
    STOP::Int32
    VOID::Int32
    BOOL::Int32
    BYTE::Int32
    DOUBLE::Int32
    I16::Int32
    I32::Int32
    I64::Int32
    STRING::Int32
    STRUCT::Int32
    MAP::Int32
    SET::Int32
    LIST::Int32
end

const TType = _enum_TTypes(0,    1,      2,      3,       4,                6,              8,              10,      11,       12,      13,    14,    15)
#                        0       1       2       3        4         5       6      7        8       9       10       11        12       13     14     15
const _TTypeNames   = ("STOP", "VOID", "BOOL", "BYTE", "DOUBLE",   "",    "I16",   "",     "I32",   "",     "I64", "STRING", "STRUCT", "MAP", "SET", "LIST")
const _TJTypes    =   (TSTOP,  TVOID,  TBOOL,  TBYTE,  TDOUBLE,  Nothing,  TI16, Nothing,  TI32,  Nothing,  TI64,  TSTRING,  TSTRUCT,  TMAP,  TSET,  TLIST)

#thrift_type_name(typ::Integer) = _TTypeNames[typ+1]

julia_type(typ::Integer) = _TJTypes[typ+1]
function julia_type(typ::Integer, narrow_typ)
    wide_typ = julia_type(typ)
    (narrow_typ <: wide_typ) && (return narrow_typ)
    error("Can not resolve type. $narrow_typ is not a subtype of $wide_typ")
end

thrift_type(::Type{TSTOP})                = Int32(0)
thrift_type(::Type{TVOID})                = Int32(1)
thrift_type(::Type{TBOOL})                = Int32(2)
thrift_type(::Type{TBYTE})                = Int32(3)
thrift_type(::Type{TDOUBLE})              = Int32(4)
thrift_type(::Type{TI16})                 = Int32(6)
thrift_type(::Type{TI32})                 = Int32(8)
thrift_type(::Type{TI64})                 = Int32(10)
thrift_type(::Type{TSTRING})              = Int32(11)
thrift_type(::Type{TUTF8})                = Int32(11)
thrift_type(::Type{TBINARY})              = Int32(11)
thrift_type(::Type{T}) where {T<:AbstractString} = Int32(11)
thrift_type(::Type{T}) where {T<:Any}            = Int32(12)
thrift_type(::Type{T}) where {T<:Dict}           = Int32(13)
thrift_type(::Type{T}) where {T<:Set}            = Int32(14)
thrift_type(::Type{T}) where {T<:Array}          = Int32(15)

const _container_type_ids   = (TType.STRUCT, TType.MAP, TType.SET, TType.LIST)
const _container_types      = (TSTRUCT, TMAP, TSET, TLIST)
const _plain_type_ids       = (TType.BOOL, TType.BYTE, TType.DOUBLE, TType.I16, TType.I32, TType.I64, TType.STRING)
const _plain_types          = (TBOOL, TBYTE, TDOUBLE, TI16, TI32, TI64, TBINARY, TUTF8)
iscontainer(typ::T)         where {T <: Integer}    = (Int32(typ) in _container_type_ids)
iscontainer(typ::Type{T})   where {T}               = iscontainer(thrift_type(typ))
isplain(typ::T)             where {T <: Integer}    = (Int32(typ) in _plain_type_ids)
isplain(typ::Type{T})       where {T}               = isplain(thrift_type(typ))

##
# base processor method
abstract type TProcessor end
process(tp::TProcessor) = nothing

##
# base transport types
abstract type TTransport end
abstract type TServerTransport <: TTransport end

##
# base server type
abstract type TServer end

##
# base protocol types
abstract type TProtocol end

##
# base protocol read and write methods
#instantiate(t::Type) = (t.abstract) ? error("can not instantiate abstract type $t") : ccall(:jl_new_struct_uninit, Any, (Any,Any...), t)
for _typ in _plain_types
    @eval begin
        write(p::TProtocol, val::$(_typ)) = 0
        read(p::TProtocol, ::Type{$(_typ)}) = nothing
        skip(p::TProtocol, ::Type{$(_typ)}) = read(p, $(_typ))
    end
end

writeMessageBegin(p::TProtocol, name::AbstractString, mtype::Int32, seqid::Integer)     = nothing
writeMessageEnd(p::TProtocol)                                                           = nothing
writeStructBegin(p::TProtocol, name::AbstractString)                                    = nothing
writeStructEnd(p::TProtocol)                                                            = nothing
writeFieldBegin(p::TProtocol, name::AbstractString, ttype::Int32, fid::Integer)         = nothing
writeFieldEnd(p::TProtocol)                                                             = nothing
writeFieldStop(p::TProtocol)                                                            = nothing
writeMapBegin(p::TProtocol, ktype::Int32, vtype::Int32, size::Integer)                  = nothing
writeMapEnd(p::TProtocol)                                                               = nothing
writeListBegin(p::TProtocol, etype::Int32, size::Integer)                               = nothing
writeListEnd(p::TProtocol)                                                              = nothing
writeSetBegin(p::TProtocol, etype::Int32, size::Integer)                                = nothing
writeSetEnd(p::TProtocol)                                                               = nothing
writeBool(p::TProtocol, val)                                                            = write(p, convert(TBOOL, val))
writeByte(p::TProtocol, val)                                                            = write(p, convert(TBYTE, val))
writeI16(p::TProtocol, val)                                                             = write(p, convert(TI16, val))
writeI32(p::TProtocol, val)                                                             = write(p, convert(TI32, val))
writeI64(p::TProtocol, val)                                                             = write(p, convert(TI64, val))
writeDouble(p::TProtocol, val)                                                          = write(p, convert(TDOUBLE, val))
writeString(p::TProtocol, val)                                                          = write(p, String(val))
writeBinary(p::TProtocol, val)                                                          = write(p, convert(TBINARY, val))

readMessageBegin(p::TProtocol)     = nothing
readMessageEnd(p::TProtocol)       = nothing
readStructBegin(p::TProtocol)      = nothing
readStructEnd(p::TProtocol)        = nothing
readFieldBegin(p::TProtocol)       = nothing
readFieldEnd(p::TProtocol)         = nothing
readMapBegin(p::TProtocol)         = nothing
readMapEnd(p::TProtocol)           = nothing
readListBegin(p::TProtocol)        = nothing
readListEnd(p::TProtocol)          = nothing
readSetBegin(p::TProtocol)         = nothing
readSetEnd(p::TProtocol)           = nothing
readBool(p::TProtocol)             = read(p, TBOOL)
readByte(p::TProtocol)             = read(p, TBYTE)
readI16(p::TProtocol)              = read(p, TI16)
readI32(p::TProtocol)              = read(p, TI32)
readI64(p::TProtocol)              = read(p, TI64)
readDouble(p::TProtocol)           = read(p, TDOUBLE)
readString(p::TProtocol)           = read(p, TUTF8)
readBinary(p::TProtocol)           = read(p, TBINARY)

skip(p::TProtocol, ::Type{T}) where {T<:TSTRUCT} = skip_container(p, T)
function skip_container(p::TProtocol, ::Type{T}) where T<:TSTRUCT
    @debug("skip TSTRUCT")
    name = readStructBegin(p)
    while true
        (name, ttype, id) = readFieldBegin(p)
        (ttype == TType.STOP) && break
        if iscontainer(ttyp)
            skip_container(p, julia_type(ttype))
        else
            skip(p, julia_type(ttype))
        end
        readFieldEnd(p)
    end
    readStructEnd(p)
    return
end

read(p::TProtocol, ::Type{T}) where {T<:TSTRUCT} = read(p, T())
read_container(p::TProtocol, ::Type{T}) where {T<:TSTRUCT} = read_container(p, T())
read(p::TProtocol, val::T) where {T<:TSTRUCT} = read_container(p, val)
function read_container(p::TProtocol, val::T) where T<:TSTRUCT
    @debug("read TSTRUCT", T)
    readStructBegin(p)

    m = meta(T)
    @debug("struct meta", meta=m)
    clear(val)
    while true
        (name, ttyp, id) = readFieldBegin(p)
        (ttyp == TType.STOP) && break

        attribs = m.numdict[Int(id)]
        jtyp = julia_type(attribs)
        fldname = attribs.fld
        if iscontainer(ttyp)
            if isdefined(val, fldname)
                @debug("reading into already defined container field", jtyp, fldname)
                read_container(p, getfield(val, fldname))
            else
                @debug("setting into container field", jtyp, fldname)
                setproperty!(val, fldname, read_container(p, jtyp))
            end
        else
            @debug("setting field", jtyp, fldname)
            setproperty!(val, fldname, read(p, jtyp))
        end
        readFieldEnd(p)
    end
    readStructEnd(p)

    # populate remaining with any defaults
    setdefaultproperties!(val)
end

function setdefaultproperties!(val::T) where T<:TSTRUCT
    m = meta(T)
    for attrib in m.ordered
        fldname = attrib.fld
        if !hasproperty(val, fldname) && !isempty(attrib.default)
            default = attrib.default[1]
            setproperty!(val, fldname, deepcopy(default))
        end
    end
    val
end

write(p::TProtocol, val::T) where {T<:TSTRUCT} = write_container(p, val)
function write_container(p::TProtocol, val::T) where T<:TSTRUCT
    m = meta(T)
    @debug("write TSTRUCT", T, meta=m)
    writeStructBegin(p, string(T))

    for attrib in m.ordered
        if !hasproperty(val, attrib.fld)
            m.symdict[attrib.fld].required && error("required field $(attrib.fld) not populated")
            continue
        end
        writeFieldBegin(p, string(attrib.fld), attrib.ttyp, attrib.fldnum)
        fld = getproperty(val, attrib.fld)
        if (attrib.ttyp == TType.STRING) && isa(fld, Vector{UInt8})
            write(p, fld, true)
        elseif attrib.ttyp == TType.BOOL
            writeBool(p, fld)
        else
            write(p, fld)
        end
        writeFieldEnd(p)
    end
    writeFieldStop(p)
    writeStructEnd(p)
    nothing
end

skip(p::TProtocol, ::Type{T}) where {T<:TMAP} = skip_container(p, T)
function skip_container(p::TProtocol, ::Type{T}) where T<:TMAP
    @debug("skip TMAP", T)
    (ktype, vtype, size) = readMapBegin(p)
    if size > 0
        jktype = julia_type(ktype)
        jvtype = julia_type(vtype)
        for i in 1:size
            skip(p, jktype)
            skip(p, jvtype)
        end
    end
    readMapEnd(p)
    return
end

read(p::TProtocol, ::Type{T}) where {T<:TMAP} = read(p, T())
read(p::TProtocol, val::T) where {T<:TMAP} = read_container(p, val)
function read_container(p::TProtocol, val::T) where T<:TMAP
    @debug("read TMAP", T)
    (ktype, vtype, size) = readMapBegin(p)
    if size > 0
        # types are valid only when size is non zero
        (_ktype, _vtype) = fieldtypes(eltype(val))
        jktype = julia_type(ktype, _ktype)
        jvtype = julia_type(vtype, _vtype)
        for i in 1:size
            k = read(p, jktype)
            v = read(p, jvtype)
            val[k] = v
        end
    end
    readMapEnd(p)
    val
end

write(p::TProtocol, val::TMAP) = write_container(p, val)
function write_container(p::TProtocol, val::TMAP)
    @debug("write TMAP", valtype=typeof(val), size=length(val))
    (ktype,vtype) = fieldtypes(eltype(val))
    writeMapBegin(p, thrift_type(ktype), thrift_type(vtype), length(val))
    for (k,v) in val
        @debug("write TMAP key")
        if ktype === Vector{UInt8}
            write(p, k, true)
        else
            write(p, k)
        end
        @debug("write TMAP value")
        if vtype === Vector{UInt8}
            write(p, v, true)
        else
            write(p, v)
        end
    end
    writeMapEnd(p)
end

skip(p::TProtocol, ::Type{T}) where {T<:TSET} = skip_container(p, T)
function skip_container(p::TProtocol, ::Type{T}) where T<:TSET
    @debug("skip TSET", T)
    (etype, size) = readSetBegin(p)
    if size > 0
        jetype = julia_type(etype)
        for i in 1:size
            skip(p, jetype)
        end
    end
    readSetEnd(p)
end

read(p::TProtocol, ::Type{T}) where {T<:TSET} = read(p, T())
read(p::TProtocol, val::T) where {T<:TSET} = read_container(p, val)
function read_container(p::TProtocol, val::T) where T<:TSET
    @debug("read TSET", T)
    (etype, size) = readSetBegin(p)
    if size > 0
        jetype = julia_type(etype, eltype(val))
        for i in 1:size
            push!(val, read(p, jetype))
        end
    end
    readSetEnd(p)
    val
end

write(p::TProtocol, val::TSET) = write_container(p, val)
function write_container(p::TProtocol, val::TSET)
    @debug("write TSET", valtype=typeof(val), size=length(val))
    jetype = eltype(val)
    tetype = thrift_type(jetype)
    writeSetBegin(p, tetype, length(val))

    if iscontainer(tetype)
        for v in val
            write(p, jetype, v)
        end
    else
        for v in val
            if jetype === Vector{UInt8}
                write(p, v, true)
            else
                write(p, v)
            end
        end
    end
    writeSetEnd(p)
end

skip(p::TProtocol, ::Type{T}) where {T<:TLIST} = skip_container(p, T)
function skip_container(p::TProtocol, ::Type{T}) where T<:TLIST
    @debug("skip TLIST", T)
    (etype, size) = readListBegin(p)
    if size > 0
        jetype = julia_type(etype)
        for i in 1:size
            skip(p, jetype)
        end
    end
    readListEnd(p)
end

read(p::TProtocol, ::Type{T}) where {T<:TLIST} = read(p, T())
read(p::TProtocol, val::T) where {T<:TLIST} = read_container(p, val)
function read_container(p::TProtocol, val::T) where T<:TLIST
    @debug("read TLIST", T)
    (etype, size) = readListBegin(p)
    if size > 0
        jetype = julia_type(etype, eltype(val))
        for i in 1:size
            push!(val, read(p, jetype))
        end
    end
    readListEnd(p)
    val
end

write(p::TProtocol, val::TLIST) = write_container(p, val)
function write_container(p::TProtocol, val::TLIST)
    @debug("write TLIST", valtype=typeof(val), size=length(val))
    etype = eltype(val)
    writeListBegin(p, thrift_type(etype), length(val))
    # TODO: need meta to convert type correctly
    for v in val
        if etype === Vector{UInt8}
            write(p, v, true)
        else
            write(p, v)
        end
    end
    writeListEnd(p)
    nothing
end

##
# Thrift Structure Metadata

mutable struct ThriftMetaAttribs
    fldnum::Int                     # the field number in the structure
    fld::Symbol
    ttyp::Int32                     # thrift type
    jtype::Type                     # Julia type
    required::Bool                  # required or optional
    default::Vector                 # the default value, empty array if none is specified, first element is used if something is specified
    elmeta::Vector                  # the ThriftMeta of a struct or element/key-value types if this is a list, set or map
end

mutable struct ThriftMeta
    jtype::Type
    symdict::Dict{Symbol,ThriftMetaAttribs}
    numdict::Dict{Int,ThriftMetaAttribs}
    ordered::Vector{ThriftMetaAttribs}

    ThriftMeta(jtype::Type, ordered::Vector{ThriftMetaAttribs}) = _setmeta(new(), jtype, ordered)
end

function _setmeta(meta::ThriftMeta, jtype::Type, ordered::Vector{ThriftMetaAttribs})
    symdict = Dict{Symbol,ThriftMetaAttribs}()
    numdict = Dict{Int,ThriftMetaAttribs}()
    for attrib in ordered
        symdict[attrib.fld] = numdict[attrib.fldnum] = attrib
    end
    meta.jtype = jtype
    meta.symdict = symdict
    meta.numdict = numdict
    meta.ordered = ordered
    meta
end

julia_type(fattr::ThriftMetaAttribs) = fattr.jtype

meta(typ::Type) = meta(typ, Symbol[], Type[], Symbol[], Int[], Dict{Symbol,Any}())
function meta(typ::Type, names::Vector{Symbol}, types::Vector{Type}, optional::Vector{Symbol}, numbers::Vector{Int}, defaults::Dict{Symbol,Any})
    m = ThriftMeta(typ, ThriftMetaAttribs[])

    attribs = ThriftMetaAttribs[]
    for fldidx in 1:length(names)
        fldtyp = types[fldidx]
        fldttyp = thrift_type(fldtyp)
        fldname = names[fldidx]
        fldnum = Int(isempty(numbers) ? fldidx : numbers[fldidx])
        fldrequired = !(fldname in optional)

        elmeta = ThriftMeta[]
        if fldttyp == TType.STRUCT
            push!(elmeta, meta(fldtyp))
        elseif fldttyp == TType.LIST
            push!(elmeta, meta(fldtyp.parameters[1]))
        elseif fldttyp == TType.SET
            push!(elmeta, meta(fldtyp.parameters[1]))
        elseif fldttyp == TType.MAP
            push!(elmeta, meta(fldtyp.parameters[1]))   # key
            push!(elmeta, meta(fldtyp.parameters[2]))   # value
        end

        default = haskey(defaults, fldname) ? Any[defaults[fldname]] : []

        push!(attribs, ThriftMetaAttribs(fldnum, fldname, fldttyp, fldtyp, fldrequired, default, elmeta))
    end
    _setmeta(m, typ, attribs)
    m
end

function show(io::IO, m::ThriftMeta)
    println(io, "ThriftMeta for $(m.jtype)")
    println(io, m.ordered)
end

@deprecate fillunset(obj) clear(obj)
@deprecate isfilled(obj, fld) hasproperty(obj, fld)
function isfilled(obj)
    fill = keys(obj.values)
    flds = obj.meta.ordered
    for fld in flds
        if fld.required
            (fld.fld in fill) || (return false)
            !isempty(fld.elmeta) && !isfilled(getproperty(obj, fld.fld)) && (return false)
        end
    end
    true
end
function isfilled(vec::Vector{T}) where T<:TMsg
    isempty(vec) || all(isfilled.(vec))
end

##
# utility methods
function copy!(to::T, from::T) where T<:TMsg
    clear(to)
    for name in propertynames(from)
        if hasproperty(from, name)
            setproperty!(to, name, getproperty(from, name))
        end
    end
    nothing
end

isinitialized(obj) = isfilled(obj)

@deprecate set_field!(obj, fld, val) setproperty!(obj, fld, val)
@deprecate get_field(obj, fld) getproperty(obj, fld)
@deprecate has_field(obj, fld) hasproperty(obj, fld)

propertynames(obj::TMsg) = collect(keys(obj.meta.symdict))
hasproperty(obj::TMsg, fld::Symbol) = haskey(obj.values, fld)
function setproperty!(obj::TMsg, fld::Symbol, val)
    symdict = obj.meta.symdict
    if fld in keys(symdict)
        fldtype = symdict[fld].jtype
        obj.values[fld] = isa(val, fldtype) ? val : convert(fldtype, val)
    else
        setfield!(obj, fld, val)
    end
end

function clear(obj::TMsg)
    empty!(obj.values)
    nothing
end

function clear(obj::TMsg, fld::Symbol)
    delete!(obj.values, fld)
    nothing
end

function thriftbuild(::Type{T}, nv::Dict{Symbol}=Dict{Symbol,Any}()) where T
    Base.depwarn("thriftbuild is deprecated; use constructor instead", :thriftbuild)
    T(; nv...)
end

function enumstr(enumname, t::Int32)
    T = typeof(enumname)
    for name in fieldnames(T)
        (getfield(enumname, name) == t) && (return string(name))
    end
    error("Invalid enum value $t for $T)")
end

##
# Exception types
"""
Exception that must be thrown by application code on any unexpected error.
"""
mutable struct TException <: TMsg
    meta::ThriftMeta
    values::Dict{Symbol,Any}

    function TException(; kwargs...)
        obj = new(__meta__TException, Dict{Symbol,Any}())
        values = obj.values
        symdict = obj.meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtype
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        Thrift.setdefaultproperties!(obj)
        obj
    end
end

const __meta__TException = meta(TException,
    Symbol[:message],
    Type[String],
    Symbol[],
    Int[],
    Dict{Symbol,Any}()
)

function Base.getproperty(obj::TException, name::Symbol)
    if name === :message
        return (obj.values[name])::String
    else
        getfield(obj, name)
    end
end

meta(::Type{TException}) = __meta__TException

struct _enum_TApplicationExceptionTypes
    UNKNOWN::Int32
    UNKNOWN_METHOD::Int32
    INVALID_MESSAGE_TYPE::Int32
    WRONG_METHOD_NAME::Int32
    BAD_SEQUENCE_ID::Int32
    MISSING_RESULT::Int32
    INTERNAL_ERROR::Int32
    PROTOCOL_ERROR::Int32
    INVALID_TRANSFORM::Int32
    INVALID_PROTOCOL::Int32
    UNSUPPORTED_CLIENT_TYPE::Int32
end

const ApplicationExceptionType = _enum_TApplicationExceptionTypes(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4), Int32(5), Int32(6), Int32(7), Int32(8), Int32(9), Int32(10))
const _appex_msgs = [
    "Default (unknown) TApplicationException",
    "Unknown method",
    "Invalid message type",
    "Wrong method name",
    "Bad sequence ID",
    "Missing result",
    "Internal error",
    "Protocol error",
    "Invalid transform",
    "Invalid protocol",
    "Unsupported client type"
]

"""
Exception thrown by Thrift.jl internals on any unexpected error.
"""
mutable struct TApplicationException <: TMsg
    meta::ThriftMeta
    values::Dict{Symbol,Any}

    function TApplicationException(; kwargs...)
        obj = new(__meta__TApplicationException, Dict{Symbol,Any}())
        values = obj.values
        symdict = obj.meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtype
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        Thrift.setdefaultproperties!(obj)
        obj
    end
end

const __meta__TApplicationException = meta(TApplicationException,
    Symbol[:typ, :message],
    Type[Int32, TUTF8],
    Symbol[],
    Int[],
    Dict{Symbol,Any}(
        :typ => ApplicationExceptionType.UNKNOWN,
        :message => ""
    )
)

function Base.getproperty(obj::TApplicationException, name::Symbol)
    if name === :typ
        return (obj.values[name])::Int32
    elseif name === :message
        return (obj.values[name])::TUTF8
    else
        getfield(obj, name)
    end
end

meta(::Type{TApplicationException}) = __meta__TApplicationException

##
# Message types
struct _enum_TMessageType
    CALL::Int32
    REPLY::Int32
    EXCEPTION::Int32
    ONEWAY::Int32
end

const MessageType = _enum_TMessageType(Int32(1), Int32(2), Int32(3), Int32(4))


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""EnvFromSource represents the source of a set of ConfigMaps

    IoK8sApiCoreV1EnvFromSource(;
        configMapRef=nothing,
        prefix=nothing,
        secretRef=nothing,
    )

    - configMapRef::IoK8sApiCoreV1ConfigMapEnvSource : The ConfigMap to select from
    - prefix::String : An optional identifier to prepend to each key in the ConfigMap. Must be a C_IDENTIFIER.
    - secretRef::IoK8sApiCoreV1SecretEnvSource : The Secret to select from
"""
mutable struct IoK8sApiCoreV1EnvFromSource <: SwaggerModel
    configMapRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ConfigMapEnvSource } # spec name: configMapRef
    prefix::Any # spec type: Union{ Nothing, String } # spec name: prefix
    secretRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1SecretEnvSource } # spec name: secretRef

    function IoK8sApiCoreV1EnvFromSource(;configMapRef=nothing, prefix=nothing, secretRef=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1EnvFromSource, Symbol("configMapRef"), configMapRef)
        setfield!(o, Symbol("configMapRef"), configMapRef)
        validate_property(IoK8sApiCoreV1EnvFromSource, Symbol("prefix"), prefix)
        setfield!(o, Symbol("prefix"), prefix)
        validate_property(IoK8sApiCoreV1EnvFromSource, Symbol("secretRef"), secretRef)
        setfield!(o, Symbol("secretRef"), secretRef)
        o
    end
end # type IoK8sApiCoreV1EnvFromSource

const _property_map_IoK8sApiCoreV1EnvFromSource = Dict{Symbol,Symbol}(Symbol("configMapRef")=>Symbol("configMapRef"), Symbol("prefix")=>Symbol("prefix"), Symbol("secretRef")=>Symbol("secretRef"))
const _property_types_IoK8sApiCoreV1EnvFromSource = Dict{Symbol,String}(Symbol("configMapRef")=>"IoK8sApiCoreV1ConfigMapEnvSource", Symbol("prefix")=>"String", Symbol("secretRef")=>"IoK8sApiCoreV1SecretEnvSource")
Base.propertynames(::Type{ IoK8sApiCoreV1EnvFromSource }) = collect(keys(_property_map_IoK8sApiCoreV1EnvFromSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1EnvFromSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1EnvFromSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1EnvFromSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1EnvFromSource[property_name]

function check_required(o::IoK8sApiCoreV1EnvFromSource)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1EnvFromSource }, name::Symbol, val)
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    SamplingMethod

A method for sampling from spatial objects.
"""
abstract type SamplingMethod end

"""
    sample(object, method)

Sample elements from `object` with `method`.
"""
function sample end

#------------------
# IMPLEMENTATIONS
#------------------
include("sampling/uniform.jl")
include("sampling/ball.jl")
include("sampling/weighted.jl")


text:
# run without startup file
# julia --startup-file=no create_sysimage.jl

using Pkg
using PackageCompiler: create_sysimage
Pkg.activate(".")
deps = [pair.second for pair in Pkg.dependencies()]
direct_deps = filter(p -> p.is_direct_dep, deps)
pkg_name_version = [(x.name, x.version) for x in direct_deps]
pkg_list = [Symbol(x.name) for x in direct_deps]
create_sysimage(
    pkg_list; sysimage_path="cnode.so", precompile_execution_file="src/ControlNeuralODE.jl"
)

# how to use sysimage to execute examples & activate current project:
# julia --sysimage cnode.so --project=.


text:
struct ∂Tag end

@inline function derivative(f::F, x::R) where {F,R<:Real}
    T = typeof(ForwardDiff.Tag(∂Tag(), R))
    out = f(ForwardDiff.Dual{T,R,1}(x, ForwardDiff.Partials((one(R),))))
    return ForwardDiff.extract_derivative(T, out)
end

@inline function gradient2(f::F, x1::R,x2::R) where {F,R<:Real}
    T = typeof(ForwardDiff.Tag(∂Tag(), R))
    _1 = one(R)
    _0 = zero(R)
    dual1 = ForwardDiff.Dual{T,R,2}(x1, ForwardDiff.Partials((_1,_0)))
    dual2 = ForwardDiff.Dual{T,R,2}(x2, ForwardDiff.Partials((_0,_1)))
    out = f(dual1,dual2)
    ∂out = ForwardDiff.partials(out)
    return SVector(∂out.values)
end

function gradient2(f::F,x1::R1,x2::R2) where{F,R1<:Real,R2<:Real}
    y1,y2 = promote(x1,x2)
    return gradient2(f,y1,y2)
end

"""
    f∂f(f,x)

returns f and ∂f/∂x evaluated in `x`, using `ForwardDiff.jl`, `DiffResults.jl` and `StaticArrays.jl` to calculate everything in one pass.
"""
@inline function f∂f(f::F, x::R) where {F,R<:Real}
    T = typeof(ForwardDiff.Tag(∂Tag(), R))
    out = f(ForwardDiff.Dual{T,R,1}(x, ForwardDiff.Partials((one(R),))))
    return ForwardDiff.value(out),  ForwardDiff.extract_derivative(T, out)
end

"""
    f∂f∂2f(f,x)

returns f,∂f/∂x,and ∂²f/∂²x and evaluated in `x`, using `ForwardDiff.jl`, `DiffResults.jl` and `StaticArrays.jl` to calculate everything in one pass.
"""
@inline function f∂f∂2f(f::F,x::R) where {F,R<:Real}
    T = typeof(ForwardDiff.Tag(∂Tag(), R))
    out = ForwardDiff.Dual{T,R,1}(x, ForwardDiff.Partials((one(R),)))
    _f,_df = f∂f(f,out)
    fx = ForwardDiff.value(_f)
    dfx = ForwardDiff.partials(_f).values[1]
    d2fx = ForwardDiff.partials(_df).values[1]
    return (fx,dfx,d2fx)
end

"""
    fgradf2(f,x1,x2)

returns f and ∇f(x),using `ForwardDiff.jl`
"""
function fgradf2(f::F,x1::R1,x2::R2) where{F,R1<:Real,R2<:Real}
    y1,y2 = promote(x1,x2)
    return fgradf2(f,y1,y2)
end

@inline function fgradf2(f::F,x1::R,x2::R) where{F,R<:Real}
    T = typeof(ForwardDiff.Tag(∂Tag(), R))
    _1 = one(R)
    _0 = zero(R)
    dual1 = ForwardDiff.Dual{T,R,2}(x1, ForwardDiff.Partials((_1,_0)))
    dual2 = ForwardDiff.Dual{T,R,2}(x2, ForwardDiff.Partials((_0,_1)))
    out = f(dual1,dual2)
    ∂out = ForwardDiff.partials(out)
    return ForwardDiff.value(out),SVector(∂out.values)
end

#Manual implementation of an hyperdual.
@inline function ∂2(f::F,x1::R,x2::R) where {F,R<:Real}
    T = typeof(ForwardDiff.Tag(∂Tag(), R))
    _1 = one(R)
    _0 = zero(R)
    dual1 = ForwardDiff.Dual{T,R,2}(x1, ForwardDiff.Partials((_1,_0)))
    dual2 = ForwardDiff.Dual{T,R,2}(x2, ForwardDiff.Partials((_0,_1)))
    _f,_df = fgradf2(f,dual1,dual2)
    fx = ForwardDiff.value(_f)
    df1,df2 = _df[1],_df[2]
    df = SVector(df1.value , df2.value)
    d2fdx2, d2fdxdy = df1.partials.values
          _, d2fdy2 = df2.partials.values
    d2f = SMatrix{2}(d2fdx2,d2fdxdy,d2fdxdy,d2fdy2)
    return (fx,df,d2f)
end

function ∂2(f::F,x1::R1,x2::R2) where{F,R1<:Real,R2<:Real}
    y1,y2 = promote(x1,x2)
    return ∂2(f,y1,y2)
end


text:
master = load("master-bench.jld")["bench"]
optim = load("optimize-bench.jld")["bench"]

expmv_t = [ratio(median(optim[i]["expmv"]), median(master[i]["expmv"])).time for i in 1:6]

expmv_m = [ratio(median(optim[i]["expmv"]), median(master[i]["expmv"])).memory for i in 1:6]

time = reduce(hcat,[[optim[i]["d"],median(optim[i]["expm"]).time,median(optim[i]["expokit"]).time,median(optim[i]["expmv"]).time] for i in 1:6])'

memory = reduce(hcat,[[optim[i]["d"],median(optim[i]["expm"]).memory,median(optim[i]["expokit"]).memory,median(optim[i]["expmv"]).memory] for i in 1:6])'


text:
using ObjFileBase
import AbstractTrees: children
import Base: start, next, done
import Base: print, show

immutable CU
    offset::UInt
    ats::AbbrevTableSet
    header::DWARFCUHeader
end

immutable LightDIERef
    io::IO
    cu::CU
    offset::UInt
end

# Small helper to seek to the start of a die ref and get the abbrev entry
function seekentry(die::LightDIERef)
    seek(die.io, die.offset)
    abbrev = UInt(read(die.io, ULEB128))
    abbrev == 0 && return zero_entry
    die.cu.ats.entries[abbrev]
end

tag(die::LightDIERef) = seekentry(die).tag

# Iterate over a DIEs attributes
immutable AttributeIterator
    die::LightDIERef
    ate::AbbrevTableEntry
end
function AttributeIterator(die::LightDIERef)
    new(die,seekentry(die))
end
function start(it::AttributeIterator)
    seekentry(it.die)
    (position(it.die.io), 1)
end
function next(it::AttributeIterator, state)
    endianness = :NativeEndian
    seek(it.die.io,state[1])
    idx = state[2]
    form, value = readorskip(it.die, it.ate.attributes[idx].form, endianness, Val{:read}())
    ret = Attribute(AttributeSpecification(it.ate.attributes[idx].name, form), value)
    (ret, (position(it.die.io), state[2]+1))
end
done(it::AttributeIterator, state) = state[2] > length(it.ate.attributes)
attributes(die::LightDIERef) = AttributeIterator(die)

# Show on a LightDIERef shows the contents of the DIE
function show(io::IO, die::LightDIERef)
    ate = seekentry(die)
    print_with_color(tag_color, io, tag_name(ate.tag))
    println(io)
    for attr in AttributeIterator(die, ate)
        print(io, "  "); show(io, attr); println(io)
    end
end
function Base.showcompact(io::IO, die::LightDIERef)
    print(io,"(cu + $(die.offset - die.cu.offset) => {$(hex(die.offset,2sizeof(die.offset)))})")
end

# Iterate over a DIEs children (not including attributes)
immutable ChildIterator
     die::LightDIERef
end
children(die::LightDIERef) = ChildIterator(die)

# State is (next_offset, next_abbrev)
function start(it::ChildIterator)
    ate = seekentry(it.die)
    ate.has_children == 0 && return (UInt(0),UInt(0))
    skip_attributes(it.die, ate)
    next_offset = position(it.die.io)
    next_abbrev = UInt(read(it.die.io, ULEB128))
    (next_offset, next_abbrev)
end
function next(it::ChildIterator, state)
    io = it.die.io
    seek(io, state[1])
    ate = it.die.cu.ats.entries[state[2]]
    ret = LightDIERef(io, it.die.cu, state[1])
    # Could just determine the size from state[2] and have state[1] be after the
    # abbrev, but so be it.
    read(io, ULEB128)
    skip_attributes(ret, ate)
    ate.has_children != 0 && skip_children(io, it.die)
    next_offset = position(io)
    next_abbrev = UInt(read(io, ULEB128))
    (ret, (next_offset, next_abbrev))
end
done(it::ChildIterator, state) = state[2] == 0

function skip_children(io::IO, cu)
    level = 1
    while level >= 1
        abbrev = UInt(read(io, ULEB128))
        if abbrev == 0
            level -= 1
            continue
        end
        ate = cu.cu.ats.entries[abbrev]
        skip_attributes(cu, ate)
        ate.has_children != 0 && (level += 1)
    end
end

# Iterate over units in a section
immutable UnitIterator{T}
    unit_section::SectionRef
    debug_abbrev::SectionRef
end
start(x::UnitIterator) = sectionoffset(x.unit_section)
function next{T}(x::UnitIterator{T}, offset)
    io = ObjFileBase.handle(x.unit_section).io
    seek(io, offset)
    endianness = :NativeEndian
    cuheader = read(io, T, endianness)
    seek(x.debug_abbrev, cuheader.debug_abbrev_offset)
    ats = read(io, AbbrevTableSet, endianness)
    die_offset = offset + sizeof(cuheader) +
      (isa(cuheader.unit_length, UInt64) ? sizeof(UInt32) : 0)
    ret = LightDIERef(io, CU(offset, ats, cuheader), die_offset)
    next_offset = die_offset + cuheader.unit_length - sizeof(cuheader) + sizeof(cuheader.unit_length)
    (ret, next_offset)
end
done(x::UnitIterator, off) = off >= sectionoffset(x.unit_section) + sectionsize(x.unit_section)
Base.iteratorsize{T<:UnitIterator}(::Type{T}) = Base.SizeUnknown()

function realize(ref::LightDIERef)
    seek(ref.io, ref.offset)
    abbrev = UInt(read(ref.io, ULEB128))
    read(ref.io, ref.cu.header, ref.cu.ats.entries[abbrev], DIE, :NativeEndian)
end

function readnthattr(ref, ate, endianness, idx)
    for i = 1:(idx-1)
        readorskip(ref, ate.attributes[i].form, :NativeEndian, Val{:skip}())
    end
    # Form could be indirect, so we need the return value
    readorskip(ref, ate.attributes[idx].form, endianness, Val{:read}())
end

function extract_attribute(ref::LightDIERef, name, endianness = :NativeEndian)
    io = ref.io
    ate = seekentry(ref)
    idx = findfirst(x->x.name==name, ate.attributes)
    if idx == 0
        idx = findfirst(x->x.name==DWARF.DW_AT_abstract_origin, ate.attributes)
        (idx == 0) && return Nullable{Attribute}()
        form, value = readnthattr(ref, ate, endianness, idx)
        return extract_attribute(value, name, endianness)
    end
    form, value = readnthattr(ref, ate, endianness, idx)
    Nullable{Attribute}(Attribute(AttributeSpecification(ate.attributes[idx].name, form), value))
end

function skip_attributes(ref::LightDIERef, ate)
    for at in ate.attributes
        readorskip(ref, at.form, :NativeEndian, Val{:skip}())
    end
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 0900dd0b-483a-404d-bf84-43c075230b73
begin
	using Pkg
	Pkg.activate("../../")
	using Revise
end

# ╔═╡ aca53457-93bb-4e6c-a3dc-68c000e8db78
Pkg.add("ImageIO")

# ╔═╡ 3d47cb3e-43b3-11ec-1b61-992ecbb3d16f
using Colors, ImageShow, TestImages, Noise, PlutoUI, Plots, FFTW

# ╔═╡ 5cdc016a-b0aa-42b4-841c-657c8bf72ca3
Pkg.add("FFTW")

# ╔═╡ e1f76c58-fa78-4a20-b30f-5ef087874899
md" ## Detectors

More bits $\rightarrow$ finer resolution of the intensity.

A lot of detectors have only 8 Bit (1 Byte)
* 8 Bit means: each pixel can measure $2^8$ different values

Some sensor have 12 Bit pixels
* 12 Bit means: each pixel can measure $2^{12}$ different values

Image sensors distributes those $N$ Bits in a certain interval equally spaced: $[0, I_\text{max}]$

## Examples
2 Bits per pixel (4 different values), the sensor can only measure: $[0, \frac13 I_\text{max}, \frac23 I_\text{max}, I_\text{max}]$

## Float64
Floating works differently. They store numbers not equally spaced, but more like exponentially spaced.
Of course, you can count how many possible values are in a certain interval $[0, 1.23]$, however that's not too easy to calculate.
"

# ╔═╡ c5e12406-97b2-450d-9fe7-b2da82452f33
2^8

# ╔═╡ 73f690da-3ee9-4a27-8e97-b062a615266e
2^12

# ╔═╡ 1f73e297-6a82-4d33-8862-5a166da4dc2f
# Float64 -> 64bits -> 2^64
img = Float64.(testimage("mandril_gray"));

# ╔═╡ fcbb40e5-87be-49f8-968b-e8bdd962ee53
# Image Sensor with only 2 Bits -> 4 values
Gray.(quantization(img, 2^2))

# ╔═╡ a564a1d9-0c5f-4391-bae4-04ffc9abc861
# Image Sensor with only 4 Bits -> 16 values
Gray.(quantization(img, 2^4))

# ╔═╡ c4087dcf-ea75-4b57-b4cc-73934f509af8
# Image Sensor with only 4 Bits -> 256 values
Gray.(quantization(img, 2^8))

# ╔═╡ 4c4199cf-a62f-4b82-ac04-0b1cb9b2543f

# ╔═╡ fc992ffe-1f4e-45e8-b2f2-7e4ec53d2188

# ╔═╡ 26eb3fff-ece9-40c9-bd9d-32cba04746a1
md"# Homework"

# ╔═╡ 09c7396d-57d4-4e26-a0bf-7b5730836135
"""
	to_equal_intensity_tuple(value::T) where T

Returns for a given scalar input a tuple of equal intensity.
There is no unique solution, we only require that the sum is always 1.
For example, below would be a valid solution.

## Examples
```julia-repl
julia> to_equal_intensity_tuple(1.0)
(0.04465819873852073, 0.33333333333333326, 0.6220084679281465)

julia> to_equal_intensity_tuple(0.0)
(0.044658198738520505, 0.3333333333333333, 0.6220084679281462)

julia> to_equal_intensity_tuple(0.3f0)
(0.26402956f0, 0.6503522f0, 0.08561832f0)
```
"""
function to_equal_intensity_tuple(value::T) where T
	# with the function signature we can access T (which is the type of `value`)
	# 0 has type Int
	# 1.0 has type Float64
	# 1.0f0 has type Float32
	if value < 0 || value > 1
		error("value must be within [0,1]")
	end

	# fix those three functions such that the property is achieved
	# maybe think about periodic functions ;)
	# only function definitions!
	red_map(x) = x^2 # red_map is a function
	# you call by: red_map(0.123)
	green_map(x) = x^3 # TODO
	blue_map(x) = x^4 # TODO

	# that's a function call of red_map
	# red_map(value)

	return (3.0, 4f0, 2) # TODO, return the correct tuple! You should call
	# red_map, green_map, blue_map
end

# ╔═╡ 2c9d77c0-a5e0-4fe6-aa26-dcb3832904ff
red_map(x) = x

# ╔═╡ fe1b69e1-966d-49fe-bcc8-63781aa334ba
red_map(0.123)

# ╔═╡ 719c42e4-44f3-4519-a107-76b3d64717e7

# ╔═╡ 0d53abcb-281a-4e29-b695-ea1d11520d09
md"## Colors are only a datatype"

# ╔═╡ 38354d23-6a58-49bc-a19b-0a4b0186b8c1
@bind g Slider(0:0.01:1, show_value=true)

# ╔═╡ 6ba803f5-bb30-4d9b-a96f-fc23fd8fa8c8
@bind r Slider(0:0.01:1, show_value=true)

# ╔═╡ 76cdfdb7-9ca9-49d3-9ec9-0d26edbd458f
@bind b Slider(0:0.01:1, show_value=true)

# ╔═╡ 3896942f-36fe-4326-a4a0-1a502610e451
t = (r, g, b)

# ╔═╡ 63a0fc2c-3171-4dac-9eff-2f7bd2dce253
RGB(t[1], t[2], t[3])

# ╔═╡ 1b2e1d30-6805-40a2-9ce8-113e5664e434
typeof(RGB(t...))

# ╔═╡ 817bbf4d-7856-4269-84bb-b3344345f9fb

# ╔═╡ fb13e2f5-62a3-4b87-bc78-44f3e8dbef7f
md"## Documentation"

# ╔═╡ 3931594b-2353-4950-83e1-5ce163e69e64
"""
	f(x)

Returns `x+1`

## Examples

```julia-repl
julia> f(10)
11

```
"""
function f(x)
	return x + 1
end

# ╔═╡ a84872cf-b435-4806-b016-2a5fa71214e9
f

# ╔═╡ e9cdd216-17dd-4e74-9f94-f62b9fb1160f

# ╔═╡ 264ae342-b0f7-4a3d-b309-377fc89f1595
md"# Sampling

## Nyquist Theorem
You need to sample a function with twice the highest occuring frequency in the signal itself.

## Example
What is the highest frequency in the function s?

How do we calculate the frequency?

1 divided by the periodicity of the function

$$f_\text{max} = \frac{1}{T}$$

The frequency of our function s is

$$f_\text{max} = \frac{1}{0.5s} = 2 Hz$$

The Nyquist frequency for s would be now $>4Hz$

(The real world consists of many, many frequencies)

"

# ╔═╡ a149d945-6ed5-40f9-b7c4-008181cd49ab
s(x) = sin(x * 2 * pi * 2)

# ╔═╡ 935556e9-a5f0-4676-b275-bb756d64979e
md"

"

# ╔═╡ b10b23fe-6aba-45b9-85e0-2fb44a6034c2
N = 1000

# ╔═╡ 648394b7-f059-443f-9308-12a26eac4f3b

# ╔═╡ 00eeb89c-9ba5-49b9-b1a5-be5a0f7c27eb
begin
	xs = range(0, 5, length=N)
	ys = s.(xs)
end

# ╔═╡ 483b5371-7e5d-42e9-abd1-64128d636090
# only a frequency Hz
bad_sampling_freq = range(0, 5, length=7)

# ╔═╡ 041b3e96-a866-41f5-a7d1-46294460c7ff
nyquist_sampling_freq = range(0, 5, length=22)

# ╔═╡ 64db368a-8fa0-4315-b27a-ad0628ff8e77
begin
	plot(xs, ys, xlabel="t in seconds")
	plot!(bad_sampling_freq, s.(bad_sampling_freq), mark="*")
	plot!(nyquist_sampling_freq, s.(nyquist_sampling_freq), mark="*")
end

# ╔═╡ cd6708f5-9a2a-4139-81b1-14381abbad73

# ╔═╡ 0d19eaca-d1f5-40e3-906e-caa82ad006b4
md"## Fourier Transform
Transforms a signal into Fourier space (Frequency space).

Fourier transform simply interprets the signal differently (or it displays it differently).
"

# ╔═╡ Cell order:
# ╠═0900dd0b-483a-404d-bf84-43c075230b73
# ╠═5cdc016a-b0aa-42b4-841c-657c8bf72ca3
# ╠═3d47cb3e-43b3-11ec-1b61-992ecbb3d16f
# ╠═aca53457-93bb-4e6c-a3dc-68c000e8db78
# ╠═e1f76c58-fa78-4a20-b30f-5ef087874899
# ╠═c5e12406-97b2-450d-9fe7-b2da82452f33
# ╠═73f690da-3ee9-4a27-8e97-b062a615266e
# ╠═1f73e297-6a82-4d33-8862-5a166da4dc2f
# ╠═fcbb40e5-87be-49f8-968b-e8bdd962ee53
# ╠═a564a1d9-0c5f-4391-bae4-04ffc9abc861
# ╠═c4087dcf-ea75-4b57-b4cc-73934f509af8
# ╠═4c4199cf-a62f-4b82-ac04-0b1cb9b2543f
# ╠═fc992ffe-1f4e-45e8-b2f2-7e4ec53d2188
# ╠═26eb3fff-ece9-40c9-bd9d-32cba04746a1
# ╠═09c7396d-57d4-4e26-a0bf-7b5730836135
# ╠═2c9d77c0-a5e0-4fe6-aa26-dcb3832904ff
# ╠═fe1b69e1-966d-49fe-bcc8-63781aa334ba
# ╠═719c42e4-44f3-4519-a107-76b3d64717e7
# ╠═0d53abcb-281a-4e29-b695-ea1d11520d09
# ╠═3896942f-36fe-4326-a4a0-1a502610e451
# ╠═38354d23-6a58-49bc-a19b-0a4b0186b8c1
# ╠═6ba803f5-bb30-4d9b-a96f-fc23fd8fa8c8
# ╠═76cdfdb7-9ca9-49d3-9ec9-0d26edbd458f
# ╠═63a0fc2c-3171-4dac-9eff-2f7bd2dce253
# ╠═1b2e1d30-6805-40a2-9ce8-113e5664e434
# ╠═817bbf4d-7856-4269-84bb-b3344345f9fb
# ╠═fb13e2f5-62a3-4b87-bc78-44f3e8dbef7f
# ╠═3931594b-2353-4950-83e1-5ce163e69e64
# ╠═a84872cf-b435-4806-b016-2a5fa71214e9
# ╠═e9cdd216-17dd-4e74-9f94-f62b9fb1160f
# ╠═264ae342-b0f7-4a3d-b309-377fc89f1595
# ╠═a149d945-6ed5-40f9-b7c4-008181cd49ab
# ╠═935556e9-a5f0-4676-b275-bb756d64979e
# ╠═b10b23fe-6aba-45b9-85e0-2fb44a6034c2
# ╠═648394b7-f059-443f-9308-12a26eac4f3b
# ╠═00eeb89c-9ba5-49b9-b1a5-be5a0f7c27eb
# ╠═483b5371-7e5d-42e9-abd1-64128d636090
# ╠═041b3e96-a866-41f5-a7d1-46294460c7ff
# ╠═64db368a-8fa0-4315-b27a-ad0628ff8e77
# ╠═cd6708f5-9a2a-4139-81b1-14381abbad73
# ╠═0d19eaca-d1f5-40e3-906e-caa82ad006b4


text:
export compute_cc_c3

"""
    compute_cc(FFT1::FFTData, FFT2::FFTData, maxlag::Float64;
               smoothing_half_win::Int=20,
               corr_type::String="cross-correlation" )

"""
function compute_cc_c3(FFT1::FFTData, FFT2::FFTData, maxlag::Float64;
                    smoothing_half_win::Int=20,
                    corr_type::String="cross-correlation")

    N = convert(Int,round(FFT1.cc_len * FFT1.fs)) # number of data points
    comp = FFT1.name[end] * FFT2.name[end]

    # assume perfect overlap of dates for this cross-correlation
    corr = correlate(FFT1.fft, FFT2.fft, N,
                     convert(Int,round(maxlag * FFT1.fs)),
                     corr_type=corr_type)
    rotated = false

    return CorrData(FFT1, FFT2, comp, rotated, corr_type,
                    maxlag, FFT1.t, corr)

end

max_stars_repo_name: kura-okubo/SeisXcorrelation

max_issues_repo_name: kura-okubo/SeisXcorrelation

max_forks_repo_name: kura-okubo/SeisXcorrelation


text:
module LifeValGenerate

using Random, Distributions

export Gen

Random.seed!(2020)

function gen_normal(meanv, sdv, minv, maxv, fact)
    if meanv==99
        val = 0
    else
        x = rand(Normal(meanv,sdv),1)
        y=x[1]
        val = round(y, digits=0)
        if val < minv
            val=minv
        elseif val > maxv
            val=maxv
        end
    end
    return (val * fact)
end

function gen_beta(alphav, betav, minv, maxv, fact)
    x = rand(Beta(alphav,betav),1)
    y=x[1]*fact
    val = round(y, digits=-3)
    if val < minv
        val=minv
    elseif val > maxv
        val=maxv
    end
    return (val)
end

#the Gen function below generates 2 sets of output, the IF book and the at inception book
#code is crudely written!  our appologies, to be revised but not top of list!
function Gen(v::Vector)
    policies = Array{Int64}(undef, v[4], 17)
    prems = Array{Int64}(undef, v[4], 17)

    for n in 1:v[4]
        #generate from Normal and Beta distributions!
        age = gen_normal(v[6],v[7],v[8],v[9],1)
        sa = gen_beta(v[11],v[12],v[13],v[14],v[23])
        term =gen_normal(v[16],v[17],v[18],v[19],1)

        if v[19] != 99  #term_if, ceiling if term applicable
            term_if =gen_normal(v[20],v[21],v[22],term-1,1) #with ceiling
        else
            term_if =gen_normal(v[20],v[21],v[22],99,1) #no ceiling
        end

        #alloc to policies in force
        policies[n,1]=n   #id
        policies[n,2]=v[1] #prod
        policies[n,3]=age #age
        policies[n,4]=sa #sumassured
        policies[n,5]=term #term
        policies[n,6]=term_if
        #current age...
        policies[n,3]=policies[n,3]+policies[n,6] #age is current age for reserve calc
        #terms are in months!
        policies[n,5]=term*12  #term converted to months
        policies[n,6]=term_if*12  #term_if converted to months
        policies[n,7]=0
        policies[n,8]=rand(0:1)  #gender, 0 males otherwise 1 females
        policies[n,9]=v[2] #spcode
        policies[n,10]=0
        policies[n,11]=0
        policies[n,12]=rand(1:12) #prem-inc month
        policies[n,13]=0
        policies[n,14]=0
        policies[n,15]=0
        policies[n,16]=0
        policies[n,17]=0

        prems[n,1]=n   #id
        prems[n,2]=v[1] #prod
        prems[n,3]=age #age
        prems[n,4]=sa #sumassured
        prems[n,5]=term*12 #term
        prems[n,6]=0
        prems[n,7]=0
        prems[n,8]=rand(0:1)  #gender, 0 males otherwise 1 females
        prems[n,9]=v[2] #spcode
        prems[n,10]=0
        prems[n,11]=0
        prems[n,12]=rand(1:12) #prem-inc month
        prems[n,13]=0
        prems[n,14]=0
        prems[n,15]=0
        prems[n,16]=0
        prems[n,17]=0
    end

    return [policies, prems]
end

end # module


text:
######################################################################
# General purpose function stubs for specialization
# -----
# Licensed under MIT License
export use
use(args...; kwargs...) = MethodError(use, (kwargs, args...))

export cancel
cancel(args...; kwargs...) = MethodError(cancel, (kwargs, args...))

export clear
clear(args...; kwargs...) = MethodError(clear, (kwargs, args...))

export update!
update!(args...; kwargs...) = MethodError(update!, (kwargs, args...))

export init
init(args...; kwargs...) = MethodError(init, (kwargs, args...))

export load
load(args...; kwargs...) = MethodError(load, (kwargs, args...))

export store
store(args...; kwargs...) = MethodError(store, (kwargs, args...))

export restore
restore(args...; kwargs...) = MethodError(restore, (kwargs, args...))


text:
export combinations,
       CoolLexCombinations,
       multiset_combinations,
       with_replacement_combinations,
       powerset

#The Combinations iterator

struct Combinations{T}
    a::T
    t::Int
end

function Base.iterate(c::Combinations, s = collect(1:c.t))
    (!isempty(s) && s[1] > length(c.a) - c.t + 1) && return

    comb = [c.a[si] for si in s]
    if c.t == 0
        # special case to generate 1 result for t==0
        return (comb, [length(c.a)+2])
    end
    s = copy(s)
    for i = length(s):-1:1
        s[i] += 1
        if s[i] > (length(c.a) - (length(s) - i))
            continue
        end
        for j = i+1:lastindex(s)
            s[j] = s[j-1]+1
        end
        break
    end
    (comb, s)
end

Base.length(c::Combinations) = binomial(length(c.a), c.t)

Base.eltype(::Type{Combinations{T}}) where {T} = Vector{eltype(T)}

"""
    combinations(a, n)

Generate all combinations of `n` elements from an indexable object `a`. Because the number
of combinations can be very large, this function returns an iterator object.
Use `collect(combinations(a, n))` to get an array of all combinations.
"""
function combinations(a, t::Integer)
    if t < 0
        # generate 0 combinations for negative argument
        t = length(a) + 1
    end
    Combinations(a, t)
end

"""
    combinations(a)

Generate combinations of the elements of `a` of all orders. Chaining of order iterators
is eager, but the sequence at each order is lazy.
"""
combinations(a) = Iterators.flatten([combinations(a, k) for k = 1:length(a)])

# cool-lex combinations iterator

"""
    CoolLexCombinations

Produce ``(n,k)``-combinations in cool-lex order.

# Reference

Ruskey, F., & Williams, A. (2009). The coolest way to generate combinations.
*Discrete Mathematics*, 309(17), 5305-5320.
"""
struct CoolLexCombinations
    n::Int
    t::Int
end

struct CoolLexIterState{T<:Integer}
    R0::T
    R1::T
    R2::T
    R3::T
end

function Base.iterate(C::CoolLexCombinations)
    if C.n < 0
        throw(DomainError(C.n))
    end
    if C.t ≤ 0
        throw(DomainError(C.t))
    end

    #What integer size should I use?
    if C.n < 8sizeof(Int)
        T = Int
    else
        T = BigInt
    end

    state = CoolLexIterState{T}(0, 0, T(1) << C.n, (T(1) << C.t) - 1)
    iterate(C, state)
end

function Base.iterate(C::CoolLexCombinations, S::CoolLexIterState)
    (S.R3 & S.R2 != 0) && return

    R0 = S.R0
    R1 = S.R1
    R2 = S.R2
    R3 = S.R3

    R0 = R3 & (R3 + 1)
    R1 = xor(R0, R0 - 1)
    R0 = R1 + 1
    R1 &= R3
    R0 = max((R0 & R3) - 1, 0)
    R3 += R1 - R0

    _cool_lex_visit(S.R3), CoolLexIterState(R0, R1, R2, R3)
end

#Converts an integer bit pattern X into a subset
#If X & 2^k == 1, then k is in the subset
function _cool_lex_visit(X::Integer)
    subset = Int[]
    n = 1
    while X != 0
        X & 1 == 1 && push!(subset, n)
        X >>= 1
        n += 1
    end
    subset
end

Base.length(C::CoolLexCombinations) = max(0, binomial(C.n, C.t))

struct MultiSetCombinations{T}
    m::T
    f::Vector{Int}
    t::Int
    ref::Vector{Int}
end

Base.eltype(::Type{MultiSetCombinations{T}}) where {T} = Vector{eltype(T)}

function Base.length(c::MultiSetCombinations)
    t = c.t
    if t > length(c.ref)
        return 0
    end
    p = [1; zeros(Int, t)]
    for i in 1:length(c.f)
        f = c.f[i]
        if i == 1
            for j in 1:min(f, t)
                p[j+1] = 1
            end
        else
            for j in t:-1:1
                p[j+1] = sum(p[max(1,j+1-f):(j+1)])
            end
        end
    end
    return p[t+1]
end

function multiset_combinations(m, f::Vector{<:Integer}, t::Integer)
    length(m) == length(f) || error("Lengths of m and f are not the same.")
    ref = length(f) > 0 ? vcat([[i for j in 1:f[i] ] for i in 1:length(f)]...) : Int[]
    if t < 0
        t = length(ref) + 1
    end
    MultiSetCombinations(m, f, t, ref)
end

"""
    multiset_combinations(a, t)

Generate all combinations of size `t` from an array `a` with possibly duplicated elements.
"""
function multiset_combinations(a, t::Integer)
    m = unique(collect(a))
    f = Int[sum([c == x for c in a]) for x in m]
    multiset_combinations(m, f, t)
end

function Base.iterate(c::MultiSetCombinations, s = c.ref)
    ((!isempty(s) && max(s[1], c.t) > length(c.ref)) || (isempty(s) && c.t > 0)) && return

    ref = c.ref
    n = length(ref)
    t = c.t
    changed = false
    comb = [c.m[s[i]] for i in 1:t]
    if t > 0
        s = copy(s)
        for i in t:-1:1
            if s[i] < ref[i + (n - t)]
                j = 1
                while ref[j] <= s[i]
                    j += 1
                end
                s[i] = ref[j]
                for l in (i+1):t
                    s[l] = ref[j+=1]
                end
                changed = true
                break
            end
        end
        !changed && (s[1] = n+1)
    else
        s = [n+1]
    end
    (comb, s)
end

struct WithReplacementCombinations{T}
    a::T
    t::Int
end

Base.eltype(::Type{WithReplacementCombinations{T}}) where {T} = Vector{eltype(T)}

Base.length(c::WithReplacementCombinations) = binomial(length(c.a) + c.t - 1, c.t)

"""
    with_replacement_combinations(a, t)

Generate all combinations with replacement of size `t` from an array `a`.
"""
with_replacement_combinations(a, t::Integer) = WithReplacementCombinations(a, t)

function Base.iterate(c::WithReplacementCombinations, s = [1 for i in 1:c.t])
    (!isempty(s) && s[1] > length(c.a) || c.t < 0) && return

    n = length(c.a)
    t = c.t
    comb = [c.a[si] for si in s]
    if t > 0
        s = copy(s)
        changed = false
        for i in t:-1:1
            if s[i] < n
                s[i] += 1
                for j in (i+1):t
                    s[j] = s[i]
                end
                changed = true
                break
            end
        end
        !changed && (s[1] = n+1)
    else
        s = [n+1]
    end
    (comb, s)
end

## Power set

"""
    powerset(a, min=0, max=length(a))

Generate all subsets of an indexable object `a` including the empty set, with cardinality
bounded by `min` and `max`. Because the number of subsets can be very large, this function
returns an iterator object. Use `collect(powerset(a, min, max))` to get an array of all
subsets.
"""
function powerset(a, min::Integer=0, max::Integer=length(a))
    itrs = [combinations(a, k) for k = min:max]
    min < 1 && append!(itrs, eltype(a)[])
    Iterators.flatten(itrs)
end


text:
struct TDVPInfo
  maxtruncerr::Float64
end


text:
######## Stan diagnose example  ###########

using StanSample, Test

ProjDir = @__DIR__

gq = "
  data {
    int<lower=0> N;
    int<lower=0> y[N];
  }
  parameters {
    real<lower=0.00001> Theta;
  }
  model {
    Theta ~ gamma(4, 1000);
    for (n in 1:N)
      y[n] ~ exponential(Theta);
  }
  generated quantities{
    real y_pred;
    y_pred = exponential_rng(Theta);
  }
";

gq_data = Dict(
  "N" => 3,
  "y" => [100, 950, 450]
);

stanmodel = SampleModel("generate_quantities", gq);

rc = stan_sample(stanmodel; data=gq_data)

if success(rc)
  # Convert to an MCMCChains.Chains object
  samples = read_samples(stanmodel)

  # Show the same output in DataFrame format
  df = StanSample.read_summary(stanmodel)

  df1 = stan_generate_quantities(stanmodel, 1)
  @test sum(df1[!, :y_pred])/size(df1, 1) ≈ 412.0 rtol=0.3
end


text:
using OrdinaryDiffEq, Test

f(u,p,t) = 2u
u0 = 0.5
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob,Tsit5())

function f(du, u, p, t)
    du .= 2.0 * u
end
prob = ODEProblem(f,u0,tspan)
@test_throws DiffEqBase.IncompatibleInitialConditionError sol = solve(prob,Tsit5())

prob = ODEProblem{false}(f,u0,tspan)
sol = solve(prob,Tsit5())

@test_throws DiffEqBase.ProblemSolverPairingError solve(prob, DFBDF())
@test_throws DiffEqBase.NoDefaultAlgorithmError solve(prob,nothing)
@test_throws DiffEqBase.NonSolverError solve(prob,5.0)

prob = ODEProblem{false}(f,u0,(nothing,nothing))
@test_throws DiffEqBase.NoTspanError solve(prob,Tsit5())

max_stars_repo_name: JuliaDiffEq/DifferentialEquationsBase

max_issues_repo_name: JuliaDiffEq/DifferentialEquationsBase

max_forks_repo_name: JuliaDiffEq/DifferentialEquationsBase