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

using Markdown
using InteractiveUtils

# ╔═╡ 7c78684c-5608-4184-b7ca-9fe42c09f3f9
begin
  using PlutoUI
  #using LaTeXStrings
  using TikzPictures
end

# ╔═╡ d8aa0d5c-273f-11ec-3d4e-d9c3dbb94278
md"""
# B.5 Trees
By definition, a _**free tree**_ is an undirected, connected, acyclic graph. At first glance, this seems to be
an object with many (`3` to be precise) adjectives. But, once we started to have a deeper understanding, we
might come to realize that a free tree is just an undirected graph that is
```math
\newcommand{\norm}[1]{\lVert{#1}\rVert}
\newcommand{\abs}[1]{\lvert{#1}\rvert}
```

- connected
  - (Exercise B.4-3) If ``G = (V, E)`` is an undirected, connected graph then ``\abs{E} \ge \abs{V} - 1``.
  - In day-to-day English, this simply says,
    > "For an undirected graph to be connected, there cannot be too few edges; there must exist at least ``\abs{V} - 1`` edges."
- acyclic
  - We prove an analogous result: If ``G = (V, E)`` is an undirected, acyclic graph then ``\abs{E} \le \abs{V} - 1``.
  - This simply says analogously,
    > "For an undirected graph to remain acyclic, there cannot be too many edges. Too many edges create cycles. There can be at most ``\abs{V} - 1`` edges."

_Proof._$(HTML("
"))
We show that for any acyclic undirected graph ``G = (V, E)``, we have ``\abs{E} \le \abs{V} - 1``.

Note that this is not saying that cyclic graphs tend to have many more edges than acyclic graphs do. Indeed, three edges btw three distinct vertices suffice to make a graph cyclic. I find it useful to picture in our head like this:
> We are given the vertex set ``V`` and an empty edge set ``E``.
> We search to add as many edges to ``E`` as possible in keeping ``G`` acyclic.

Again, we shall prove this by mathematical induction on ``\abs{V}``. It is easy to
verify the validity of the statement when ``\abs{V} = 1, 2``. When ``\abs{V} = 3``, assume that ``G`` is acyclic. If ``\abs{E} = 3``, then ``G`` contains a cycle; in fact, in this case, ``G`` is itself a cycle. So we can conclude ``\abs{E} \le 2 = \abs{V} - 1.`` (Three vertices admit at most ``\begin{pmatrix} 3 \\ 2\end{pmatrix} = 3`` edges, and,  since ``3`` is excluded, we must have ``\abs{E} \le 2``.)

Now, assume that the statement holds for ``\abs{V} = 1, 2, 3, \ldots, k``. Let ``G = (V, E)`` be an acyclic undirected graph with ``\abs{V} = k+1``.
First, we would like to
01. Find some vertex ``u \in V``
02. Exclude ``u`` from ``V`` and ``E``. That is, obtain a new graph ``G_{1} = (V_{1}, E_{1})``, subgraph of ``G``, such that ``u`` is neither in ``V_{1}`` nor in any edge from ``E_{1}``
03. Argue that ``\abs{E} = \abs{E_{1}} + 1 \le (\abs{V_{1}} - 1) + 1 = \abs{V} - 1``

We soon realize that not every vertex of ``V`` can play the role of ``u`` above. Indeed, for the third argument above to work, we need ``\deg(u) = 1``.
In other words, we would like to pick ``u`` such that it is adjent to only one other vertex. Now, we give the argument that such ``u`` exists. More precisely, we show that
```math
\exists\, u \in V \;\text{s.t.}\; \deg(u) \le 1.
```
This can be shown using proof by contradiction: Assume the contrary, i.e. that every vertex in ``V`` has degree greater than or equal to ``2``. Then starting from any vertex ``v`` of ``V``, which we put at the left of the following figure, we will always arrive at the conclusion that there exists a cycle in ``V``.

Once the existence of ``u`` s.t. ``\deg(u) \le 1`` is established, it is easy to prove the rest.

- If ``\deg(u) = 0``, then removing ``u`` makes the aforementioned ``G_{1} = (V_{1}, E_{1})`` an acyclic graph of ``k`` vertices, whence
  ```math
  \abs{E} = \abs{E_{1}} \le (\abs{V} - 1) - 1 \le \abs{V} - 1.
  ```
- Similarly, if ``\deg(u) = 1``, then removing ``u`` makes ``G_{1}`` an acyclic graphs of ``k`` vertices, whence
  ```math
  \abs{E} = \abs{E_{1}} + 1 \le (\abs{V_{1}} - 1) + 1 = \abs{V_{1}} = \abs{V} - 1.
  ```
This concludes the proof.
"""

# ╔═╡ 684a24c0-9ca9-46d6-923f-f35c9129d221

# ╔═╡ 1a24de70-4203-475d-b183-086c6e0e392c

# ╔═╡ 4e5a7a11-397a-4f02-835f-cf56ff8b2ead

# ╔═╡ 01a0665d-794c-4f72-80cc-237c59753f94

# ╔═╡ 638c65ab-3eb4-416b-8adf-110649dac5c4
md"""
## Theorem of The Section
Since a tree (i.e. a free tree. It seems that we often drop the adjective and speak of trees as free trees) is at the same time connected and acyclic, from the above relation btw number of edges and number of vertices, we should have
```math
\abs{E} = \abs{V} - 1.
```

**Theorem B.2 (Properties of free trees)**$(HTML("
"))
The followings are equiv:
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
TikzPictures = "37f6aa50-8035-52d0-81c2-5a1d08754b2d"

[compat]
PlutoUI = "~0.7.14"
TikzPictures = "~3.4.1"
"""

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[HypertextLiteral]]
git-tree-sha1 = "f6532909bf3d40b308a0f360b6a0e626c0e263a8"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.1"

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

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

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

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

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

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

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

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

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

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

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

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

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

[[LittleCMS_jll]]
deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pkg"]
git-tree-sha1 = "110897e7db2d6836be22c18bffd9422218ee6284"
uuid = "d3a379c0-f9a3-5b72-a4c0-6bf4d2e8af0f"
version = "2.12.0+0"

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

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

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

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

[[OpenJpeg_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libtiff_jll", "LittleCMS_jll", "Pkg", "libpng_jll"]
git-tree-sha1 = "76374b6e7f632c130e78100b166e5a48464256f8"
uuid = "643b3616-a352-519d-856d-80112ee9badc"
version = "2.4.0+0"

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

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

[[Pixman_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b4f5d02549a10e20780a24fce72bea96b6329e29"
uuid = "30392449-352a-5448-841d-b1acce4e97dc"
version = "0.40.1+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 = ["Base64", "Dates", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "d1fb76655a95bf6ea4348d7197b22e889a4375f4"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.14"

[[Poppler_jll]]
deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "OpenJpeg_jll", "Pkg", "libpng_jll"]
git-tree-sha1 = "e11443687ac151ac6ef6699eb75f964bed8e1faa"
uuid = "9c32591e-4766-534b-9725-b71a8799265b"
version = "0.87.0+2"

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

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

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

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

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

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

[[Tectonic]]
deps = ["Pkg"]
git-tree-sha1 = "acf12eccb390a78653ee805cd527898f01f78a85"
uuid = "9ac5f52a-99c6-489f-af81-462ef484790f"
version = "0.6.1"

[[TikzPictures]]
deps = ["LaTeXStrings", "Poppler_jll", "Requires", "Tectonic"]
git-tree-sha1 = "a08671c0979063a437378f6410bb75a465f3cd1c"
uuid = "37f6aa50-8035-52d0-81c2-5a1d08754b2d"
version = "3.4.1"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# ╔═╡ Cell order:
# ╠═7c78684c-5608-4184-b7ca-9fe42c09f3f9
# ╠═d8aa0d5c-273f-11ec-3d4e-d9c3dbb94278
# ╠═684a24c0-9ca9-46d6-923f-f35c9129d221
# ╠═1a24de70-4203-475d-b183-086c6e0e392c
# ╠═4e5a7a11-397a-4f02-835f-cf56ff8b2ead
# ╠═01a0665d-794c-4f72-80cc-237c59753f94
# ╟─638c65ab-3eb4-416b-8adf-110649dac5c4
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
import NLsolve: IsFiniteException

@testset "throws" begin

function f_inf!(F, x)
    copyto!(F, x)
    F[1] = Inf
    return F
end

function f_nan!(F, x)
    copyto!(F, x)
    F[1] = NaN
    return F
end

@test_throws IsFiniteException nlsolve(f_inf!, [ -0.5; 1.4], method = :trust_region, autodiff=true)
@test_throws IsFiniteException nlsolve(f_inf!, [ -0.5; 1.4], method = :newton, autodiff=true)

@test_throws IsFiniteException nlsolve(f_nan!, [ -0.5; 1.4], method = :trust_region, autodiff=true)
@test_throws IsFiniteException nlsolve(f_nan!, [ -0.5; 1.4], method = :newton, autodiff=true)

end # testset

ext: jl


text:
include("../../matplotlib_settings.jl")

set_config(:tab_compensated_summation, false)

using ChargedParticleDynamics.GuidingCenter4d: plot_integral_error, plot_poincare_loop, plot_poincare_trajectories

const coll_title = "Guiding Center 4D Symmetric 1st Poincaré Integral Invariant"

const XMIN = -0.5
const XMAX = +0.5
const YMIN = -0.5
const YMAX = +0.5

include("../guiding_center_4d_poincare_invariant_1st.jl")


text:
import Base.sync_varname
import Base.@async

macro _async(expr)
    thunk = esc(:(()->($expr)))
    var = esc(sync_varname)
    quote
        local task = Task($thunk)
        if $(Expr(:isdefined, var))
            push!($var, task)
        end
        schedule(task)
        task
    end
end

macro async(expr)
    quote
        @_async try
            $expr
        catch err
            @warn "error within async" exception=err
        end
    end
end


text:
# Generative model using the change of variables formula
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using LinearAlgebra, InvertibleNetworks, PyPlot, Flux, Random
import Flux.Optimise.update!

# Random seed
Random.seed!(11)

####################################################################################################

# Define network
nx = 1
ny = 1
n_in = 2
n_hidden = 64
batchsize = 20
depth = 4
AN = Array{ActNorm}(undef, depth)
L = Array{CouplingLayerGlow}(undef, depth)
Params = Array{Parameter}(undef, 0)

# Create layers
for j=1:depth
    AN[j] = ActNorm(n_in; logdet=true)
    L[j] = CouplingLayerGlow(n_in, n_hidden; k1=1, k2=1, p1=0, p2=0, logdet=true)

    # Collect parameters
    global Params = cat(Params, get_params(AN[j]); dims=1)
    global Params = cat(Params, get_params(L[j]); dims=1)
end

# Forward pass
function forward(X)
    logdet = 0f0
    for j=1:depth
        X_, logdet1 = AN[j].forward(X)
        X, logdet2 = L[j].forward(X_)
        logdet += (logdet1 + logdet2)
    end
    return X, logdet
end

# Backward pass
function backward(ΔX, X)
    for j=depth:-1:1
        ΔX_, X_ = L[j].backward(ΔX, X)
        ΔX, X = AN[j].backward(ΔX_, X_)
    end
    return ΔX, X
end

####################################################################################################

# Loss
function loss(X)
    Y, logdet = forward(X)
    f = -log_likelihood(Y) - logdet
    ΔY = -∇log_likelihood(Y)
    ΔX = backward(ΔY, Y)[1]
    return f, ΔX
end

# Training
maxiter = 1000
opt = Flux.ADAM(1f-3)
fval = zeros(Float32, maxiter)

for j=1:maxiter

    # Evaluate objective and gradients
    X = sample_banana(batchsize)
    fval[j] = loss(X)[1]
    mod(j, 10) == 0 && (print("Iteration: ", j, "; f = ", fval[j], "\n"))

    # Update params
    for p in Params
        update!(opt, p.data, p.grad)
    end
    clear_grad!(Params)
end

####################################################################################################

# Testing
test_size = 500
X = sample_banana(test_size)
Y_ = forward(X)[1]
Y = randn(Float32, 1, 1, 2, test_size)
X_ = backward(Y, Y)[2]

# Plot
figure(figsize=[8,8])
ax1 = subplot(2,2,1); plot(X[1, 1, 1, :], X[1, 1, 2, :], "."); title(L"Data space: $x \sim \hat{p}_X$")
ax1.set_xlim([-3.5,3.5]); ax1.set_ylim([0,50])
ax2 = subplot(2,2,2); plot(Y_[1, 1, 1, :], Y_[1, 1, 2, :], "g."); title(L"Latent space: $z = f(x)$")
ax2.set_xlim([-3.5, 3.5]); ax2.set_ylim([-3.5, 3.5])
ax3 = subplot(2,2,3); plot(X_[1, 1, 1, :], X_[1, 1, 2, :], "g."); title(L"Data space: $x = f^{-1}(z)$")
ax3.set_xlim([-3.5,3.5]); ax3.set_ylim([0,50])
ax4 = subplot(2,2,4); plot(Y[1, 1, 1, :], Y[1, 1, 2, :], "."); title(L"Latent space: $z \sim \hat{p}_Z$")
ax4.set_xlim([-3.5, 3.5]); ax4.set_ylim([-3.5, 3.5])


text:
function __init__()
    RI_INFO_ROOT[] = joinpath(artifact"refractiveindex.info",
                              "refractiveindex.info-database-2020-01-19", "database")

    lib = YAML.load_file(joinpath(RI_INFO_ROOT[], "library.yml"), dicttype=Dict{String, Any})
    for shelf in lib
        shelfname = shelf["SHELF"]
        for book in shelf["content"]
            haskey(book, "DIVIDER") && continue
            bookname = book["BOOK"]
            for page in book["content"]
                haskey(page, "DIVIDER") && continue
                pagename = string(page["PAGE"])
                RI_LIB[(shelfname, bookname, pagename)] = (name = page["name"], path=page["data"])
            end
        end
    end

end


text:
using Erjulix
using Documenter

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

makedocs(;
    modules=[Erjulix],
    authors="Paul Bayer",
    repo="https://github.com/pbayer/erjulix/blob/{commit}{path}#{line}",
    sitename="erjulix",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://pbayer.github.io/erjulix",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/pbayer/erjulix",
)

max_stars_repo_name: pbayer/erjulix

max_issues_repo_name: pbayer/erjulix

max_forks_repo_name: pbayer/erjulix


text:
# ArgMacros
# handleast.jl

"Extract name of macro and change from `ArgMacros.@...` to `@...`"
function _get_macroname(arg::Expr)::Symbol
    return arg.args[1] isa Expr && arg.args[1].head == :. && arg.args[1].args[1] == :ArgMacros ?
                arg.args[1].args[2].value :
                arg.args[1]
end

"Get the macrocall expressions from a block as a generator"
function _getmacrocalls(block::Expr)
    Iterators.filter(arg -> arg isa Expr && arg.head == :macrocall, block.args)
end

"""
Enforce argument declaration ordering:
    Flagged → Required Positional → Optional Positional

Throw `ArgumentError` if ordering violates this rule.
"""
function _validateorder(block::Expr)
    encountered_positional = false
    encoundered_optional_positional = false

    for arg in _getmacrocalls(block)
        # Fix namespace issues
        macroname::Symbol = _get_macroname(arg)

        if macroname in FLAGGED_SYMBOLS
            if encountered_positional
                throw(ArgumentError(
                    "Positional arguments must be declared after all flagged arguments.\nFrom: $arg"
                ))
            end
        elseif macroname == POSITIONAL_REQUIRED_SYMBOL
            encountered_positional = true

            if encoundered_optional_positional
                throw(ArgumentError(
                    "Required positional arguments must be declared in order before all optional positional arguments.\nFrom: $arg"
                ))
            end
        elseif macroname in POSITIONAL_OPTIONAL_SYMBOLS
            encountered_positional = true
            encoundered_optional_positional = true
        end
    end
end

"Extract name => type pair from argument macrocall expression"
function _getargumentpair(arg::Expr)::Union{Expr, Nothing}
    macroname::Symbol = _get_macroname(arg)

    if macroname in (ARGUMENT_REQUIRED_SYMBOL, POSITIONAL_REQUIRED_SYMBOL)
        :($(arg.args[4])::$(arg.args[3]))
    elseif macroname in (ARGUMENT_DEFAULT_SYMBOL, POSITIONAL_DEFAULT_SYMBOL)
        :($(arg.args[5])::$(arg.args[3]))
    elseif macroname in (ARGUMENT_OPTIONAL_SYMBOL, POSITIONAL_OPTIONAL_SYMBOL)
        :($(arg.args[4])::Union{$(arg.args[3]), Nothing})
    elseif macroname == ARGUMENT_FLAG_SYMBOL
        :($(arg.args[3])::Bool)
    elseif macroname == ARGUMENT_COUNT_SYMBOL
        :($(arg.args[3])::Int)
    end
end

"Extract name => type pairs for all argument macros in block"
function _getargumentpairs(block::Expr)
    Iterators.filter(!isnothing, _getargumentpair(arg) for arg in _getmacrocalls(block))
end


text:
"""
    Fermi.MollerPlesset

module for running MP2 energies on restricted and unrestricted HF references.
## methods
    do_rmp2 -> see docstring ?Fermi.MollerPlesset.do_rmp2
    do_ump2 -> see docstring ?Fermi.MollerPlesset.do_ump2
"""
module MollerPlesset

using Fermi
#using Fermi.Direct
using Fermi.Output
using TensorOperations

export RMP2
export RMP3

function print_header()
    @output "================================================================================\n"
    @output "|   Moller-Plesset Perturbation Theory                                         |\n"
    @output "|       module written by M.M. Davis                                           |\n"
    @output "================================================================================\n"
end

abstract type MP2Algorithm end

struct Conventional <: MP2Algorithm end
struct DF           <: MP2Algorithm end
struct Direct       <: MP2Algorithm end

abstract type AbstractMPWavefunction <: Fermi.AbstractCorrelatedWavefunction end

include("RMP2.jl")
include("RMP3.jl")

end #module


text:
function save(
	var   :: AbstractArray{Float32,3},
	dt    :: TimeType,
	npd   :: NASAPrecipitationDataset,
	geo   :: GeoRegion,
	ginfo :: RegionGrid,
)

	@info "$(modulelog()) - Saving $(npd.lname) data in the $(geo.name) GeoRegion for $(dt)"

	fnc = npdfnc(npd,geo,dt)
	if !isdir(dirname(fnc)); mkpath(dirname(fnc)) end
	if isfile(fnc)
		@info "$(modulelog()) - Overwrite stale NetCDF file $(fnc) ..."
        rm(fnc)
	end
	@info "$(modulelog()) - Creating NetCDF file $(fnc) ..."
	pds = makenpdnc(npd,fnc)

	pds.dim["longitude"] = length(ginfo.glon)
	pds.dim["latitude"]  = length(ginfo.glat)
	pds.dim["time"]      = size(var,3)

	nclon = defVar(pds,"longitude",Float32,("longitude",),attrib = Dict(
	    "units"     => "degrees_east",
	    "long_name" => "longitude",
	))

	nclat = defVar(pds,"latitude",Float32,("latitude",),attrib = Dict(
	    "units"     => "degrees_north",
	    "long_name" => "latitude",
	))

	ncvar = defVar(pds,"prcp_rate",Float32,("longitude","latitude","time"),attrib = Dict(
	    "units"     => "kg m**-2 s**-1",
	    "long_name" => "precipitation_rate",
		"full_name" => "Precipitation Rate",
	))

	nclon[:] = ginfo.glon
	nclat[:] = ginfo.glat
	ncvar[:] = var

	close(pds)

	@info "$(modulelog()) - $(npd.lname) data in the $(geo.name) GeoRegion for $(dt) has been saved into $(fnc)"

end

function makenpdnc(
	npd :: IMERGHalfHourly,
	fnc :: AbstractString
)

	return NCDataset(fnc,"c",attrib = Dict(
		"doi"              => npd.doi,
		"AlgorithmID"      => "3IMERGHH",
		"AlgorithmVersion" => "3IMERGH_6.3"
	))

end

function makenpdnc(
	npd :: IMERGDaily,
	fnc :: AbstractString
)

	return NCDataset(fnc,"c",attrib = Dict(
		"doi"         => npd.doi,
		"AlgorithmID" => "3IMERGD",
	))

end

function makenpdnc(
	npd :: IMERGMonthly,
	fnc :: AbstractString
)

	return NCDataset(fnc,"c",attrib = Dict(
		"doi"              => npd.doi,
		"AlgorithmID"      => "3IMERGM",
		"AlgorithmVersion" => "3IMERGM_6.3",
	))

end

function makenpdnc(
	npd :: TRMM3Hourly,
	fnc :: AbstractString
)

	return NCDataset(fnc,"c",attrib = Dict(
		"doi"              => npd.doi,
		"AlgorithmID"      => "3B42",
		"AlgorithmVersion" => "3B42_7.0"
	))

end

function makenpdnc(
	npd :: TRMMDaily,
	fnc :: AbstractString
)

	return NCDataset(fnc,"c",attrib = Dict("doi"=>npd.doi))

end

function makenpdnc(
	npd :: TRMMMonthly,
	fnc :: AbstractString
)

	return NCDataset(fnc,"c",attrib = Dict(
		"doi"              => npd.doi,
		"AlgorithmID"      => "3B43",
		"AlgorithmVersion" => "3B43_7.0",
	))

end

max_stars_repo_name: natgeo-wong/NASAPrecipitation

max_issues_repo_name: natgeo-wong/NASAPrecipitation

max_forks_repo_name: natgeo-wong/NASAPrecipitation


text:
#=
 Train a simple network to predict formation energy per atom (downloaded from Materials Project). =#
using CSV, DataFrames
using Random, Statistics
using Flux
using Flux: @epochs, glorot_uniform
using ChemistryFeaturization
using AtomicGraphNets
using Formatting
using Distributions
using CalibrationErrors
using CalibrationErrorsDistributions
using CalibrationTests

include("models.jl")
include("utils.jl")

function train_formation_energy(; num_epochs = 25, num_samples = 10, verbose = true)
    # model hyperparameters – keeping it pretty simple for now
    num_conv = 3 # how many convolutional layers?
    crys_fea_len = 32 # length of crystal feature vector after pooling (keep node dimension constant for now)
    num_hidden_layers = 1 # how many fully-connected layers after convolution and pooling?
    num_features, train_data, test_data = get_data()

    model = MC_CGCNN(
        num_features,
        num_conv = num_conv,
        atom_conv_feature_length = crys_fea_len,
        pooled_feature_length = (Int(crys_fea_len / 2)),
        num_hidden_layers = num_hidden_layers,
    )

    ps = Flux.params(model)
    opt = ADAM(0.001) # optimizer

    # define loss function and a callback to monitor progress
    loss(x, y) = Flux.Losses.mse(model(x), y)

    # train
    if verbose
        println("Training!")
    end

    function get_preds(data)
        predictions = []
        targets = []
        samples = []
        for (x, y) in data
            for sample in num_samples
                preds = model(x)
                push!(samples, preds[1])
            end
            samples = Flux.batch(samples)
            mean_predictions = mean(samples, dims = ndims(samples))
            mean_predictions = dropdims(mean_predictions, dims = ndims(samples))
            push!(predictions, mean_predictions[1])
            push!(targets, y)
        end
        return predictions, targets
    end

    # define kernel
    kernel = WassersteinExponentialKernel() ⊗ SqExponentialKernel()

    for epoch = 1:num_epochs
        train_loss = 0
        ntot = 0
        Flux.train!(loss, ps, train_data, opt)
        test_preds, test_targets = get_preds(test_data)
        train_preds, train_targets = get_preds(train_data)
        test_loss = Flux.Losses.mse(test_preds, test_targets) |> round4
        train_loss = Flux.Losses.mse(train_preds, train_targets) |> round4
        sigma = sqrt(train_loss)
        predictions = [Normal(mean, sigma) for mean in test_preds]

        # unbiased estimator of SKCE
        unbiased_estimator = UnbiasedSKCE(kernel)
        skce = calibrationerror(unbiased_estimator, predictions, test_targets) |> round4

        # biased estimator of SKCE
        biased_estimator = BiasedSKCE(kernel)
        biased_skce =
            calibrationerror(biased_estimator, predictions, test_targets) |> round4

        println(format("Train Loss: {}", train_loss))
        println(format("Test Loss: {}", test_loss))
        println(format("Unbiased SKCE: {}", skce))
        println(format("Biased SKCE: {}", biased_skce))
        println("===============")
    end

    return model
end

train_formation_energy()


text:
using SafeTestsets

@safetestset "Catch22" begin
using Catch22
import Catch22.featurenames, Catch22.testdata, Catch22.testoutput, Catch22.testnames
using Test

function isnearlyequalorallnan(a::AbstractArray, b::AbstractArray)
    replace!(a, NaN=>0.0)
    replace!(b, NaN=>0.0)
    all(map((x, y) -> isapprox(x, y, rtol=1e-5), a, b))
end
function isnearlyequalorallnan(a::Real, b::Real)
    isapprox(a, b, rtol=1e-5) || (isnan(a) && isnan(b))
end

# Test features one by one
println("Testing individual features")
@testset "Feature $(getname(f))" for f ∈ catch22
        @test isnearlyequalorallnan(f(testdata[:test]), testoutput[:test][f])
end;

# Test catch22, time series by time series
catch22(testdata[:test]) # To avoid compilation in test @time
println("Testing sample datasets")
function testFeatures(t::Symbol)
    @time f = catch22(testdata[t])
    ff = testoutput[t]
    isnearlyequalorallnan(Array(f), ff)
end
@testset "Dataset $f" for f in testnames
    @test testFeatures(f)
end;

# Test catch22 on a matrix
println("Testing 1000×100 array input")
catch22(randn(10, 10))
X = randn(1000, 100)
@testset "Matrices" begin
    @test @time catch22(X) isa FeatureMatrix
end;

println("Testing 1000×20×20 array input")
catch22(randn(10, 10, 10))
X = randn(1000, 20, 20)
@testset "Arrays" begin
    @test @time catch22(X) isa FeatureArray{T, 3} where {T}
end;

println("Testing FeatureSet operations")

@testset "FeatureSet" begin
    𝒇₁ = FeatureSet([sum, length], [:sum, :length], [["distribution"], ["sampling"]], ["∑x¹", "∑x⁰"])
    𝒇₂ = catch22[1:2]
    X = randn(100, 2)
    𝒇₃ = 𝒇₁ + 𝒇₂
    @test_nowarn 𝒇₁(X)
    @test_nowarn 𝒇₃(X)
    @test getnames(𝒇₃) == [:sum, :length , :DN_HistogramMode_5, :DN_HistogramMode_10]
    @test 𝒇₁ == 𝒇₃ \ 𝒇₂ == setdiff(𝒇₃, 𝒇₂)
    @test 𝒇₃ == 𝒇₁ ∪ 𝒇₂
    @test 𝒇₂ == 𝒇₃ ∩ 𝒇₂
end;

println("Testing CovarianceImage")
@testset "CovarianceImage" begin
    using Plots
    using Clustering
    X = hcat(randn(100, 100), 1:100)
    F = catch22(X)
    verbose = false
    @test covarianceimage(F; colormode=:top, verbose) isa Plots.Plot
    @test covarianceimage(F; colormode=:all, verbose) isa Plots.Plot
    @test covarianceimage(F; colormode=:raw, verbose, colorbargrad=:viridis) isa Plots.Plot
end

end


text:
"""
TeukolskyQNMFunctions.jl computes the quasinormal modes and eigenfunctions
for the spin s Teukolsky equation
using a horizon penetrating, hyperboloidally compactified coordinate system.
The main advantage of using these coordinates is that the quasinormal
wavefunctions are finite valued from the black hole to future null infinity.
"""
module TeukolskyQNMFunctions

export F, compute_om

include("CustomTypes.jl")
include("Spheroidal.jl")
include("RadialODE.jl")

using .CustomTypes

import .Spheroidal as SH
import .RadialODE  as RD

"""
    F(
      nr::myI,
      nl::myI,
      s::myI,
      l::myI,
      m::myI,
      a::myF,
      om::myC
    )::myF

Absolute difference of Lambda seperation constant
for radial and angular ODEs.
"""
function F(
      nr::myI,
      nl::myI,
      s::myI,
      l::myI,
      m::myI,
      a::myF,
      om::myC
   )::myF

   lmin = max(abs(s),abs(m))
   neig = l - lmin + 1 # number of eigenvalues

   la_s, _ = SH.eig_vals_vecs(nl, neig, s, m,  a*om)
   la_r, _ = RD.eig_vals_vecs_c(nr,     s, m, a, om)

   ## The Lambdas are ordered in size of smallest magnitude
   ## to largest magnitude, we ASSUME this is the same as the
   ## ordering of the l-angular indexing (this may not
   ## hold when a->1; need to check).

   return abs(la_s[l-lmin+1] - la_r)
end

"""
    compute_om(
      nr::myI,
      nl::myI,
      s::myI,
      l::myI,
      m::myI,
      a::myF,
      om::myC;
      tolerance::myF=tomyF(1e-6),
      epsilon::myF=tomyF(1e-6),
      gamma::myF=tomyF(1),
      verbose::Bool=false
   )::Tuple{
            myC,
            myC,
            Vector{myC},
            Vector{myC},
            Vector{myF}}

Search for quasinormal mode frequency in the complex plane
using Newton's method.

# Arguments

* `nr`: number of radial Chebyshev collocation points
* `nl`: number of spherical harmonic terms
* `s`:  spin of the field in the Teukolsky equation
* `l`:  l angular number
* `m`:  m angular number
* `a`:  black hole spin
* `om`: guess for the initial quasinormal mode
"""
function compute_om(
      nr::myI,
      nl::myI,
      s::myI,
      l::myI,
      m::myI,
      a::myF,
      om::myC;
      tolerance::myF=tomyF(1e-6),
      epsilon::myF=tomyF(1e-6),
      gamma::myF=tomyF(1),
      verbose::Bool=false
   )::Tuple{
            myC,
            myC,
            Vector{myC},
            Vector{myC},
            Vector{myF}}

   om_n   = tomyC(-1000)
   om_np1 = om

   f(omega) = F(nr,nl,s,l,m,a,omega)

   ## Newton search with 2nd order finite differences

   newgamma = gamma
   iterations = 1
   while abs(om_np1 - om_n) > tolerance

      ## if too many iterations, reduce search step size
      iterations += 1
      if iterations%100 == 0
         newgamma /= tomyF(2)
      end

      om_n = om_np1

      df = (
            (f(om_n+epsilon)    - f(om_n-epsilon)   )/(tomyF(2)*epsilon)
            +
            (f(om_n+im*epsilon) - f(om_n-im*epsilon))/(tomyF(2)*im*epsilon)
           )
      om_np1 = om_n - newgamma*f(om_n)/df
      if verbose==true
         println("om_np1=$om_np1\tdf=$df\tdiff=$(abs(om_np1-om_n))")
      end
   end

   lmin = max(abs(s),abs(m))
   neig = l - lmin + 1 # number of eigenvalues

   la_s, v_s        = SH.eig_vals_vecs(  nl, neig, s, m,  a*om_np1)
   la_r, v_r, rvals = RD.eig_vals_vecs_c(nr,       s, m, a, om_np1)

   return om_np1, la_r, v_s[:,l-lmin+1], v_r, rvals
end

end # module


text:
#=
This is Julia implementation of Gomoku Board by Tejan Karmali.
=#

import Base: ==, deepcopy, show
using IterTools: chain

# Cell occupants
#@enum CELL_OCC WHITE = -1 EMPTY BLACK  FILL KO UNKNOWN
WHITE, EMPTY, BLACK, FILL, KO, UNKNOWN = collect(-1:4)

# Represents "group not found" in the LibertyTracker object
MISSING_GROUP_ID = -1
#=
struct PlayerMove
  color
  move
end

# making a (x,y) coordinate for board from flattened index x
coordify(x; N = 19) = 1 + (x - 1) % N, 1 + (x - 1) ÷ N

get_first_n(x::Array{Int8, 3}, n) = size(x, 3) < n ? x : x[:, :, 1:n]
=#
mutable struct GomokuPosition <: Position
  #=
  board: an array
  n: an int representing moves played so far
  recent: a tuple of PlayerMoves, such that recent[end] is the last move.
  board_deltas: an array of shape (N, N, n) representing changes
      made to the board at each move (played move and captures).
      Should satisfy next_pos.board - next_pos.board_deltas[:, :, 0] == pos.board
  to_play: BLACK or WHITE
  =#
  env::GomokuEnv
  board::Array{Int8, 2}
  n::Int
  recent::Vector{PlayerMove}
  board_deltas::Array{Int8, 3}
  to_play::Int
  done::Bool
  winner::Int8
  function GomokuPosition(env::GomokuEnv; board = nothing, n = 0, recent = Vector{PlayerMove}(),
	board_deltas = nothing, to_play = BLACK)

    b = board != nothing ? board : deepcopy(env.EMPTY_BOARD)
    bd = board_deltas != nothing ? board_deltas : zeros(Int8, env.N, env.N, 0)
    done, winner = has_game_ended(b, env)

    new(env, b, n, recent, bd, to_play, done, winner)
  end
end

function deepcopy(pos::GomokuPosition)
  new_board = deepcopy(pos.board)
  new_recent = deepcopy(pos.recent)
  return GomokuPosition(pos.env; board = new_board, n = pos.n, recent = new_recent,
                  board_deltas = pos.board_deltas, to_play = pos.to_play)
end

function show(io::IO, pos::GomokuPosition)
  pretty_print_map = Dict{Int, String}([
            WHITE => "O",
            EMPTY => ".",
            BLACK => "X",
            FILL => "#",
        ])
  board = deepcopy(pos.board)
  N = size(board, 1)

  raw_board_contents = []
  for i = 1:N
    row = []
    for j = 1:N
      appended = length(pos.recent) != 0 && (i, j) == pos.recent[end].move ? "<" : " "
      push!(row, pretty_print_map[board[i,j]] * appended)
    end
    push!(raw_board_contents, join(row))
  end

  row_labels_fmt = [@sprintf("%2d ", i) for i = N:-1:1]
  row_labels = [string(i) for i = N:-1:1]
  annotated_board_contents = [join(r) for r in zip(row_labels_fmt, raw_board_contents, row_labels)]
  header_footer_rows = ["   " * join("ABCDEFGHJKLMNOPQRST"[1:N], " ") * "   "]
  annotated_board = join(collect(chain(header_footer_rows, annotated_board_contents,
                                        header_footer_rows)), "\n")
  details = "\nMove: $(pos.n)\n"
  turn = pos.done ? "GAME OVER\n" : "To Play: " * (pos.to_play == BLACK ? "X(BLACK)\n" : "O(WHITE)\n")
  println(annotated_board * details * turn)
end

is_move_legal(pos::GomokuPosition, move) = pos.board[move...] == EMPTY

all_legal_moves(pos::GomokuPosition) = vec(pos.board .== EMPTY)

has_game_ended(pos::GomokuPosition) = has_game_ended(pos.board, pos.env)

function has_game_ended(board::Array{Int8, 2}, env)
  n_in_row = env.n_in_row
  dim = env.N

  for h = 1:dim
    for w = 1:dim
      if w in 1:dim-n_in_row+1 && board[h,w]!=EMPTY
        pieces = Set([board[h, i] for i = w:w+n_in_row-1])
        length(pieces) == 1 && return true, collect(pieces)[1]
      end

      if h in 1:dim-n_in_row+1 && board[h,w]!=EMPTY
        pieces = Set([board[i, w] for i = h:h+n_in_row-1])
	length(pieces) == 1 && return true, collect(pieces)[1]
      end

      if w in 1:dim-n_in_row+1 && h in 1:dim-n_in_row+1 && board[h,w]!=EMPTY
        pieces = Set([board[h+i,w+i] for i=0:n_in_row-1])
        length(pieces) == 1 && return true, collect(pieces)[1]
      end

      if w in n_in_row:dim && h in 1:dim-n_in_row+1 && board[h,w]!=EMPTY
        pieces = Set([board[h+i,w-i] for i=0:n_in_row-1])
        length(pieces) == 1 && return true, collect(pieces)[1]
      end
    end
  end

  sum(board .== EMPTY) == 0 && return true, EMPTY

  false, EMPTY
end

function flip_playerturn!(pos::GomokuPosition; mutate = false)
  new_pos = mutate ? pos : deepcopy(pos)
  new_pos.to_play *= -1
  return new_pos
end

function play_move!(pos::GomokuPosition, c; color = nothing, mutate = false)
  if color == nothing
    color = pos.to_play
  end

  new_pos = mutate ? pos : deepcopy(pos)

  @assert !new_pos.done

  if !is_move_legal(pos, c)
    throw(IllegalMove())
  end

  new_pos.board[c...] = color
  opp_color = -color

  N = size(pos.board, 1)

  new_board_delta = zeros(Int8, N, N)
  new_board_delta[c...] = color

  new_pos.n += 1
  push!(new_pos.recent, PlayerMove(color, c))
  # keep a rolling history of last 7 deltas - that's all we'll need to
  # extract the last 8 board states.

  new_pos.board_deltas = cat(dims=3, reshape(new_board_delta, N, N, 1),
			     get_first_n(new_pos.board_deltas, new_pos.env.planes-2))
  new_pos.to_play *= -1
  new_pos.done, new_pos.winner = has_game_ended(new_pos)
  return new_pos
end

score(pos::GomokuPosition) = pos.winner

function result(pos::GomokuPosition)
  points = score(pos)
  if points > 0
    return 1
  elseif points < 0
    return -1
  else
    return 0
  end
end

function result_string(pos::GomokuPosition)
  points = score(pos)
  if points > 0
    return "B"
  elseif points < 0
    return "W"
  else
    return "DRAW"
  end
end

function replay_position(pos::GomokuPosition, result)
  #=
  Wrapper for a go.GoPosition which replays its history.
  Assumes an empty start position! (i.e. no handicap, and history must be exhaustive.)

  Result must be passed in, since a resign cannot be inferred from position
  history alone.

  for position_w_context in replay_position(position):
    print(position_w_context.position)
  =#
  pos.n == length(pos.recent) ? nothing : throw(AssertionError("GomokuPosition history is incomplete"))
  replay_buffer = Vector{PositionWithContext}()
  dummy_pos = GomokuPosition(pos.env)

  for player_move in pos.recent
    color, next_move = player_move.color, player_move.move
    push!(replay_buffer, PositionWithContext(dummy_pos, next_move, result))
    dummy_pos = play_move!(dummy_pos, next_move, color = color)
  end
  return replay_buffer
end


text:
using IncrementalInference
using Test

## a new factor that is broken
struct BrokenFactor{T<: SamplableBelief} <: AbstractRelativeRoots
    Z::T
end

function IIF.getSample(cf::CalcFactor{<:BrokenFactor})
  return rand(cf.factor.Z, 1)
end

function (s::CalcFactor{<:BrokenFactor})(z,
                                        wxi,
                                        wxj)
    #
    error("User factor has a bug.")
end

# # FIXME consolidate with CalcFactor according to #467
# function (s::BrokenFactor{<:IIF.ParametricTypes})(X1::AbstractArray{<:Real},
#                                                   X2::AbstractArray{<:Real};
#                                                   userdata::Union{Nothing,FactorMetadata}=nothing )
#     error("User factor has a bug -- USE NEW CalcFactor API INSTEAD, v0.21.")
# end

##

@testset "Test CSM monitor/watchdog on errors" begin

##

# create a factor graph
fg = generateGraph_LineStep(10;
                            poseEvery=1,
                            landmarkEvery=5,
                            posePriorsAt=[0],
                            sightDistance=4,
                            solverParams=SolverParams(algorithms=[:default, :parametric]))
#

initAll!(fg)

##

#TODO test FSM watchdog
# add a broken factor - mid
addFactor!(fg, [:x9, :lm10], BrokenFactor(Normal()); graphinit=false)
smtasks = Task[]
@test_throws CompositeException tree = IIF.solveTree!(fg; smtasks=smtasks);
sleep(0.1)

## Test parametric solve also

addFactor!(fg, [:x9, :lm10], BrokenFactor(Normal()); graphinit=false)
##
# IIF.solveTree!(fg; smtasks=smtasks, algorithm = :parametric)
##
@test_throws CompositeException tree2 = IIF.solveTree!(fg; smtasks=smtasks, algorithm = :parametric)
sleep(0.1)

deleteFactor!(fg, :x9lm10f2)

## add a broken factor - leave

addFactor!(fg, [:x10, :lm10], BrokenFactor(Normal()); graphinit=false)
@test_throws CompositeException tree2 = IIF.solveTree!(fg; smtasks=smtasks, algorithm = :parametric)
sleep(0.1)

deleteFactor!(fg, :x10lm10f2)

## add a broken factor - root
addFactor!(fg, [:x7, :lm10], BrokenFactor(Normal()); graphinit=false)
@test_throws CompositeException tree2 = IIF.solveTree!(fg; smtasks=smtasks, algorithm = :parametric)

##

end

@testset "test CSM debug options work" begin

## create a factor graph

fg = generateGraph_LineStep(10;
                            poseEvery=1,
                            landmarkEvery=5,
                            posePriorsAt=[0],
                            sightDistance=4 )
#

smtasks = Task[];
solveTree!(fg, smtasks=smtasks, recordcliqs=ls(fg));

## make sure we fetch the results

hists = fetchCliqHistoryAll!(smtasks);

printCSMHistoryLogical(hists)

printCSMHistorySequential(hists)

## test CSM resolve steps

@info "test repeatCSMStep"
csmc_ = repeatCSMStep!(hists, 1, 1);

##

end


text:
using SimpleGradientDescent
using Base.Test

@testset "naive solver" begin
    num_vars = 10

    # Let's define a very simple cost function:
    xstar = [1., 2, 3, 4, 5, 6, 7, 8, 9, 10]
    cost(x) = sum((x .- xstar) .^ 4)

    # getmodel() returns a MathProgBase model which will allow us to perform
    # the gradient descent. We can specify a particular solver if we want to,
    # or it will try to use MathProgBase's default NLP solver if we don't.
    model = getmodel(cost, num_vars, NaiveSolver(num_vars))

    # Set the initial point from which the solver will start:
    setwarmstart!(model, [0.0 for i in 1:num_vars])

    # Optimize and extract the solution:
    # (note that setwarmstart!, optimize!, status, and getsolution are all
    # just re-exported from MathProgBase for convenience)
    optimize!(model)
    @test status(model) == :Optimal
    x = getsolution(model)
    @show x xstar
    @test all(isapprox(x, xstar, atol=1e-3))

    # We can specify lower and upper bounds on the decision variables:
    model = getmodel(cost, num_vars, NaiveSolver(num_vars),
                    [-1.0 for i in 1:num_vars],
                    [1.0 for i in 1:num_vars])
    setwarmstart!(model, [0.0 for i in 1:num_vars])
    optimize!(model)
    x = getsolution(model)
    @show x xstar
    # Since we've constrained x to lie in [0, 1], it won't reach xstar
    @test all(isapprox(x, [1.0 for i in 1:num_vars], atol=1e-3))
end

@testset "naive preconditioned" begin
    num_vars = 10
    xstar = [1., 2, 3, 4, 5, 6, 7, 8, 9, 10]
    cost(x) = sum((x .- xstar) .^ 4)
    unconditioned_model = getmodel(cost, num_vars, NaiveSolver(num_vars; precondition_divisors=ones(num_vars),
        iteration_limit=2))
    setwarmstart!(unconditioned_model, [0.0 for i in 1:num_vars])
    optimize!(unconditioned_model)
    @test status(unconditioned_model) == :UserLimit

    conditioned_model = getmodel(cost, num_vars, NaiveSolver(num_vars; precondition_divisors=collect(1:num_vars),
        iteration_limit=2))
    setwarmstart!(conditioned_model, [0.0 for i in 1:num_vars])
    optimize!(conditioned_model)
    @test status(conditioned_model) == :Optimal
    x = getsolution(conditioned_model)
    @test all(isapprox(x, xstar, atol=1e-3))

end


text:
include("../src/RandomField.jl")
using .RandomField

# choose parameters:
## remember, we need that 4|n_i for i=1,...,d
N = [400,400]		# number pixels per dimension
L = [20.0,80.0]		# real space box is [-l,l]^d
f(V::Vector{T} where T<:Real) = gaus_rf(V)

# derived parameters
d = length(N) 		# dimension
H = (pi/2) .* N ./L

# allocate memory
V = Cloud(d,N)
X = Cloud(d,N)
F = Cscalar(V)
G = Cscalar(X)
G_tru = Cscalar(X)

# initialize grids
init!(X,L)
init!(V,H)

sample!(F,f)

ift!(G,F)

sample!(G_tru,gaus_rf)

tocsv(slice(G,2), "../data/x.g2.csv")
tocsv(slice(G_tru,2), "../data/x.gtrue2.csv")


text:
export size_criterion, sqrt_criterion, change_criterion, log_criterion, epsilon_criterion

"""
Stops when the distance between far items achieves the given `e`
"""
function epsilon_criterion(e)
   (dmaxlist, database) -> dmaxlist[end] < e
end

"""
Stops when the number of far items are equal or larger than the given `maxsize`
"""
function size_criterion(maxsize)
   (dmaxlist, database) -> length(dmaxlist) >= maxsize
end

"""
Stops when the number of far items are equal or larger than the square root of the size of the database
"""
function sqrt_criterion()
   (dmaxlist, database) -> length(dmaxlist) >= Int(length(database) |> sqrt |> round)
end

"""
Stops when the number of far items are equal or larger than logarithm-2 of the size of the database
"""
function log_criterion()
   (dmaxlist, database) -> length(dmaxlist) >= Int(length(database) |> log2 |> round)
end

"""
Stops the process whenever the maximum distance converges, i.e., after `window` far items the maximum distance
change is below or equal to the allowed tolerance `tol`
"""
function change_criterion(tol=0.001, window=3)
    mlist = Float64[]
    count = 0.0
    function stop(dmaxlist, database)
        count += dmaxlist[end]

        if length(dmaxlist) % window != 1
            return false
        end
        push!(mlist, count)
        count = 0.0
        if length(dmaxlist) < 2
            return false
        end

        s = abs(mlist[end] - mlist[end-1])
        return s <= tol
    end

    return stop
end

"""
It nevers stops by side, it explores the entire dataset making a full farthest first traversal
"""
function salesman_criterion()
    function stop(dmaxlist, dataset)
        return false
    end
end


text:
function run!{T,NI}(method::Union{ODENEB, StaticNEB}, E, dE, x0::Path{T,NI})
   # read all the parameters
   @unpack k, interp, tol, maxtol, maxnit, precon_scheme, path_traverse, fixed_ends,
            verbose = method
   @unpack precon, precon_prep! = precon_scheme
   @unpack direction = path_traverse
   # initialise variables
   x = x0.x

   nit = 0
   numdE, numE = 0, 0
   log = PathLog()
   # and just start looping
   if verbose >= 2
       @printf("SADDLESEARCH:         k  =  %1.2e        <- parameters\n", k)
   end
   file = []
   if verbose >= 4
       dt = Dates.format(now(), "d-m-yyyy_HH:MM")
       file = open("log_$(dt).txt", "w")
       strlog = @sprintf("SADDLESEARCH:         k  =  %1.2e        <- parameters\n", k)
       write(file, strlog)
       flush(file)
   end
   xout, log, alpha = odesolve(solver(method),
               (X, P, nit) -> forces(P, typeof(x0), X, dE, precon_scheme,
                                       direction(NI, nit), k, interp, fixed_ends),
               vec(x), log; file = file,
               tol = tol, maxtol = maxtol, maxnit = maxnit,
               P = precon,
               precon_prep! = (P, X) -> precon_prep!(P, convert(typeof(x0), X)),
               method = "$(typeof(method))",
               verbose = verbose)

   x_return = verbose < 4 ? convert(typeof(x0), xout[end]) : [convert(typeof(x0), xout_n) for xout_n in xout]
   return x_return, log, alpha
end

function run!{T,NI}(method::AccelNEB, E, dE, ddE, x0::Path{T,NI})
   # read all the parameters
   @unpack k, interp, tol, maxtol, maxnit, precon_scheme, path_traverse, fixed_ends,
            verbose = method
   @unpack precon, precon_prep! = precon_scheme
   @unpack direction = path_traverse
   # initialise variables
   x = x0.x

   nit = 0
   numdE, numE = 0, 0
   log = PathLog()
   # and just start looping
   if verbose >= 2
       @printf("SADDLESEARCH:         k  =  %1.2e        <- parameters\n", k)
   end
   file = []
   if verbose >= 4
       dt = Dates.format(now(), "d-m-yyyy_HH:MM")
       file = open("log_$(dt).txt", "w")
       strlog = @sprintf("SADDLESEARCH:         k  =  %1.2e        <- parameters\n", k)
       write(file, strlog)
       flush(file)
   end

   xout, log, alpha = odesolve(solver(method),
               (X, P, nit) -> forces(P, typeof(x0), X, dE, precon_scheme,
                                       direction(NI, nit), k, interp, fixed_ends),
               (X, P) -> jacobian(P, typeof(x0), X, dE, ddE, k),
               vec(x), log; file = file,
               tol = tol, maxtol=maxtol, maxnit=maxnit,
               P = precon,
               precon_prep! = (P, X) -> precon_prep!(P, convert(typeof(x0), X)),
               method = "$(typeof(method))",
               verbose = verbose)

   x_return = verbose < 4 ? convert(typeof(x0), xout[end]) : [convert(typeof(x0), xout_n) for xout_n in xout]
   return x_return, log, alpha
end

# forcing term for NEB method
function forces{T,NI}(precon, path_type::Type{Path{T,NI}}, X::Vector{Float64}, dE, precon_scheme,
                  direction, k::Float64, interp::Int, fixed_ends::Bool)

   x = convert(path_type, X)
   dxds = deepcopy(x)

   # preconditioner
   Np = size(precon, 1); N = length(x)
   P(i) = precon[mod(i-1,Np)+1, 1]
   P(i, j) = precon[mod(i-1,Np)+1, mod(j-1,Np)+1]

   # interpolate path to find tangents and 2nd derivatives
   if interp == 1
      # central finite differences
      dxds = [[zeros(x[1])]; [0.5*(x[i+1]-x[i-1]) for i=2:N-1]; [zeros(x[1])]]
      d²xds² = [[zeros(x[1])]; [x[i+1] - 2*x[i] + x[i-1] for i=2:N-1];
               [zeros(x[1])]]
   elseif interp > 1
      # splines
      ds = [dist(precon_scheme, P, x, i) for i=1:length(x)-1]
      param = [0; [sum(ds[1:i]) for i in 1:length(ds)]]
      param /= param[end]; param[end] = 1.
      d²xds² = parametrise!(dxds, x, ds, parametrisation = param)
      k *= (1/(N*N))
   else
      error("SADDLESEARCH: invalid `interpolate` parameter")
   end
   dxds ./= point_norm(precon_scheme, P, dxds)
   dxds[1] = zeros(dxds[1]); dxds[end] = zeros(dxds[1])

   # elastic interactions between adjacent images
   Fk = elastic_force(precon_scheme, P, k*N*N, dxds, d²xds²)

   # potential gradient
   dE0_temp = []
   if !fixed_ends
      dE0_temp = [dE(x[i]) for i in direction]
      cost = N
   else
      dE0_temp = [[zeros(x[1])]; [dE(x[i]) for i in direction[2:end-1]];
                  [zeros(x[1])]]
      cost = N - 2
   end
   dE0 = [dE0_temp[i] for i in direction]

   # projecting out tangent term of potential gradient
   dE0⟂ = proj_grad(precon_scheme, P, dE0, dxds)

   # collecting force term
   F = forcing(precon_scheme, precon, dE0⟂-Fk)

   # residual error
   res = maxres(precon_scheme, P, dE0⟂)

   return F, res, cost, (X, Y) -> dot_P(precon_scheme, convert(path_type, X), P, convert(path_type, Y))
end

function jacobian{T,NI}(precon, path_type::Type{Path{T,NI}}, X::Vector{Float64},
   dE, ddE, k::Float64)
   x = convert(path_type, X)

   # preconditioner
   Np = size(precon, 1)
   P(i) = precon[mod(i-1,Np)+1, 1]
   P(i, j) = precon[mod(i-1,Np)+1, mod(j-1,Np)+1]

   hessian = ddE.precon; hessian_prep! = ddE.precon_prep!
   hessian = hessian_prep!(hessian, x)
   H(i) = hessian[mod(i-1,Np)+1, 1]
   H(i, j) = hessian[mod(i-1,Np)+1, mod(j-1,Np)+1]

   N = length(x); M = length(x[1])
   O = zeros(M, M); J = fill(O,(N, N))
   [J[n,n-1] = ∂Fⁿ⁻(x, n, dE, P) + ∂Sⁿ⁻(k, x, n, P) for n=2:N-1]
   if Np==1 && P(1)==I
      [J[n,n] = δFⁿ(x, n, H, P, H) + ∂Sⁿ(k, x, n, P) for n=1:N]
   else
      [J[n,n] = δFⁿ(x, n, H, P, n -> I) + ∂Sⁿ(k, x, n, P) for n=1:N]
   end
   [J[n,n+1] = ∂Fⁿ⁺(x, n, dE, P) + ∂Sⁿ⁺(k, x, n, P) for n=2:N-1]

   return ref(J)
end


text:
#
# Exponential fit
#

struct SingleExponential
  A :: Float64
  b :: Float64
  C :: Float64
  R :: Float64
  x :: Vector{Float64}
  y :: Vector{Float64}
  ypred :: Vector{Float64}
  residues :: Vector{Float64}
end

struct MultipleExponential
  n :: Int64
  A :: Vector{Float64}
  b :: Vector{Float64}
  C :: Float64
  R :: Float64
  x :: Vector{Float64}
  y :: Vector{Float64}
  ypred :: Vector{Float64}
  residues :: Vector{Float64}
end

function sum_of_exps(x :: Numbers, p :: AbstractVector)
  n = round(Int64,(length(p)-1)/2)
  f = p[length(p)]
  for i in 1:n
    f = f + p[i]*exp(x/p[n+i])
  end
  return f
end

function exp_model(X :: AbstractVector, p :: AbstractVector)
  f = Vector{Float64}(undef,length(X))
  for i in 1:length(X)
    f[i] = sum_of_exps(X[i],p)
  end
  return f
end

function fitexponential(X :: AbstractVector, Y :: AbstractVector, options :: Options; n :: Int = 1)
  X, Y = checkdata(X,Y)
  model(x,p) = exp_model(x,p)
  fit = find_best_fit(model,X,Y,2*n+1,options)
  R = pearson(X,Y,model,fit)
  x, y, ypred = finexy(X,options.fine,model,fit)
  if n == 1
    return SingleExponential(fit.param[1],fit.param[2],fit.param[3],
                             R,x,y,ypred,fit.resid)
  else
    ind = collect(1:n)
    sort!(ind,by=i->fit.param[n+i])
    A = fit.param[1:n][ind]
    b = fit.param[n+1:2*n][ind]
    return MultipleExponential(n,A,b,fit.param[2*n+1],
                               R,x,y,ypred,fit.resid)
  end
end
fitexponential(X :: AbstractVector, Y :: AbstractVector; n :: Int = 1) = fitexponential(X,Y,Options(),n=n)
fitexp(X :: AbstractVector, Y :: AbstractVector; n :: Int = 1) = fitexponential(X,Y,Options(),n=n)
fitexp(X :: AbstractVector, Y :: AbstractVector, options :: Options; n :: Int = 1) = fitexponential(X,Y,options,n=n)

function Base.show( io :: IO, fit :: SingleExponential )
  println("")
  println(" ------------ Single Exponential fit ----------- ")
  println("")
  println(" Equation: y = A exp(x/b) + C")
  println("")
  println(" With: A = ", fit.A)
  println("       b = ", fit.b)
  println("       C = ", fit.C)
  println("")
  println(" Pearson correlation coefficient, R = ", fit.R)
  println("")
  println(" Predicted Y: ypred = [",fit.ypred[1],", ",fit.ypred[2],"...")
  println(" residues = [", fit.residues[1],", ",fit.residues[2],"...")
  println("")
  println(" ----------------------------------------------- ")
end

function Base.show( io :: IO, fit :: MultipleExponential )
  println("")
  println(" -------- Multiple-exponential fit ------------- ")
  println("")
  println(" Equation: y = sum(A[i] exp(x/b[i]) for i in 1:$(fit.n)) + C ")
  println("")
  println(" With: A = ", fit.A)
  println("       b = ", fit.b)
  println("       C = ", fit.C)
  println("")
  println(" Pearson correlation coefficient, R = ", fit.R)
  println("")
  println(" Predicted Y: ypred = [",fit.ypred[1],", ",fit.ypred[2],"...")
  println(" residues = [", fit.residues[1],", ",fit.residues[2],"...")
  println("")
  println(" ----------------------------------------------- ")
end

export fitexp, fitexponential


text:
struct SimulationInputs
    dynamic_injectors::Vector{DynamicWrapper{<:PSY.DynamicInjection}}
    static_injectors::Vector
    static_loads::Vector
    dynamic_branches::Vector{BranchWrapper}
    injection_n_states::Int
    branches_n_states::Int
    variable_count::Int
    inner_vars_count::Int
    bus_count::Int
    ode_range::UnitRange{Int}
    ybus_rectangular::SparseArrays.SparseMatrixCSC{Float64, Int}
    dyn_lines::Bool
    total_shunts::SparseArrays.SparseMatrixCSC{Float64, Int}
    lookup::Dict{Int, Int}
    DAE_vector::Vector{Bool}
    mass_matrix::LinearAlgebra.Diagonal{Float64}
    global_vars_update_pointers::Dict{Int, Int}
    global_state_map::MAPPING_DICT

    function SimulationInputs(
        sys::PSY.System,
        ::Type{T},
    ) where {T <: Union{ConstantFrequency, ReferenceBus}}
        n_buses = get_n_buses(sys)
        Ybus, lookup = _get_ybus(sys)

        TimerOutputs.@timeit BUILD_TIMER "Wrap Branches" begin
            wrapped_branches = _wrap_dynamic_branches(sys, lookup)
            has_dyn_lines = !isempty(wrapped_branches)
            branch_state_counts = 2 * length(wrapped_branches)
            injection_start = 2 * n_buses + branch_state_counts + 1
        end

        TimerOutputs.@timeit BUILD_TIMER "Wrap Dynamic Injectors" begin
            wrapped_injectors = _wrap_dynamic_injector_data(sys, lookup, injection_start)
            var_count = wrapped_injectors[end].ix_range[end]
        end

        TimerOutputs.@timeit BUILD_TIMER "Calculate MM, DAE_vector, Total Shunts" begin
            mass_matrix = _make_mass_matrix(wrapped_injectors, var_count, n_buses)
            DAE_vector = _make_DAE_vector(mass_matrix, var_count, n_buses)
            total_shunts = _make_total_shunts(wrapped_branches, n_buses)
        end

        TimerOutputs.@timeit BUILD_TIMER "Wrap Static Injectors" begin
            wrapped_loads = _wrap_loads(sys, lookup)
            wrapped_static_injectors = _wrap_static_injectors(sys, lookup)
        end

        _adjust_states!(
            DAE_vector,
            mass_matrix,
            total_shunts,
            n_buses,
            PSY.get_frequency(sys),
        )

        global_vars =
            _make_global_variable_index(wrapped_injectors, wrapped_static_injectors, T)

        inner_vars_count = 0
        for i in length(wrapped_injectors):-1:1
            if length(wrapped_injectors[i].inner_vars_index) > 0
                inner_vars_count = wrapped_injectors[i].inner_vars_index[end]
                break
            end
        end
        IS.@assert_op inner_vars_count > 0
        new(
            wrapped_injectors,
            wrapped_static_injectors,
            wrapped_loads,
            wrapped_branches,
            var_count - 2 * n_buses - branch_state_counts,
            branch_state_counts,
            var_count,
            inner_vars_count,
            n_buses,
            injection_start:var_count,
            Ybus,
            has_dyn_lines,
            total_shunts,
            lookup,
            DAE_vector,
            mass_matrix,
            global_vars,
            MAPPING_DICT(),
        )
    end
end

get_dynamic_injectors(inputs::SimulationInputs) = inputs.dynamic_injectors
get_dynamic_branches(inputs::SimulationInputs) = inputs.dynamic_branches
get_static_injectors(inputs::SimulationInputs) = inputs.static_injectors
get_static_loads(inputs::SimulationInputs) = inputs.static_loads
get_ybus(inputs::SimulationInputs) = inputs.ybus_rectangular
get_total_shunts(inputs::SimulationInputs) = inputs.total_shunts
get_lookup(inputs::SimulationInputs) = inputs.lookup
get_DAE_vector(inputs::SimulationInputs) = inputs.DAE_vector
get_mass_matrix(inputs::SimulationInputs) = inputs.mass_matrix
has_dyn_lines(inputs::SimulationInputs) = inputs.dyn_lines
get_global_vars_update_pointers(inputs::SimulationInputs) =
    inputs.global_vars_update_pointers

get_injection_n_states(inputs::SimulationInputs) = inputs.injection_n_states
get_branches_n_states(inputs::SimulationInputs) = inputs.branches_n_states
get_variable_count(inputs::SimulationInputs) = inputs.variable_count
get_inner_vars_count(inputs::SimulationInputs) = inputs.inner_vars_count
get_ode_ouput_range(inputs::SimulationInputs) = inputs.ode_range
get_bus_count(inputs::SimulationInputs) = inputs.bus_count
get_bus_range(inputs::SimulationInputs) = 1:(2 * inputs.bus_count)

# Utility function not to be used for performance sensitive operations
function get_voltage_buses_ix(inputs::SimulationInputs)
    n_buses = get_bus_count(inputs::SimulationInputs)
    return findall(get_DAE_vector(inputs)[1:n_buses])
end

# Utility function not to be used for performance sensitive operations
function get_current_buses_ix(inputs::SimulationInputs)
    n_buses = get_bus_count(inputs::SimulationInputs)
    return findall(.!get_DAE_vector(inputs)[1:n_buses])
end

"""
SimulationInputs build function for MassMatrixModels
"""
function SimulationInputs(::Type{MassMatrixModel}, sys::PSY.System, frequency_reference)
    return SimulationInputs(sys, frequency_reference)
end

"""
SimulationInputs build function for ResidualModels
"""
function SimulationInputs(::Type{ResidualModel}, sys::PSY.System, frequency_reference)
    return SimulationInputs(sys, frequency_reference)
end

function _wrap_dynamic_injector_data(sys::PSY.System, lookup, injection_start::Int)
    injector_data = get_injectors_with_dynamics(sys)
    isempty(injector_data) && error("System doesn't contain any DynamicInjection devices")
    # TODO: Needs a better container that isn't parametrized on an abstract type
    wrapped_injector = Vector(undef, length(injector_data))
    injection_count = 1
    inner_vars_count = 1
    sys_base_power = PSY.get_base_power(sys)
    sys_base_freq = PSY.get_frequency(sys)
    @assert !isempty(injector_data)
    for (ix, device) in enumerate(injector_data)
        @debug "Wrapping $(PSY.get_name(device))"
        dynamic_device = PSY.get_dynamic_injector(device)
        n_states = PSY.get_n_states(dynamic_device)
        ix_range = range(injection_start, length = n_states)
        ode_range = range(injection_count, length = n_states)
        bus_n = PSY.get_number(PSY.get_bus(device))
        bus_ix = lookup[bus_n]
        inner_vars_range =
            range(inner_vars_count, length = get_inner_vars_count(dynamic_device))
        @debug "ix_range=$ix_range ode_range=$ode_range inner_vars_range= $inner_vars_range"
        wrapped_injector[ix] = DynamicWrapper(
            device,
            bus_ix,
            ix_range,
            ode_range,
            inner_vars_range,
            sys_base_power,
            sys_base_freq,
        )
        injection_count += n_states
        injection_start += n_states
        inner_vars_count =
            length(inner_vars_range) > 0 ? inner_vars_range[end] : inner_vars_count
    end
    return wrapped_injector
end

function _wrap_dynamic_branches(sys::PSY.System, lookup::Dict{Int, Int})
    branches_start = 2 * get_n_buses(sys) + 1
    sys_base_power = PSY.get_base_power(sys)
    sys_base_freq = PSY.get_frequency(sys)
    dynamic_branches = get_dynamic_branches(sys)
    wrapped_branches = Vector{BranchWrapper}(undef, length(dynamic_branches))
    if !isempty(wrapped_branches)
        branches_count = 1
        for (ix, br) in enumerate(dynamic_branches)
            @debug "Wrapping Branch $(PSY.get_name(br))"
            arc = PSY.get_arc(br)
            n_states = PSY.get_n_states(br)
            from_bus_number = PSY.get_number(arc.from)
            to_bus_number = PSY.get_number(arc.to)
            bus_ix_from = lookup[from_bus_number]
            bus_ix_to = lookup[to_bus_number]
            ix_range = range(branches_start, length = n_states)
            ode_range = range(branches_count, length = n_states)
            @debug "ix_range=$ix_range ode_range=$ode_range"
            wrapped_branches[ix] = BranchWrapper(
                br,
                bus_ix_from,
                bus_ix_to,
                ix_range,
                ode_range,
                sys_base_power,
                sys_base_freq,
            )
            branches_count += n_states
            branches_start += n_states
        end
    else
        @debug("System doesn't contain Dynamic Branches")
    end
    return wrapped_branches
end

function _wrap_static_injectors(sys::PSY.System, lookup::Dict{Int, Int})
    static_injection_data = get_injection_without_dynamics(sys)
    container = Vector{StaticWrapper}(undef, length(static_injection_data))
    for (ix, ld) in enumerate(static_injection_data)
        bus_n = PSY.get_number(PSY.get_bus(ld))
        bus_ix = lookup[bus_n]
        container[ix] = StaticWrapper(ld, bus_ix)
    end
    return container
end

function _wrap_loads(sys::PSY.System, lookup::Dict{Int, Int})
    # This needs to change if we implement dynamic load models
    static_loads =
        PSY.get_components(PSY.ElectricLoad, sys, x -> !isa(x, PSY.FixedAdmittance))
    container = Vector(undef, length(static_loads))
    for (ix, ld) in enumerate(static_loads)
        bus_n = PSY.get_number(PSY.get_bus(ld))
        bus_ix = lookup[bus_n]
        container[ix] = StaticWrapper(ld, bus_ix)
    end
    return container
end

function _get_ybus(sys::PSY.System)
    n_buses = length(PSY.get_components(PSY.Bus, sys))
    dyn_lines = PSY.get_components(PSY.DynamicBranch, sys)
    if !isempty(PSY.get_components(PSY.ACBranch, sys))
        Ybus_ = PSY.Ybus(sys)
        ybus = Ybus_[:, :]
        lookup = Ybus_.lookup[1]
        ybus_rectangular = transform_ybus_to_rectangular(ybus)
        for br in dyn_lines
            ybus_update!(ybus_rectangular, br, lookup, -1.0)
        end
    else
        ybus_rectangular =
            SparseArrays.SparseMatrixCSC{Float64, Int}(zeros(2 * n_buses, 2 * n_buses))
        lookup = Dict{Int.Int}()
    end
    return ybus_rectangular, lookup
end

function _static_injection_inputs!(inputs::SimulationInputs, ::Vector{Int}, sys::PSY.System)
    for s in PSY.get_components(PSY.StaticInjection, sys)
        # Temporary condition needs to be removed for other static injections
        PSY.get_dynamic_injector(s) !== nothing && continue
        index_static_injection(inputs, s)
        _attach_control_refs!(s)
    end
    return
end

function _make_mass_matrix(wrapped_injectors, var_count::Int, n_buses::Int)
    mass_matrix = LinearAlgebra.Diagonal(ones(var_count))
    mass_matrix[1:(2 * n_buses), 1:(2 * n_buses)] .= 0.0
    for d in wrapped_injectors
        device_mass_matrix_entries!(mass_matrix, d)
    end
    return mass_matrix
end

function _init_DAE_vector(var_count::Int, n_buses::Int)
    DAE_vector = trues(var_count)
    DAE_vector[1:(n_buses * 2)] .= false
    return DAE_vector
end

function _make_DAE_vector(mass_matrix::AbstractArray, var_count::Int, n_buses::Int)
    DAE_vector = _init_DAE_vector(var_count, n_buses)
    bus_range = 1:(2 * n_buses)
    ode_range = (2 * n_buses + 1):var_count
    mass_buses = @view mass_matrix[bus_range, bus_range]
    DAE_ode = @view DAE_vector[ode_range]
    mass_ode = @view mass_matrix[ode_range, ode_range]
    for i in eachindex(DAE_vector[bus_range])
        IS.@assert_op DAE_vector[bus_range][i] == (mass_buses[i, i] > 0.0)
    end
    for i in eachindex(DAE_ode)
        DAE_ode[i] = (mass_ode[i, i] > 0.0)
    end
    return DAE_vector
end

function _make_total_shunts(wrapped_branches, n_buses::Int)
    shunts = SparseArrays.SparseMatrixCSC{Float64, Int}(zeros(2 * n_buses, 2 * n_buses))
    for br in wrapped_branches
        bus_ix_from = get_bus_ix_from(br)
        bus_ix_to = get_bus_ix_to(br)
        b_from = PSY.get_b(br).from
        b_to = PSY.get_b(br).to
        shunts[bus_ix_from, bus_ix_from + n_buses] += b_from
        shunts[bus_ix_to, bus_ix_to + n_buses] += b_to
        shunts[bus_ix_from + n_buses, bus_ix_from] -= b_from
        shunts[bus_ix_to + n_buses, bus_ix_to] -= b_to
    end
    return shunts
end

function _adjust_states!(
    DAE_vector::BitVector,
    mass_matrix::LinearAlgebra.Diagonal{Float64},
    total_shunts::SparseArrays.SparseMatrixCSC{Float64, Int},
    n_buses::Int,
    sys_f::Float64,
)
    all(iszero.(total_shunts)) && return
    line_constant = 1 / (2.0 * π * sys_f)
    # Takes the upper quadrant of the rectangular shunts matrix
    shunts = LinearAlgebra.diag(total_shunts[1:n_buses, (n_buses + 1):end])
    for (ix, val) in enumerate(shunts)
        if val > 0
            @debug "Found shunt with value $val in bus index $ix"
            mass_matrix[ix, ix] =
                mass_matrix[ix + n_buses, ix + n_buses] = val * line_constant
            DAE_vector[ix] = DAE_vector[ix + n_buses] = true
        end
    end
    return
end

function _make_global_variable_index(
    wrapped_injectors::Vector,
    static_injection_data::Vector,
    frequency_reference::Type{T},
) where {T <: Union{ConstantFrequency, ReferenceBus}}
    global_vars_dict = get_vars_ix()
    global_vars_dict[GLOBAL_VAR_SYS_FREQ_INDEX] = get_frequency_reference!(
        frequency_reference,
        wrapped_injectors,
        static_injection_data,
    )
    return global_vars_dict
end

function get_dynamic_wrapper(inputs::SimulationInputs, name::AbstractString)
    for w in get_dynamic_injectors(inputs)
        if PSY.get_name(w) == name
            return w
        end
    end
    println("DynamicWrapper with name $(name) was not found.")
    return nothing
end

function get_setpoints(inputs::SimulationInputs)
    dic = Dict{String, Dict{String, Float64}}()
    for w in get_dynamic_injectors(inputs)
        dic_w = Dict{String, Float64}()
        dic_w["P_ref"] = get_P_ref(w)
        dic_w["Q_ref"] = get_Q_ref(w)
        dic_w["ω_ref"] = get_ω_ref(w)
        dic_w["V_ref"] = get_V_ref(w)
        dic[PSY.get_name(w)] = dic_w
    end
    return dic
end


text:
# Functions used for 1d plotting, i.e. plotting data depending on only one
# parameter

function plot(gr::Graph, y::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.plot(gr.ptr, Data(y).ptr, stl, opt)
end

function plot(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.plot_xy(gr.ptr, Data(x).ptr, Data(y).ptr, stl, opt)
end

function plot(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.plot_xyz(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function radar(gr::Graph, a::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.radar(gr.ptr, Data(a).ptr, stl, opt)
end

function step(gr::Graph, y::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.step(gr.ptr, Data(y).ptr, stl, opt)
end

function step(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.step_xy(gr.ptr, Data(x).ptr, Data(y).ptr, stl, opt)
end

function step(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.step_xyz(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function tens(gr::Graph, y::Array{cmgl.Float}, c::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.tens(gr.ptr, Data(y).ptr, Data(c).ptr, stl, opt)
end

function tens(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, c::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.tens_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(c).ptr, stl, opt)
end

function tens(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}, c::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.tens_xyz(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, Data(c).ptr, stl, opt)
end

function tape(gr::Graph, y::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.tape(gr.ptr, Data(y).ptr, stl, opt)
end

function tape(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.tape_xy(gr.ptr, Data(x).ptr, Data(y).ptr, stl, opt)
end

function tape(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.tape_xyz(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function area(gr::Graph, y::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.area(gr.ptr, Data(y).ptr, stl, opt)
end

function area(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.area_xy(gr.ptr, Data(x).ptr, Data(y).ptr, stl, opt)
end

function area(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.area_xyz(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function region(gr::Graph, y1::Array{cmgl.Float}, y2::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.region(gr.ptr, Data(y1).ptr, Data(y2).ptr, stl, opt)
end

function region(gr::Graph, x::Array{cmgl.Float}, y1::Array{cmgl.Float}, y2::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.region_xy(gr.ptr, Data(x).ptr, Data(y1).ptr, Data(y2).ptr, stl, opt)
end

function region(gr::Graph, x1::Array{cmgl.Float}, x2::Array{cmgl.Float}, y1::Array{cmgl.Float}, y2::Array{cmgl.Float}, z1::Array{cmgl.Float}, z2::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.region_3d(gr.ptr, Data(x1).ptr, Data(x2).ptr, Data(y1).ptr, Data(y2).ptr, Data(z1).ptr, Data(z2).ptr, stl, opt)
end

function bars(gr::Graph, y::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.bars(gr.ptr, Data(y).ptr, stl, opt)
end

function bars(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.bars_xy(gr.ptr, Data(x).ptr, Data(y).ptr, stl, opt)
end

function bars(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.bars_xyz(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function stem(gr::Graph, y::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.stem(gr.ptr, Data(y).ptr, stl, opt)
end

function stem(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.stem_xy(gr.ptr, Data(x).ptr, Data(y).ptr, stl, opt)
end

function stem(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.stem_xyz(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function barh(gr::Graph, y::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.barh(gr.ptr, Data(y).ptr, stl, opt)
end

function barh(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.barh_xy(gr.ptr, Data(x).ptr, Data(y).ptr, stl, opt)
end

function cones(gr::Graph, y::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.cones(gr.ptr, Data(y).ptr, stl, opt)
end

function cones(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.cones_xy(gr.ptr, Data(x).ptr, Data(y).ptr, stl, opt)
end

function cones(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.cones_xyz(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, stl, opt)
end

function chart(gr::Graph, a::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.chart(gr.ptr, Data(a).ptr, stl, opt)
end

function boxplot(gr::Graph, a::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.boxplot(gr.ptr, Data(y).ptr, stl, opt)
end

function boxplot(gr::Graph, x::Array{cmgl.Float}, a::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.boxplot_xy(gr.ptr, Data(x).ptr, Data(a).ptr, stl, opt)
end

function candle(gr::Graph, v1::Array{cmgl.Float}, y1::Array{cmgl.Float}, y2::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.candle(gr.ptr, Data(v1).ptr, Data(y1).ptr, Data(y2).ptr, stl, opt)
end

function candle(gr::Graph, v1::Array{cmgl.Float}, v2::Array{cmgl.Float}, y1::Array{cmgl.Float}, y2::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.candle_yv(gr.ptr, Data(v1).ptr, Data(v2).ptr, Data(y1).ptr, Data(y2).ptr, stl, opt)
end

function candle(gr::Graph, x::Array{cmgl.Float}, v1::Array{cmgl.Float}, v2::Array{cmgl.Float}, y1::Array{cmgl.Float}, y2::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.candle_xyv(gr.ptr, Data(x).ptr, Data(v1).ptr, Data(v2).ptr, Data(y1).ptr, Data(y2).ptr, stl, opt)
end

function ohlc(gr::Graph, o::Array{cmgl.Float}, h::Array{cmgl.Float}, l::Array{cmgl.Float}, c::Array{cmgl.Float}, stl::String="", opt::String="")
    cmgl.ohcl(gr.ptr, Data(o).ptr, Data(h).ptr, Data(l).ptr, Data(c).ptr, stl, opt)
end

function ohlc(gr::Graph, x::Array{cmgl.Float}, o::Array{cmgl.Float}, h::Array{cmgl.Float}, l::Array{cmgl.Float}, c::Array{cmgl.Float}, stl::String="", opt::String="")
    cmgl.ohcl_x(gr.ptr, Data(x).ptr, Data(o).ptr, Data(h).ptr, Data(l).ptr, Data(c).ptr, stl, opt)
end

function error(gr::Graph, y::Array{cmgl.Float}, ey::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.error(gr.ptr, Data(y).ptr, Data(ey).ptr, pen, opt)
end

function error(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, ey::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.error_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(ey).ptr, pen, opt)
end

function error(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, ex::Array{cmgl.Float}, ey::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.error_exy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(ex).ptr, Data(ey).ptr, pen, opt)
end

#=function mark(gr::Graph, y::Array{cmgl.Float}; r::cmgl.Float=1, stl::String="", opt::String="")=#
    #=r_ = zeros(length(y))=#
    #=fill!(r_, r)=#
    #=cmgl.mark_y(gr.ptr, Data(y).ptr, Data(r_).ptr, pen, opt)=#
#=end=#

#=function mark(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}; r::cmgl.Float=1, stl::String="", opt::String="")=#
    #=r_ = zeros(length(y))=#
    #=fill!(r_, r)=#
    #=cmgl.mark_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(r_).ptr, pen, opt)=#
#=end=#

#=function mark(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}; r::cmgl.Float=ones(length(y)), stl::String="", opt::String="")=#
    #=r_ = zeros(length(y))=#
    #=fill!(r_, r)=#
    #=cmgl.mark_xyz(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, Data(r_).ptr, pen, opt)=#
#=end=#

function mark(gr::Graph, y::Array{cmgl.Float}, r::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.mark_y(gr.ptr, Data(y).ptr, Data(r).ptr, pen, opt)
end

function mark(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, r::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.mark_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(r).ptr, pen, opt)
end

function mark(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}, r::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.mark_xyz(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, Data(r).ptr, pen, opt)
end

#=function textmark(gr::Graph, y::Array{cmgl.Float}, txt::String; r::cmgl.Float=1, stl::String="", opt::String="")=#
    #=r_ = zeros(length(y))=#
    #=fill!(r_, r)=#
    #=cmgl.textmark_yr(gr.ptr, Data(y).ptr, Data(r_).ptr, txt, stl, opt)=#
#=end=#

function textmark(gr::Graph, y::Array{cmgl.Float}, txt::String, r::Array{cmgl.Float}=ones(length(y)); stl::String="", opt::String="")
    cmgl.textmark_yr(gr.ptr, Data(y).ptr, Data(r).ptr, txt, stl, opt)
end

function textmark(gr::Graph, y::Array{cmgl.Float}, txt::String, r::cmgl.Float; stl::String="", opt::String="")
    r_ = zeros(length(y))
    fill!(r_, r)
    cmgl.textmark_yr(gr.ptr, Data(y).ptr, Data(r_).ptr, txt, stl, opt)
end

function textmark(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, txt::String, r::Array{cmgl.Float}=ones(length(y)); stl::String="", opt::String="")
    cmgl.textmark_xyr(gr.ptr, Data(x).ptr, Data(y).ptr, Data(r).ptr, txt, stl, opt)
end

function textmark(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, txt::String, r::cmgl.Float; stl::String="", opt::String="")
    r_ = zeros(length(y))
    fill!(r_, r)
    cmgl.textmark_xyr(gr.ptr, Data(x).ptr, Data(y).ptr, Data(r_).ptr, txt, stl, opt)
end

function textmark(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}, txt::String, r::Array{cmgl.Float}=ones(length(y)); stl::String="", opt::String="")
    cmgl.textmark_xyzr(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, Data(r).ptr, txt, stl, opt)
end

function textmark(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}, txt::String, r::cmgl.Float; stl::String="", opt::String="")
    r_ = zeros(length(y))
    fill!(r_, r)
    cmgl.textmark_xyzr(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, Data(r_).ptr, txt, stl, opt)
end

function label(gr::Graph, y::Array{cmgl.Float}, txt::String; stl::String="", opt::String="")
    cmgl.label(gr.ptr, Data(y).ptr, txt, stl, opt)
end

function label(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, txt::String; stl::String="", opt::String="")
    cmgl.label_xy(gr.ptr, Data(x).ptr, Data(y).ptr, txt, stl, opt)
end

#=function label(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, txt::UTF8String; stl::String="", opt::String="")=#
    #=cmgl.labelw_xy(gr.ptr, Data(x).ptr, Data(y).ptr, txt, stl, opt)=#
#=end=#

function label(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}, txt::String; stl::String="", opt::String="")
    cmgl.label_xyz(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, txt, stl, opt)
end

#=function label(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}, txt::UTF8String; stl::String="", opt::String="")=#
    #=cmgl.labelw_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, txt, stl, opt)=#
#=end=#

function table(gr::Graph, val::Array{cmgl.Float}, txt::Array{String, 1}; stl::String="", opt::String="")
    cmgl.table(gr.ptr, 0, 0, Data(val).ptr, join(txt, '\n'), stl)
end

#=function table(gr::Graph, val::Array{cmgl.Float}, txt::Array{UTF8String, 1}; stl::String="", opt::String="")=#
    #=cmgl.tablew(gr.ptr, 0, 0, Data(val).ptr, join(txt, '\n'), stl)=#
#=end=#

function table(gr::Graph, x::cmgl.Float, y::cmgl.Float, val::Array{cmgl.Float}, txt::Array{String, 1}; stl::String="", opt::String="")
    cmgl.table(gr.ptr, x, y, Data(val).ptr, join(txt, '\n'), stl)
end

#=function table(gr::Graph, x::cmgl.Float, y::cmgl.Float, val::Array{cmgl.Float}, txt::Array{UTF8String, 1}; stl::String="", opt::String="")=#
    #=cmgl.tablew(gr.ptr, x, y, Data(val).ptr, join(txt, '\n'), stl)=#
#=end=#

function tube(gr::Graph, y::Array{cmgl.Float}, r::cmgl.Float; stl::String="", opt::String="")
    cmgl.tube(gr.ptr, Data(y).ptr, r, stl, opt)
end

function tube(gr::Graph, y::Array{cmgl.Float}, r::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.tube_r(gr.ptr, Data(y).ptr, Data(r).ptr, stl, opt)
end

function tube(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, r::cmgl.Float; stl::String="", opt::String="")
    cmgl.tube_xy(gr.ptr, Data(x).ptr, Data(y).ptr, r, stl, opt)
end

function tube(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, r::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.tube_xyr(gr.ptr, Data(x).ptr, Data(y).ptr, Data(r).ptr, stl, opt)
end

function tube(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}, r::cmgl.Float; stl::String="", opt::String="")
    cmgl.tube_xyz(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, r, stl, opt)
end

function tube(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}, r::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.tube_xyzr(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, Data(r).ptr, stl, opt)
end

function torus(gr::Graph, r::Array{cmgl.Float}, z::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.torus(gr.ptr, Data(r).ptr, Data(z).ptr, stl, opt)
end


text:
module Chemostat_EColi

    import Plots
    using ProjAssistant
    @gen_top_proj

    # ---------------------------------------------------------------------------------
    # Dat
    DATfile(src) = procdir(Chemostat_EColi, "dat", (;src), ".jls")
    load_DAT(src) = ldat(DATfile(src))

    LP_DATfile(src) = procdir(Chemostat_EColi, "lp_dat", (;src), ".jls")
    load_LP_DAT(src) = ldat(LP_DATfile(src))

    ## ------------------------------------------------------
    # ustils
    function lims(d; marginf = 0.1)
        m = minimum(d)
        M = maximum(d)
        margin = (M - m) * marginf
        (m - margin, M + margin)
    end

    function lims(d, oldl; marginf = 0.1)
        isempty(oldl) && lims(d; marginf)
        newl = lims(d; marginf)
        [min(first(newl), first(oldl)), max(last(newl), last(oldl))]
    end

end


text:
#=
fbp-par.jl
Test basic parallel-beam FBP
=#

using Sinograms: fbp
using Test: @test, @testset, @inferred
#using MIRTjim: jim # debug only

@testset "fbp-par" begin
    nr = 64
    dr = 10 / nr
    r = (-(nr-1)/2):((nr-1)/2) * dr
    na = 128
    ϕ = (0:(na-1))/na * π

    proj1 = r -> abs(r) < 1 ? sqrt(1 - r^2) : 0.
    proj2 = (r, ϕ, x, y) -> proj1(r - (x * cos(ϕ) + y * sin(ϕ)))
    sino = proj2.(r, ϕ', 3, 1)
    #jim(r, ϕ, sino)

    tmp = @inferred fbp(sino)
    @test tmp isa Matrix
end


text:
# This file is a part of DataDeps.jl. License is MIT.

"""
    progress_update_period()

Returns the period between updated being logged on the progress.
This is used by the default `fetch_method` and is generally a good idea
to use it in any custom fetch method, if possible
"""
function progress_update_period()
    envvar = get(ENV, "DATADEP_PROGRESS_UPDATE_PERIOD") do
        if haskey(ENV, "DATADEP_ALWAYS_ACCEPT")
            # Running in a script, probably want minimal updates
            "Inf"
        else # default
            "5" # seconds
        end
    end
    parse(Float32, envvar)
end

"""
    fetch_http(remotepath, localdir; update_period=5)

Pass in a HTTP[/S] URL  and a directory to save it to,
and it downloads that file, returing the local path.
This is using the HTTP protocol's method of defining filenames in headers,
if that information is present.
Returns the localpath that it was donwloaded to.

`update_period` controls how often to print the download progress to the log.
It is expressed in seconds. It is printed at `@info` level in the log.
By default it is once per second, though this depends on configuration

"""
function fetch_http(remotepath, localdir; update_period=progress_update_period())
    @assert(isdir(localdir))
    HTTP.download(remotepath, localdir; update_period=update_period)
end


text:
"""
$(TYPEDEF)

Stores  a  `NamedTuple`  of  [`Parameter`](@ref).   Encapsulates  an  autotuning
problem definition.

$(TYPEDFIELDS)
"""
struct Configuration{names, T}
    parameters::NamedTuple{names, T}
end

getindex(c::Configuration, symbol::Symbol) = getindex(c.parameters, symbol)
keys(c::Configuration) = keys(c.parameters)
values(c::Configuration) = values(c.parameters)
pairs(c::Configuration) = pairs(c.parameters)

"""
$(TYPEDEF)

A  configuration  parameter  with value  of  type  `T`  and  valid range  in  an
`Array{T,1}`. Should be created as part of a [`Configuration`](@ref).

$(TYPEDFIELDS)
"""
struct Parameter{T}
    current::T
    values::Array{T, 1}
end

values(p::Parameter) = Tuple(p.values)
getindex(p::Parameter, index::Int64) = getindex(p.values, index)

"""
$(TYPEDEF)

Primitive type, `Integer` wrapper.
"""
primitive type PowerOfTwo <: Integer 64 end

"""
$(TYPEDSIGNATURES)

Creates a  `PowerOfTwo` instance using  an `Int64`.
"""
PowerOfTwo(x::Int64) = reinterpret(PowerOfTwo, 2 ^ x)

Int64(x::PowerOfTwo) = reinterpret(Int64, x)

Base.show(io::IO, x::PowerOfTwo) = print(io, Int64(x))

"""
$(TYPEDSIGNATURES)

Receives a [`Parameter`](@ref) `p` and  returns a new [`Parameter`](@ref) with a
new `current_value`.   The new value  is chosen from a `Distribution`.

The  caller must  create the  appropriate  distribution to  control valid  value
ranges and perturbation magnitude. For `Float64` parameters, `distribution` must
be a valid continuous distribution or a `MixtureModel`.
"""
function perturb(p::Parameter{Float64}, distribution::Distribution)
    return Parameter(rand(distribution), p.values)
end

"""
$(TYPEDSIGNATURES)

Perturbing  categorical  or discrete  numerical  parameters  must be  done  with
`Distributions.Categorical` distributions.
"""
function perturb(p::Union{Parameter{PowerOfTwo}, Parameter{Int64}, Parameter{String}},
                 distribution::Distribution)
    return Parameter(p.values[rand(distribution)], p.values)
end

"""
$(TYPEDSIGNATURES)

Perturbs a [`Configuration`](@ref) by calling  [`perturb`](@ref) for each of its
[`Parameter`](@ref)  with  the  distributions   specified  in  the  `NamedTuple`
`distributions`.
"""
function perturb(c::Configuration, distribution::NamedTuple)
    Configuration(NamedTuple{keys(c)}([perturb(v, distribution[k]) for (k, v) in pairs(c)]))
end


text:
function add5(n)
    n + 5
end

function return15()
    a::Int64 = 9
    add5(a + 1)
end

function square(n)
    n * n
end

function square_string()
    square("1")
end

function add_kr_li(n, k, l)
    n + Complex(k, l)
end

println("add5, Int64:")
code_native(add5, (Int64,), syntax = :intel)

println("add5, Float64:")
code_native(add5, (Float64,), syntax = :intel)

println("return15:")
code_native(return15, (), syntax = :intel)

# println("square, Int64:")
# code_native(square, (Int64,), syntax = :intel)
#
# println("square, Float64:")
# code_native(square, (Float64,), syntax = :intel)
#
# println("square, String:")
# code_native(square, (String,), syntax = :intel)
#
# println("square_string:")
# code_native(square_string, (), syntax = :intel)

max_stars_repo_name: phrb/PDP

max_issues_repo_name: phrb/PDP

max_forks_repo_name: phrb/PDP