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

using Markdown
using InteractiveUtils

# ╔═╡ b18da173-12b7-43cf-bcaa-abac32d4390c
using Unitful, CommonMark

# ╔═╡ 8ccd67d7-193e-4a5b-9274-e1386c301c5c
begin
	using Lazy

	import MarkdownLiteral: @mdx

	macro css_str(stylesheet_url)
		the_css = open(download(stylesheet_url)) do io
			read(io, String)
		end
		return the_css
	end

end

# ╔═╡ 853f2df9-e266-4f4f-91f4-2bb9e844b74c
begin
	using LaTeXStrings,
		PlutoUI,
		HypertextLiteral,
		MarkdownLiteral,
		Hyperscript,
		Random

	macro title_str(title::String, bookID::String="")
		cm_parser = CommonMark.Parser()
		enable!(cm_parser, AttributeRule)

		bookID == "" && return title |> cm_parser

		bib = Dict(
			"ms6" 	=> 	"""
						Moran, MJ & HN Shapiro (2008)
						Fundamentals of engineering thermodynamics 6e,
						978-0471787358
						https://gitlab.com/Good-Trouble/youCAN/doThermodynamics/-/raw/main/images/ms6cover.jpg
						""",
			"ms8" 	=> 	"""
						Moran, MJ Ed (2014)
						Fundamentals of engineering thermodynamics 8e,
						978-1118412930
						https://gitlab.com/Good-Trouble/youCAN/doThermodynamics/-/raw/main/images/ms8cover.jpg
						""",
			"ms9" 	=> 	"""
						Moran, MJ, _et al_ (2018)
						Fundamentals of engineering thermodynamics 9e,
						978-1119391388
						https://gitlab.com/Good-Trouble/youCAN/doThermodynamics/-/raw/main/images/ms9cover47x60.png
						""",
			"cb7" 	=> 	"""
						Çengel, YA & MA Boles (2011)
						Thermodynamics: An engineering approach 7e,
						978-0073529325
						https://gitlab.com/Good-Trouble/youCAN/doThermodynamics/-/raw/main/images/cb7cover.jpg
						""",
			"cb8" 	=> 	"""
						Çengel, YA & MA Boles (2015)
						Thermodynamics: An engineering approach 8e,
						978-0073398174
						https://gitlab.com/Good-Trouble/youCAN/doThermodynamics/-/raw/main/images/cb8cover.jpg
						""",
			"cb9" 	=> 	"""
						Çengel, YA, Boles, MA, & M Kanoglu (2019)
						Thermodynamics: An engineering approach 9e,
						978-1259822674
						https://gitlab.com/Good-Trouble/youCAN/doThermodynamics/-/raw/main/images/cb9cover48x60.jpg
						""",
		)

		haskey(bib, bookID) || return KeyError(bookID)

		the_bibliography, the_thumb =
			@as bb bib[bookID] begin
				split(bb, '\n', keepempty=false)
				strip.(bb)
				(
					join(bb[1:3], " _", "_ "),
					m("img", class="no-dark-invert", src=bb[4], title=join(bb[2:3]," "))

				)
			end

		out = @mdx("""

\n\n
			\n\n# $title\n\n{.source}\n$the_bibliography\n

			$the_thumb

			""")

		return out

	end

	  # ═══ load css ════════════════════════════════════════════*═

	@mdx("""

Load stylesheets

	{.appendix-item}
	Load Julia \"infrastructure\" pkgs
	""")

end

# ╔═╡ 916c35d9-2007-4a4d-8e54-6c8a60583e7e
title"Problem 3.044"ms9

# ╔═╡ eab9d3c3-79ac-454b-8849-03e2f5ad1b44
# Given

begin

	T1 = -10u"°C"
	x1 =  0.8
	T2 =  40u"°C"

	V = 0.01u"m^3"

	Volt = 12u"V"
	I = 5u"A"
	t = 5u"minute"

end;

# ╔═╡ 83d5fd68-3e02-4fa0-a0a2-9a0876befa82
cm"""

## Given

As shown in **Fig. P3.44**, a closed, rigid tank fitted with a fine‐wire electric resistor is filled with Refrigerant 22, initially at $T1, a quality of $(Int(x1*100))%, and a volume of $V. A $Volt battery provides a current of $I  to the resistor for $(t)s. The final temperature of the refrigerant is $T2.

{.w-80 .center title="Textbook 'Figure P3.44' " }
$(Resource("https://imgur.com/if7XPCk.png"))

{.fig_caption }
"Figure P3.44" from the textbook.

"""

# ╔═╡ 3bd8532c-a167-4400-9e98-7e407032a834
cm"""

```math
\begin{gather}
  E_{in} &-& E_{out} &=& \Delta E_{sys} \\[6pt]

  W_{elec,in} &-& Q_{out} &=& U_2 - U_1 \\[6pt]

  (VI)t  &-& Q_{out} &=& m\left( u_2 - u_1 \right) \\[6pt] {}

\end{gather}

```

```math
\begin{equation}
   Q_{out} = (VI)t  + m\left( u_1 - u_2 \right) \\

\end{equation}

```

"""

# ╔═╡ 2a1d3c88-4381-4690-93c9-e300adaddbdb
cm"""

Take inventory. _Note that ``m`` was replaced with ``\displaystyle\frac{V}{v}`` and the properties we need are marked red._

```math
\begin{gather}
   Q_{out} = (VI)t  + \frac{V}{\color{red}v}\left(\color{red} u_1 - u_2 \right) \\[6pt]{}\\
\end{gather}

```

**State ①**. Get ``v, u`` at ``T = \pu{-10 °C}, x = 0.8``

"""

# ╔═╡ 5aa61edf-aada-4303-9598-b18b1df63521
begin

  vf_⁻10°C = 0.7606e-3u"m^3/kg"  # A-7
  vg_⁻10°C = 0.0652u"m^3/kg"     # A-7

# If you prefer, recall 1 L = 0.001 m³
  vf_⁻10°C = 0.7606u"L/kg"

  uf_⁻10°C = 33.27u"kJ/kg"       # A-7
  ug_⁻10°C = 223.02u"kJ/kg"      # A-7

end;

# ╔═╡ 30b360c2-bde1-4f1d-81bd-cf9a4abbaac1

# ╔═╡ fcc88e47-f88e-4151-8ee7-4ba328dbd5e1
v = vf_⁻10°C + x1 * (vg_⁻10°C - vf_⁻10°C)

# ╔═╡ b72069ef-76b2-44b9-9439-6251f8cb3d51

# ╔═╡ 414b5a13-6dc1-4bce-a8b4-85a48ba647cd
u1 = uf_⁻10°C + x1 * (ug_⁻10°C - uf_⁻10°C)

# ╔═╡ 1636d935-0673-4a66-84a1-651542bb420d

# ╔═╡ a20e9fd2-9540-4f15-bfa3-f32c61dd37ba
cm"""

**State ②**. Get ``u`` at ``T = \pu{-10 °C}, v =`` $(round(v |> typeof, v, digits=4) |> ustrip) ``\rm m^3/kg``.

We find that State ② is a superheated vapor, and a quick look shows we'll have to interpolate to resolve our state:

At 40 °C,

```math
\begin{gather}
   \quad \frac{v - v_{5bar}}{v_{5.5bar} - v_{5bar}} =
             \frac{u_2 - u_{5bar}}{u_{5.5bar} - u_{5bar}}  \\{}
\end{gather}
```

"""

# ╔═╡ 6f2db943-228f-4905-b3f0-97462c69e23c
# interpolate superheated vapor table

begin

  v_5bar_40°C   = 0.05636u"m^3/kg"  # A-9
  v_5ˌ5bar_40°C = 0.05086u"m^3/kg"  # A-9

  u_5bar_40°C   = 250.70u"kJ/kg"    # A-9
  u_5ˌ5bar_40°C = 250.20u"kJ/kg"    # A-9

end;

# ╔═╡ de0e0c33-4239-4050-928a-298510224c29
u2 = u_5bar_40°C +
       (v - v_5bar_40°C) / (v_5ˌ5bar_40°C - v_5bar_40°C) *
	   (u_5ˌ5bar_40°C - u_5bar_40°C)

# ╔═╡ 30c658f0-2736-4a87-a7de-bd539f919b37

# ╔═╡ a5bd2cc7-1fc0-4b14-9e86-fe2a0e97b15b

cm"""

Now we have everything we need to solve for the heat transfer, and close out the problem:

```math
\begin{equation}
   Q_{out} = (VI)t  + \frac{V}{v}\left( u_1 - u_2 \right) \\
\end{equation}

```

"""

# ╔═╡ 838a1526-67ef-4d89-84d1-a96af8321041
Q_out = Volt*I*t + V/v * (u1 - u2) |> u"kJ"

# ╔═╡ b0bd133f-7b36-4a8c-9c5c-18e1287e8e63
begin
	Q_out_rounded = round(Q_out |> typeof, Q_out, digits=4)

cm"""

## _finis_

{.message .success}
> et voilà
>
>``\qquad Q_{out}`` = $(Q_out_rounded)

"""
end

# ╔═╡ b190af88-07f0-47e6-b534-a4561ab31bb3
cm"""

## Find

Determine the heat transfer, in kJ, from the refrigerant.

{.message .found}
> found
> ``\qquad Q_{out} = `` $Q_out_rounded

"""

# ╔═╡ 11078015-d9e2-4215-a5db-f7f224a62e32

# ╔═╡ 7634fd1d-7c03-4c63-a005-9db29759154e

# ╔═╡ a4d9de7f-edac-4b4d-acda-ca656a83cee0
@mdx """

{style="border-color: var(--mcs-color-appendix); border-width: 0.5px;"}
---

### Appendix

_One of the major benefits of reactive notebooks_{.newthought} like Pluto is that allow large, unwieldy chunks of code to be placed outside the reading flow without removing them from their proper place in the computation flow. Their results can be used in calculations or spliced (interpolated) into the narrative.

For interpolation examples, see the
$(Resource("https://raw.githubusercontent.com/jupyterhub/binderhub/43ee103410b0ee9e5761171b37580281c3936eb5/binderhub/static/images/markdown-icon.svg", :height => 24))  markup throughout the notebook, perhaps even here 😉. In particular, the values in the problem statement are actually sourced in the `Julia` code under the image. Play with the values and check the results.

"""

# ╔═╡ a3a8cc93-6812-4b2e-acb6-30627a3c0c2f

# ╔═╡ 23270456-5512-4b87-9345-62f6e6ba9fb1
isochoric_process_svg = """
	<svg width="400px" viewBox="0 0 265 166" version="1.1" xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink">
	    <g id="Page-1" stroke="none" stroke-width="1" fill="none" fill-rule="evenodd">
	        <g id="Mobile" transform="translate(-16.000000, -98.000000)">
	            <g id="Group" transform="translate(16.000000, 100.000000)">
	                <line x1="131.5" y1="63.5" x2="131.5" y2="128.5" id="Line-3" stroke="#979797" stroke-width="0.5" stroke-linecap="square" stroke-dasharray="3,1"></line>
	                <path id="Line" d="M254,140 L263,144.5 L254,149 L254,145 L13,145 L13,144 L254,144 L254,140 Z" fill="#2A2D38" fill-rule="nonzero"></path>
	                <path id="Line" d="M13.5,-0.128205128 L18,8.87179487 L14,8.871 L14,145.128205 L13,145.128205 L13,8.871 L9,8.87179487 L13.5,-0.128205128 Z" fill="#2A2D38" fill-rule="nonzero"></path>
	                <path d="M29,126 C40.2466732,83.6345528 58.9681265,0 82.4098041,0 C105.851482,0 165.01311,126 263,126" id="Path" stroke="#F2590C"></path>
	                <text id="v" font-family="CMUClassicalSerif-Italic, CMU Classical Serif" font-size="12" font-style="italic" font-weight="normal" fill="#000000">
	                    <tspan x="247" y="161">v</tspan>
	                </text>
	                <text id="P" font-family="CMUClassicalSerif-Italic, CMU Classical Serif" font-size="12" font-style="italic" font-weight="normal" fill="#000000">
	                    <tspan x="0" y="17">P</tspan>
	                </text>
	                <text id="1" font-family="Helvetica-Light, Helvetica" font-size="9" font-weight="300" fill="#000000">
	                    <tspan x="123" y="74">1</tspan>
	                </text>
	                <path id="Line-2" d="M132,71.8714286 L132,108 L136,108 L131.5,117 L127,108 L131,108 L131,71.8714286 L132,71.8714286 Z" fill="#0000FF" fill-rule="nonzero"></path>
	            </g>
	        </g>
	    </g>
	</svg>
	""";

# ╔═╡ 27d9d932-e3b4-4917-b3f1-7a819af310c7
cm"""

## Get Organized

These things you should see:

1. "...closed, rigid tank..." means we have an
   _isochoric process_{.addblue}. Our system
   is also  constant mass; thus, we have
   _constant specific volume_{.killerorange}.

1. We're losing heat. Now we should have an idea of our process's flow.

   {.center .invertible .dimmer}
   $isochoric_process_svg

   {.fig_caption }
   Our process.

   _Get in the habit of sketching your processes._{.newthought} A sketch
   always tells a great story. This one tells us
   **no boundary work**{.killerorange}
   without us having to think about anything.

1. We have one work interaction, in the form of electricity.

1. We have enough information to resolve State ①.

1. Our technique for State ② is straightforward, one you'll repeat many,
   many times: resolve a state using a known adjacent state and the process
   connecting them.

   ```math
   \left.
   \begin{align}
     T_1 &= \pu{-10 °C} \\
     x_1 &= 0.8
   \end{align}
   \right\}
   ~\ce{① ->[~{\rm d}V = 0~~] ②}~
   \left\{
   \begin{aligned}
     T_2 &= \pu{40 °C} \\
     v_2 &= v_1
   \end{aligned}
   \right.
   ```

"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CommonMark = "a80b9123-70ca-4bc0-993e-6e3bcb318db6"
Hyperscript = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
HypertextLiteral = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
Lazy = "50d2b5c4-7a5e-59d5-8109-a42b560f39c0"
MarkdownLiteral = "736d6165-7244-6769-4267-6b50796e6954"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"

[compat]
CommonMark = "~0.8.3"
Hyperscript = "~0.0.4"
HypertextLiteral = "~0.9.3"
LaTeXStrings = "~1.3.0"
Lazy = "~0.15.1"
MarkdownLiteral = "~0.1.1"
PlutoUI = "~0.7.16"
Unitful = "~1.9.0"
"""

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

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

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

[[CommonMark]]
deps = ["Crayons", "JSON", "URIs"]
git-tree-sha1 = "393ac9df4eb085c2ab12005fc496dae2e1da344e"
uuid = "a80b9123-70ca-4bc0-993e-6e3bcb318db6"
version = "0.8.3"

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

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

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

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

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

[[HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

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

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

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

[[Lazy]]
deps = ["MacroTools"]
git-tree-sha1 = "1370f8202dac30758f3c345f9909b97f53d87d3f"
uuid = "50d2b5c4-7a5e-59d5-8109-a42b560f39c0"
version = "0.15.1"

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

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

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

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

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

[[MarkdownLiteral]]
deps = ["CommonMark", "HypertextLiteral"]
git-tree-sha1 = "0d3fa2dd374934b62ee16a4721fe68c418b92899"
uuid = "736d6165-7244-6769-4267-6b50796e6954"
version = "0.1.1"

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

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

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

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

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

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

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

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

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

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

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

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

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

[[Unitful]]
deps = ["ConstructionBase", "Dates", "LinearAlgebra", "Random"]
git-tree-sha1 = "a981a8ef8714cba2fd9780b22fd7a469e7aaf56d"
uuid = "1986cc42-f94f-5a68-af5c-568840ba703d"
version = "1.9.0"

[[libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
"""

# ╔═╡ Cell order:
# ╟─916c35d9-2007-4a4d-8e54-6c8a60583e7e
# ╟─83d5fd68-3e02-4fa0-a0a2-9a0876befa82
# ╠═b18da173-12b7-43cf-bcaa-abac32d4390c
# ╠═eab9d3c3-79ac-454b-8849-03e2f5ad1b44
# ╟─b190af88-07f0-47e6-b534-a4561ab31bb3
# ╟─27d9d932-e3b4-4917-b3f1-7a819af310c7
# ╟─3bd8532c-a167-4400-9e98-7e407032a834
# ╟─2a1d3c88-4381-4690-93c9-e300adaddbdb
# ╠═5aa61edf-aada-4303-9598-b18b1df63521
# ╟─30b360c2-bde1-4f1d-81bd-cf9a4abbaac1
# ╠═fcc88e47-f88e-4151-8ee7-4ba328dbd5e1
# ╟─b72069ef-76b2-44b9-9439-6251f8cb3d51
# ╠═414b5a13-6dc1-4bce-a8b4-85a48ba647cd
# ╟─1636d935-0673-4a66-84a1-651542bb420d
# ╟─a20e9fd2-9540-4f15-bfa3-f32c61dd37ba
# ╠═6f2db943-228f-4905-b3f0-97462c69e23c
# ╠═de0e0c33-4239-4050-928a-298510224c29
# ╟─30c658f0-2736-4a87-a7de-bd539f919b37
# ╟─a5bd2cc7-1fc0-4b14-9e86-fe2a0e97b15b
# ╠═838a1526-67ef-4d89-84d1-a96af8321041
# ╟─b0bd133f-7b36-4a8c-9c5c-18e1287e8e63
# ╟─11078015-d9e2-4215-a5db-f7f224a62e32
# ╟─7634fd1d-7c03-4c63-a005-9db29759154e
# ╟─a4d9de7f-edac-4b4d-acda-ca656a83cee0
# ╟─a3a8cc93-6812-4b2e-acb6-30627a3c0c2f
# ╟─23270456-5512-4b87-9345-62f6e6ba9fb1
# ╟─8ccd67d7-193e-4a5b-9274-e1386c301c5c
# ╟─853f2df9-e266-4f4f-91f4-2bb9e844b74c
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: mdAshford/youCANdoThermodynamics

max_issues_repo_name: mdAshford/youCANdoThermodynamics

max_forks_repo_name: mdAshford/youCANdoThermodynamics


text:
# This file was generated, do not modify it. # hide
confidence_intervals(e_tuned_iterated_pipe)

ext: jl

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using Downloads: download
using Flux
using BSON: load
network = joinpath(@__DIR__, "modelnnpp-0.904-0.0005.bson")
download("https://gist.github.com/JuliaOd/97728c2c15d6a7255ced6e46e3a605b6/raw/modelnnpp-0.904-0.0005.bson", network)
model2d = load(network, @__MODULE__)[:model2d]

using Random: seed!
seed!(0)

using OrdinaryDiffEq
using Trixi

# This elixir was one of the setups used in the following master thesis:
# - Julia Odenthal (2021)
#   Shock capturing with artifical neural networks
#   University of Cologne, advisors: Gregor Gassner, Michael Schlottke-Lakemper
# This motivates the particular choice of fluxes, mesh resolution etc.

###############################################################################
# semidiscretization of the compressible Euler equations
gamma = 1.4
equations = CompressibleEulerEquations2D(gamma)

"""
    initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D)

A version of the classical Kelvin-Helmholtz instability based on
https://rsaa.anu.edu.au/research/established-projects/fyris/2-d-kelvin-helmholtz-test.
"""
function initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D)
  # change discontinuity to tanh
  # typical resolution 128^2, 256^2
  # domain size is [-0.5,0.5]^2
  dens0 = 1.0 # outside density
  dens1 = 2.0 # inside density
  velx0 = -0.5 # outside velocity
  velx1 = 0.5 # inside velocity
  slope = 50 # used for tanh instead of discontinuous initial condition
  # pressure equilibrium
  p     = 2.5
  #  y velocity v2 is only white noise
  v2  = 0.01*(rand(Float64,1)[1]-0.5)
  # density
  rho = dens0 + (dens1-dens0) * 0.5*(1+(tanh(slope*(x[2]+0.25)) - (tanh(slope*(x[2]-0.25)) + 1)))
  #  x velocity is also augmented with noise
  v1 = velx0 + (velx1-velx0) * 0.5*(1+(tanh(slope*(x[2]+0.25)) - (tanh(slope*(x[2]-0.25)) + 1)))+0.01*(rand(Float64,1)[1]-0.5)
  return prim2cons(SVector(rho, v1, v2, p), equations)
end
initial_condition = initial_condition_kelvin_helmholtz_instability

surface_flux = flux_lax_friedrichs
volume_flux  = flux_chandrashekar
polydeg = 3
basis = LobattoLegendreBasis(polydeg)
indicator_sc = IndicatorNeuralNetwork(equations, basis,
                                      indicator_type=NeuralNetworkPerssonPeraire(),
                                      alpha_max=0.002,
                                      alpha_min=0.0001,
                                      alpha_smooth=true,
                                      alpha_continuous=true,
                                      alpha_amr=false,
                                      variable=density_pressure,
                                      network=model2d)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                 volume_flux_dg=volume_flux,
                                                 volume_flux_fv=surface_flux)
solver = DGSEM(basis, surface_flux, volume_integral)

coordinates_min = (-0.5, -0.5)
coordinates_max = ( 0.5,  0.5)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=5,
                n_cells_max=100_000)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 5.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

amr_indicator = IndicatorNeuralNetwork(semi,
                                       indicator_type=NeuralNetworkPerssonPeraire(),
                                       alpha_max=1.0,
                                       alpha_min=0.0001,
                                       alpha_smooth=false,
                                       alpha_continuous=true,
                                       alpha_amr=true,
                                       variable=density_pressure,
                                       network=model2d)
amr_controller = ControllerThreeLevel(semi, amr_indicator,
                                      base_level=4,
                                      med_level=6, med_threshold=0.3, # med_level = current level
                                      max_level=7, max_threshold=0.5)
amr_callback = AMRCallback(semi, amr_controller,
                           interval=1,
                           adapt_initial_condition=true,
                           adapt_initial_condition_only_refine=true)

stepsize_callback = StepsizeCallback(cfl=1.3)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        amr_callback, stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
module Metadata

export outdir, projectdir

outdir = "D:/MyTemp/tmp_exp_out/"
println("metadata.jl > outdir = ", outdir)

projectdir = "D:/MyTemp/gitlocal/bayesian-pccf-code/"
println("metadata.jl > projectdir = ", projectdir)

end

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
proto_default(::Type) = missing
proto_default(::Type{T}) where {T<:Number} = zero(T)
proto_default(::Type{T}) where {T<:AbstractString} = ""
proto_default(::Type{Bool}) = false
proto_default(::Type{Vector{T}}) where T = T[]

# this method constructs a UIntXX from a string of 0's and 1's and
# spaces (which are skipped)
function frombits_str(s, ::Type{T}) where T <: Unsigned
	i = zero(T)
	for c in s
		if c == ' '
			continue
		end
		i = i << 1
		if c == '1'
			i = i | one(T)
		end
	end
	i
end

# some string macros for diffent sizes of UInts

macro frombits8_str(s)
    frombits_str(s, UInt8)
end

macro frombits64_str(s)
    frombits_str(s, UInt64)
end


text:
### extirpolation.jl
#
# Copyright (C) 2016 Astropy Developers
# Copyright (C) 2016, 2017 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: periodogram, lomb scargle, extirpolation
#
# This file is a part of LombScargle.jl.
#
# License is BSD 3-clause "New" or "Revised".
#
### Commentary:
#
# Utility functions used for implementing the fast algorithm proposed here:
# * Press, W. H., Rybicki, G. B. 1989, ApJ, 338, 277 (URL:
#   http://dx.doi.org/10.1086/167197,
#   Bibcode: http://adsabs.harvard.edu/abs/1989ApJ...338..277P)
# The actual implementation is adapted from Astropy project.  See
# stats/lombscargle/implementations/utils.py in Astropy source code.
#
### Code:

function add_at!(arr::AbstractVector{T1}, ind::AbstractVector{<:Integer},
                 vals::AbstractVector{T2}) where {T1,T2}
    @inbounds for i in eachindex(ind)
        arr[ind[i]] += vals[i]
    end
end

function extirpolate!(result, x::AbstractVector{<:Real}, y::AbstractVector{NU},
                      N::Integer, M::Integer=4) where {NU<:Number}
    # Now use legendre polynomial weights to populate the results array; This is
    # an efficient recursive implementation (See Press et al. 1989)
    fill!(result, zero(NU))
    # first take care of the easy cases where x is an integer
    integers = findall(isinteger, x)
    add_at!(result, mod.(trunc.(Int, x[integers]), N) .+ 1, y[integers])
    deleteat!(x, integers)
    deleteat!(y, integers)
    # For each remaining x, find the index describing the extirpolation range.
    # i.e. ilo[i] < x[i] < ilo[i] + M with x[i] in the center, adjusted so that
    # the limits are within the range 0...N
    ilo = clamp.(trunc.(Int, x .- div(M, 2)), 0, N - M)
    v = collect(0:(M - 1))
    numerator = [y[j] * prod(x[j] - ilo[j] .- v) for j in eachindex(x)]
    denominator = gamma(M)
    ilo .+= M .+ 1
    @inbounds for j in v
        if j > 0
            denominator *= j/(j - M)
        end
        ilo .-= 1
        add_at!(result, ilo, numerator ./ (denominator .* (x .- ilo .+ 1)))
    end
    return result
end

function trig_sum!(grid, fftgrid, bfft_vec, bfft_plan, t::AbstractVector{<:Real},
                   h::AbstractVector{<:Real}, df::Real, N::Integer,
                   Nfft::Integer, t0::Real, f0::Real=0.0,
                   freq_factor::Integer=1, Mfft::Integer=4)
    df *= freq_factor
    f0 *= freq_factor
    if f0 > 0
        H = h .* cis.(f0 .* (t .- t0) .* 2 .* pi)
    else
        H = complex(h)
    end
    tnorm = mod.(((t .- t0) .* Nfft .* df), Nfft)
    extirpolate!(grid, tnorm, H, Nfft, Mfft)
    mul!(bfft_vec, bfft_plan, grid)
    fftgrid .= @view(bfft_vec[1:N])
    if t0 != 0
        fftgrid .*= cis.(t0 .* (f0 .+ df .* (0:(N - 1))) .* 2 .* pi)
    end
    return real(fftgrid), imag(fftgrid)
end


text:
# routines related to the PETSc backend

export initialize_backend, petsc_data

"""
    Stores key PETSc information with the `grid` structure
"""
struct petsc_data
    petsclib
    da
end
petsc_data() = petsc_data(nothing,nothing)

"""
    check_backend(backend{BackendPETSc}; Scalar=Float64, dim=1)

checks PETSc backend for the `Scalar`, which loads PETSc (and, if requested, MPI)
"""
function check_backend(b::Backend{BackendPETSc}; dim=1, Scalar=Float64)

    if !isdefined(Main, :PETSc)
        error("PETSc is not loaded; ensure it is loaded first with: using PETSc")
    end
    if (b.mpi==true && !isdefined(Main, :MPI))
        error("MPI is not loaded; ensure it is loaded first with: using MPI")
    end
    @eval using PETSc

    # get the PETSc lib with our chosen `PetscScalar` type
    petsclib = PETSc.getlib(; PetscScalar = Scalar)

    # Initialize PETSc
    PETSc.initialize(petsclib)

    return petsclib
end

# Internal function that translates the CompGrids topology to PETSc nomenclature
function bcs_translate(bound)
    bound_petsc = PETSc.LibPETSc.DM_BOUNDARY_NONE
    if bound == Ghosted
        bound_petsc = PETSc.LibPETSc.DM_BOUNDARY_GHOSTED
    elseif bound == Periodic
        bound_petsc = PETSc.DM_BOUNDARY_PERIODIC
    end
    return bound_petsc
end

function globalfromlocalsize(size, localsize, opts, stencilwidth, topology, b::Backend{BackendPETSc, FT}) where FT

    if !b.mpi & !isempty(localsize)
        size = localsize   # not using MPI
    elseif b.mpi & !isempty(localsize)
        error("The PETSc backend currently doesn't allow specifying the local grid size")
    else
        size=size
    end

    return size
end


text:
using FactCheck
f(a) = error("Doge")
facts() do
  @fact f(1) --> 1 "TestName"
end


text:
using DataStructures

part1, part2 = 0, 0
for line in readlines("input.txt")
    range, char, password = split(line)

    l, u = parse.(Int, split(range, "-"))
    c = first(char)

    cntr = counter(collect(password))

    part1 += l ≤ cntr[c] ≤ u
    part2 += (password[l] == c) ⊻ (password[u] == c)
end

@info "Solutions:" part1 part2

max_stars_repo_name: dglmoore/advent-of-code

max_issues_repo_name: dglmoore/advent-of-code

max_forks_repo_name: dglmoore/advent-of-code


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

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using QuadGK, Plots
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.3"

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
begin
	sup = (-1,1)
	f1(x) = 3/4*(1-x^2)
	f2(x) = x < 0 ? x+1 : 1-x

	expect(f,support) = quadgk((x) -> x*f(x),support...)[1]
end

# ╔═╡ 3a837612-1f94-11eb-315f-0d9296c7ba90
begin
	x = -1:0.01:1
	plot(x, [f1.(x), f2.(x)], lab=["f1" "f2"])
end

# ╔═╡ aa0bdedc-1f93-11eb-12a7-f3ad5370c280
Text("Mean 1: $(expect(f1,sup)), Mean 2: $(expect(f2,sup))")

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.3"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╠═3a837612-1f94-11eb-315f-0d9296c7ba90
# ╠═aa0bdedc-1f93-11eb-12a7-f3ad5370c280
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
#Permutations

export
    levicivita,
    multiset_permutations,
    nthperm!,
    nthperm,
    parity,
    permutations

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

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

Base.length(p::Permutations) = (0 <= p.t <= length(p.a)) ? factorial(length(p.a), length(p.a)-p.t) : 0

"""
    permutations(a)

Generate all permutations of an indexable object `a` in lexicographic order. Because the number of permutations
can be very large, this function returns an iterator object.
Use `collect(permutations(a))` to get an array of all permutations.
"""
permutations(a) = Permutations(a, length(a))

"""
    permutations(a, t)

Generate all size `t` permutations of an indexable object `a`.
"""
function permutations(a, t::Integer)
    if t < 0
        t = length(a) + 1
    end
    Permutations(a, t)
end

function Base.iterate(p::Permutations, s = collect(1:length(p.a)))
    (!isempty(s) && max(s[1], p.t) > length(p.a) || (isempty(s) && p.t > 0)) && return
    nextpermutation(p.a, p.t ,s)
end

function nextpermutation(m, t, state)
    perm = [m[state[i]] for i in 1:t]
    n = length(state)
    if t <= 0
        return(perm, [n+1])
    end
    s = copy(state)
    if t < n
        j = t + 1
        while j <= n &&  s[t] >= s[j]; j+=1; end
    end
    if t < n && j <= n
        s[t], s[j] = s[j], s[t]
    else
        if t < n
            reverse!(s, t+1)
        end
        i = t - 1
        while i>=1 && s[i] >= s[i+1]; i -= 1; end
        if i > 0
            j = n
            while j>i && s[i] >= s[j]; j -= 1; end
            s[i], s[j] = s[j], s[i]
            reverse!(s, i+1)
        else
            s[1] = n+1
        end
    end
    return (perm, s)
end

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

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

function Base.length(c::MultiSetPermutations)
    t = c.t
    if t > length(c.ref)
        return 0
    end
    if t > 20
        g = [factorial(big(i)) for i in 0:t]
    else
        g = [factorial(i) for i in 0:t]
    end
    p = [g[t+1]; zeros(Float64,t)]
    for i in 1:length(c.f)
        f = c.f[i]
        if i == 1
            for j in 1:min(f, t)
                p[j+1] = g[t+1]/g[j+1]
            end
        else
            for j in t:-1:1
                q = 0
                for k in (j+1):-1:max(1,j+1-f)
                    q += p[k]/g[j+2-k]
                end
                p[j+1] = q
            end
        end
    end
    return round(Int, p[t+1])
end

"""
    multiset_permutations(m, f, t)

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

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

function Base.iterate(p::MultiSetPermutations, s = p.ref)
    (!isempty(s) && max(s[1], p.t) > length(p.ref) || (isempty(s) && p.t > 0)) && return
    nextpermutation(p.m, p.t, s)
end

"""
    nthperm!(a, k)

In-place version of [`nthperm`](@ref); the array `a` is overwritten.
"""
function nthperm!(a::AbstractVector, k::Integer)
    k -= 1 # make k 1-indexed
    k < 0 && throw(ArgumentError("permutation k must be ≥ 0, got $k"))
    n = length(a)
    n == 0 && return a
    f = factorial(oftype(k, n-1))
    for i=1:n-1
        j = div(k, f) + 1
        k = k % f
        f = div(f, n-i)

        j = j+i-1
        elt = a[j]
        for d = j:-1:i+1
            a[d] = a[d-1]
        end
        a[i] = elt
    end
    a
end

"""
    nthperm(a, k)

Compute the `k`th lexicographic permutation of the vector `a`.
"""
nthperm(a::AbstractVector, k::Integer) = nthperm!(collect(a), k)

"""
    nthperm(p)

Return the integer `k` that generated permutation `p`. Note that
`nthperm(nthperm([1:n], k)) == k` for `1 <= k <= factorial(n)`.
"""
function nthperm(p::AbstractVector{<:Integer})
    isperm(p) || throw(ArgumentError("argument is not a permutation"))
    k, n = 1, length(p)
    for i = 1:n-1
        f = factorial(n-i)
        for j = i+1:n
            k += ifelse(p[j] < p[i], f, 0)
        end
    end
    return k
end

# Parity of permutations

const levicivita_lut = cat([0 0  0;  0 0 1; 0 -1 0],
                           [0 0 -1;  0 0 0; 1  0 0],
                           [0 1  0; -1 0 0; 0  0 0]; dims=3)

"""
    levicivita(p)

Compute the Levi-Civita symbol of a permutation `p`. Returns 1 if the permutation
is even, -1 if it is odd, and 0 otherwise.

The parity is computed by using the fact that a permutation is odd if and
only if the number of even-length cycles is odd.
"""
function levicivita(p::AbstractVector{<:Integer})
    n = length(p)

    if n == 3
        @inbounds valid = (0 < p[1] <= 3) * (0 < p[2] <= 3) * (0 < p[3] <= 3)
        return valid ? levicivita_lut[p[1], p[2], p[3]] : 0
    end

    todo = trues(n)
    first = 1
    cycles = flips = 0

    while cycles + flips < n
        first = coalesce(findnext(todo, first), 0)
        (todo[first] = !todo[first]) && return 0
        j = p[first]
        (0 < j <= n) || return 0
        cycles += 1
        while j ≠ first
            (todo[j] = !todo[j]) && return 0
            j = p[j]
            (0 < j <= n) || return 0
            flips += 1
        end
    end

    return iseven(flips) ? 1 : -1
end

"""
    parity(p)

Compute the parity of a permutation `p` using the [`levicivita`](@ref) function,
permitting calls such as `iseven(parity(p))`. If `p` is not a permutation then an
error is thrown.
"""
function parity(p::AbstractVector{<:Integer})
    epsilon = levicivita(p)
    epsilon == 0 && throw(ArgumentError("Not a permutation"))
    epsilon == 1 ? 0 : 1
end


text:
using Random, JuMP, JuDGE, Test

include("solvers/setup_gurobi.jl")
#include("solvers/setup_cplex.jl")
#include("solvers/setup_coin.jl")
#include("solvers/setup_glpk.jl")

function cutting_stock(;seed::Int=200,L=10,sizes=[2,3,5])
	mytree = narytree(0,0)

	Random.seed!(seed)

	n=length(sizes)
	demand=rand(10:200,n)
	num_sizes=length(sizes)
	max_repeat=L÷minimum(sizes)
	max_patterns=2

	mytree.ext[:sum_max]=sum(ceil(demand[i]/(L ÷ sizes[i])) for i in 1:n)#max_patterns
	mytree.ext[:sum_min]=1
	mytree.ext[:col_max]=maximum(ceil(demand[i]/(L ÷ sizes[i])) for i in 1:n)

	function initialise_columns(judge::JuDGEModel)
		patterns=[]
		for i in 1:n
			z=zeros(n)
			z[i]=L ÷ sizes[i]
			push!(patterns,z)
		end

		for p in patterns
			sp=sub_problems(judge.tree,pattern=p,mipgap=0.0)
			sp.ext[:form]=:binary
			sp.ext[:objective]=@variable(sp, obj)
			sp.ext[:objective_con]=@constraint(sp, sp.ext[:objective]-objective_function(sp) == 0)
			sp.ext[:discrete]=VariableRef[]
			sp.ext[:discrete_branch]=VariableRef[]
			optimize!(sp)
			col=JuDGE.add_column(judge.master_problem, sp, judge.tree)
		end
	end

	function sub_problems(node;pattern=nothing,mipgap=0.00)
		sp=JuMP.Model(JuDGE_SP_Solver)
		max_branch=L*4
		@shutdown(sp, 0<=made[1:num_sizes]<=L÷minimum(sizes), Int,lb=0.0,ub=1000)

		@variable(sp, dif[1:num_sizes,1:max_branch]>=0)
		@variable(sp, val[1:num_sizes,1:max_branch])
		@variable(sp, switch[1:num_sizes,1:max_branch], Bin)

		@constraint(sp,dif1[i in 1:num_sizes, j in 1:max_branch], dif[i,j]<=switch[i,j])
		@constraint(sp,dif2[i in 1:num_sizes, j in 1:max_branch], dif[i,j]<=made[i]-val[i,j]+(L÷minimum(sizes))*(1-switch[i,j]))

		@constraint(sp,block[j in 1:max_branch],sum(dif[i,j] for i in 1:num_sizes)>=1)

		@constraint(sp,sum(sizes[i]*made[i] for i in 1:num_sizes)<=L)

		if typeof(pattern)!=Nothing
			@constraint(sp, set[i=1:num_sizes],made[i]==pattern[i])
		end

		@objective(sp,Min,1.0)
	    sp
	end

	function meet_demand(model,tree)
		@constraint(model,meet_demand[i in 1:num_sizes],made[tree][i]>=demand[i])
	end

	function pattern_branch(jmodel::JuDGEModel, inttol::Float64)
		master=jmodel.master_problem
		subproblems=jmodel.sub_problems
		tree=jmodel.tree

		index=nothing
		maxvalue=inttol
		for col in master.ext[:columns][tree]
			val=value(col.var)+0.5
			v=val > ceil(val)-inttol ? 0 : val-floor(val)
			if v>maxvalue
				maxvalue=v
				index=col
			end
		end

		if index!=nothing
			temp=index.coeffs[:made]

			branch1=JuDGE.BranchConstraint(index.var,:ge,ceil(value(index.var)),master)

			branch2 = JuDGE.BranchConstraint[]
			push!(branch2,JuDGE.BranchConstraint(index.var,:le,floor(value(index.var)),master))
			j=length(jmodel.ext[:branches])+1
			for i in 1:num_sizes
				if haskey(temp,i)
					push!(branch2,JuDGE.BranchConstraint(subproblems[tree][:val][i,j],:eq,temp[i],tree))
				else
					push!(branch2,JuDGE.BranchConstraint(subproblems[tree][:val][i,j],:eq,0.0,tree))
				end
			end
			f=function f(col::JuDGE.Column)
				at_least_one_less=false

				for i in 1:num_sizes
					if !haskey(col.coeffs[:made],i) && !haskey(temp,i)
						continue
					end
					if !haskey(temp,i) || (haskey(col.coeffs[:made],i) && col.coeffs[:made][i]>temp[i])
						return
					elseif !haskey(col.coeffs,i) || (haskey(temp,i) && col.coeffs[:made][i]<temp[i])
						at_least_one_less=true
					end
				end
				if at_least_one_less
					return :ban
				else
					return
				end
			end
			println(f(index))
			return [JuDGE.Branch(branch1),JuDGE.Branch(branch2,f)]
		end
	end

	function print_patterns(jmodel::JuDGEModel)
		active = JuDGE.get_active_columns(judy,inttol=10^-6)
		println("")
		print("Demands ")
		for i in 1:n
			print(" ")
			print(sizes[i])
			print(": ")
			print(demand[i])
			if i!=n
				print(", ")
			end
		end
		println("\n")
		for (col,count) in active[mytree]
			print("[")
			for (index,number) in col.coeffs[:made]
				for n in 1:number
					print(" ")
					print(sizes[index])
					print(" ")
				end
			end
			print("]: ")
			println(count)
		end
	end

	judy = JuDGEModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_MP_Solver,sideconstraints=meet_demand,check=false)
	initialise_columns(judy)

	judy=JuDGE.branch_and_price(judy,termination=Termination(abstol=0.99),branch_method=pattern_branch,verbose=1)
	print_patterns(judy)

	JuDGE.get_objval(judy)
end

@test cutting_stock() == 70.0
@test cutting_stock(seed=12,L=100,sizes=[2,3,5,7,11,13,17,19,24,37,47,62]) == 238.0


text:
## User should just use TransferFunction
immutable SisoRational <: SisoTf
    num::Poly{Float64}
    den::Poly{Float64}
    function SisoRational(num::Poly{Float64}, den::Poly{Float64})
        if all(den == zero(den))
            error("Zero denominator")
        elseif all(num == zero(num))
            # The numerator is zero, make the denominator 1
            den = one(den)
        end
        new(num, den)
    end
end
SisoRational(num::Vector, den::Vector) = SisoRational(Poly(map(Float64,num)), Poly(map(Float64,den)))

function minreal(sys::SisoRational, eps::Real=sqrt(eps()))
    return SisoRational(minreal(SisoZpk(sys), eps))
end

function print_siso(io::IO, t::SisoRational, var=:s)
    # Convert the numerator and denominator to strings
    numstr = sprint(print_poly, t.num, var)
    denstr = sprint(print_poly, t.den, var)

    # Figure out the length of the separating line
    len_num = length(numstr)
    len_den = length(denstr)
    dashcount = max(len_num, len_den)

    # Center the numerator or denominator
    if len_num < dashcount
        numstr = "$(repeat(" ", div(dashcount - len_num, 2)))$numstr"
    else
        denstr = "$(repeat(" ", div(dashcount - len_den, 2)))$denstr"
    end
    println(io, numstr)
    println(io, repeat("-", dashcount))
    println(io, denstr)
end

function print_compact(io::Base.IO, t::SisoRational, var=:s)
    numstr = sprint(print_poly, t.num, var)
    denstr = sprint(print_poly, t.den, var)
    println(io, "($numstr)/($denstr)")
end

Base.promote_rule{T<:Real}(::Type{SisoRational}, ::Type{T}) = SisoRational
Base.convert(::Type{SisoRational}, b::Real) = SisoRational([b], [1])

Base.zero(::Type{SisoRational}) = SisoRational(zero(Poly{Float64}), one(Poly{Float64}))
Base.zero(::SisoRational) = Base.zero(SisoRational)

Base.length(t::SisoRational) = max(length(t.num), length(t.den))

function Base.num(t::SisoRational)
    lt = length(t)
    n = zeros(lt)
    n[(lt - length(t.num) + 1):end] = t.num[:]
    return n
end

function Base.den(t::SisoRational)
    lt = length(t)
    d = zeros(lt)
    d[(lt - length(t.den) + 1):end] = t.den[:]
    return d
end

denpoly(G::SisoRational) = Poly(den(G))
numpoly(G::SisoRational) = Poly(num(G))

function evalfr(sys::SisoRational, s::Number)
    S = promote_type(typeof(s), Float64)
    den = polyval(sys.den, s)
    if den == zero(S)
        convert(S, Inf)
    else
        polyval(sys.num, s)/den
    end
end

tzero(t::SisoRational) = roots(t.num)
pole(t::SisoRational) = roots(t.den)

==(t1::SisoRational, t2::SisoRational) = (t1.num == t2.num && t1.den == t2.den)
# We might want to consider alowing scaled num and den as equal
function isapprox(t1::SisoRational, t2::SisoRational; rtol::Real=sqrt(eps()), atol::Real=0)
    isapprox(t1.num,t2.num, rtol=rtol, atol=atol) && isapprox(t1.den, t2.den, rtol=rtol, atol=atol)
end

+(t1::SisoRational, t2::SisoRational) = SisoRational(t1.num*t2.den + t2.num*t1.den, t1.den*t2.den)
+(t::SisoRational, n::Real) = SisoRational(t.num + n*t.den, t.den)
+(n::Real, t::SisoRational) = t + n
.+(t::SisoRational, n::Real) = t + n
.+(n::Real, t::SisoRational) = t + n

-(t1::SisoRational, t2::SisoRational) = SisoRational(t1.num*t2.den - t2.num*t1.den, t1.den*t2.den)
-(n::Real, t::SisoRational) = SisoRational(n*t.den - t.num, t.den)
-(t::SisoRational, n::Real) = +(t, -n)
.-(t::SisoRational, n::Real) = -(t, n)
.-(n::Real, t::SisoRational) = -(n, t)

-(t::SisoRational) = SisoRational(-t.num, t.den)

*(t1::SisoRational, t2::SisoRational) = SisoRational(t1.num*t2.num, t1.den*t2.den)
*(t::SisoRational, n::Real) = SisoRational(t.num*n, t.den)
*(n::Real, t::SisoRational) = *(t, n)
.*(t1::SisoRational, t2::SisoRational) = *(t1, t2)
.*(t::SisoRational, n::Real) = *(t, n)
.*(n::Real, t::SisoRational) = *(t, n)

/(n::Real, t::SisoRational) = SisoRational(n*t.den, t.num)
/(t::SisoRational, n::Real) = t*(1/n)
/(t1::SisoRational, t2::SisoRational) = t1*(1/t2)


text:
println("Imports")
using AdventOfCode2018
using AdventOfCode2018.Day3

println("getting input")
input = get_aoc_input(3)

println("start stuff")
rectangles = parse_input(input)
used_patches = build_used_patches(rectangles)

println("Solution 1: ", part_1(used_patches))
println("Solution 2: ", part_2(rectangles, used_patches))


text:
function pairs(from_memory::Bool)
    from_memory ? String["EUR_TRY","USD_JPY","CHF_ZAR","NZD_USD","USD_THB","GBP_CHF","EUR_SEK","USD_SGD","EUR_SGD","AUD_CHF","CHF_HKD","SGD_HKD","CHF_JPY","HKD_JPY","EUR_HUF","USD_CAD","USD_CHF","CAD_SGD","AUD_CAD","USD_SEK","GBP_SGD","GBP_HKD","NZD_CAD","AUD_USD","EUR_CAD","NZD_JPY","EUR_NZD","GBP_PLN","CAD_JPY","NZD_SGD","CAD_HKD","EUR_NOK","USD_SAR","GBP_CAD","USD_HUF","GBP_AUD","USD_PLN","GBP_USD","USD_MXN","AUD_JPY","EUR_CHF","AUD_SGD","EUR_CZK","GBP_ZAR","EUR_USD","EUR_ZAR","USD_ZAR","SGD_CHF","USD_CNH","GBP_JPY","USD_TRY","ZAR_JPY","USD_DKK","USD_HKD","EUR_HKD","AUD_NZD","AUD_HKD","TRY_JPY","EUR_AUD","EUR_PLN","EUR_DKK","USD_NOK","SGD_JPY","NZD_HKD","GBP_NZD","USD_CZK","EUR_JPY","EUR_GBP","NZD_CHF","CAD_CHF"] : pairs()
end

function pairs()
    data = Requests.get("https://api-fxpractice.oanda.com/v3/accounts/" * ACCOUNT * "/instruments"; headers = HEADERS)
    [i["name"] for i in Requests.json(data)["instruments"]]
end

max_stars_repo_name: cadwallader/Oanda

max_issues_repo_name: cadwallader/Oanda

max_forks_repo_name: cadwallader/Oanda


text:
SCORE_PARAMS = """  normalize::Function
A function to reduce a matrix, where each row represents an instance and each column a score of specific detector, to a
vector of scores for each instance. See [`scale_minmax`](@ref) for a specific implementation.

    combine::Function
A function to reduce a matrix, where each row represents an instance and each column represents the score of specific
detector, to a vector of scores for each instance. See [`combine_mean`](@ref) for a specific implementation."""

"""
    ScoreTransformer(combine = combine,
                     normalize = normalize)

Transform the results of a single or multiple outlier detection models to combined and normalized scores.

Parameters
----------
$SCORE_PARAMS
"""
@detector mutable struct ScoreTransformer <: MLJ.Static
    normalize::Function = scale_minmax
    combine::Function = combine_mean
end

"""
    ProbabilisticTransformer(combine = combine,
                             normalize = normalize)

Transform the results of a single or multiple outlier detection models to combined univariate finite distributions.

Parameters
----------
$SCORE_PARAMS
"""
@detector mutable struct ProbabilisticTransformer <: MLJ.Static
    normalize::Function = scale_minmax
    combine::Function = combine_mean
end

default_percentile_threshold = classify_quantile(DEFAULT_THRESHOLD)
"""
    DeterministicTransformer(combine = combine,
                             normalize = normalize,
                             classify = classify_quantile(DEFAULT_THRESHOLD))

Transform the results of a single or multiple outlier detection models to combined categorical values.

Parameters
----------
$SCORE_PARAMS
"""
@detector mutable struct DeterministicTransformer <: MLJ.Static
    normalize::Function = scale_minmax
    combine::Function = combine_mean
    classify::Function = default_percentile_threshold
end

const StaticTransformer = Union{
    ScoreTransformer,
    DeterministicTransformer,
    ProbabilisticTransformer
}

# returns the augmented train/test scores
function MLJ.transform(ev::StaticTransformer, _, scores::Tuple{Scores, Scores}...) # _ because there is no fitresult
    _, scores_test = to_scores(ev.normalize, ev.combine, scores...)
    scores_test
end

function MLJ.predict(ev::ProbabilisticTransformer, _, scores::Tuple{Scores, Scores}...) # _ because there is no fitresult
    _, scores_test = to_scores(ev.normalize, ev.combine, scores...)
    to_univariate_finite(scores_test)
end

function MLJ.predict(ev::DeterministicTransformer, _, scores::Tuple{Scores, Scores}...) # _ because there is no fitresult
    scores = to_scores(ev.normalize, ev.combine, scores...)
    _, classes_test = ev.classify(scores)
    to_categorical(classes_test)
end

function to_scores(normalize::Function, combine::Function, scores::Tuple{Scores, Scores}...)::Tuple{Scores, Scores}
    # make sure that we have at least one tuple of train and test scores
    @assert length(scores) > 0

    # [(train1, test1), (train2, test2)] -> matrix of scores for train [train1 train2]' and test [test1 test2]'
    # Note: The matrices contain one observation per row!
    scores_train, scores_test = combine(map(normalize, scores)...)
    scores_train, scores_test
end


text:
# these graph layout functions are experimental and mostly for decorative purposes
# adapted from various places, modified for use with Luxor

"""
    layoutgraph(am::Array{T,2}; kwargs...)

Layout the graph in adjacency matrix `am`, using a spring-based method. Returns a tuple of `([xcoords], [ycoords])`, which are scaled according to the supplied BoundingBox.

### Keyword arguments

Lay out the graph in `adjmatrix` using the spring/repulsion model of Fruchterman and Reingold (1991). Returns a tuple of `([xcoords], [ycoords])`.

## Arguments

`am`
    An adjacency matrix of some type. Non-zero of the eltype of the matrix is used to determine if a link exists, but currently no sense of magnitude

`densityconstant`
    Constant to adjust the density of resulting layout. The default value is 2.0.

`maxiterations`
    how many iterations for applying the forces

`initialtemperature`
    the initial temperature controls movement per iteration. The default value is 2.0. Each iteration uses a lower temperature (TEMP = initialtemperature / iter). The idea is that the displacements of vertices are limited to some maximum temperature value, and this decreases over time. As the layout becomes better, the amount of adjustment becomes smaller.

`locs_x`
    the starting x values, between -1 and 1. The default values are randomly selected.

`locs_y`
    the starting y values, between -1 and 1. The default values are randomly selected.

`boundingbox`
    The Luxor BoundingBox into which the graph's coordinates will fit. Defaults to 500 × 500.

"""
layoutgraph(am::Array{T,2} where T; kwargs...) = layout_spring(am; kwargs...)

"""
    layoutgraph(al::AbstractVector; kwargs)

Layout the graph in adjacency list `al`, assuming the graph is a tree layout with a single root node. Return a list of Luxor points that can be used to lay out the nodes of the tree.

### Keyword arguments

```nodesize=ones(size(adjlist, 1)))```

`nodesize` is a vector of values that define the area occupied by each node. No idea how it works yet...

Return a list of Luxor points that can be used to lay out a tree defined in adjacency list `adjlist`.
"""
layoutgraph(al::AbstractVector; kwargs...) = layout_tree(al; kwargs...)

"""
    layout_spring(adjmatrix::Array{T,2} where T;
        densityconstant = 2.0,
        maxiterations = 100,
        initialtemperature = 2.0,
        boundingbox=BoundingBox(O - (250, 250), O + (250, 250)))

"""
function layout_spring(adjmatrix::Array{T,2} where T;
    densityconstant = 2.0,
    maxiterations = 100,
    initialtemperature = 2.0,
    locs_x = 2rand(size(adjmatrix, 1)) .- 1.0,
    locs_y = 2rand(size(adjmatrix, 1)) .- 1.0,
    boundingbox =  BoundingBox(O - (250, 250), O + (250, 250)))

    N = size(adjmatrix, 1)
    if N != size(adjmatrix, 2)
        error("Adjacency matrix must be square but is $(size(adjmatrix)).")
    end

    # randomize results arrays, normalized to between -1 and 1

    # the optimal distance bewteen normalized vertices
    K = densityconstant * sqrt(4.0 / N)

    # store forces and apply at end of iteration all at once
    force_x = zeros(N)
    force_y = zeros(N)

    @inbounds for iter in 1:maxiterations
        # Calculate forces
        for i in 1:N
            force_vec_x = 0.0
            force_vec_y = 0.0
            for j in 1:N
                i == j && continue
                d_x = locs_x[j] - locs_x[i]
                d_y = locs_y[j] - locs_y[i]
                d   = sqrt(d_x^2 + d_y^2)
                if adjmatrix[i,j] != zero(eltype(adjmatrix)) || adjmatrix[j,i] != zero(eltype(adjmatrix))
                    F_d = d^2 / K
                else
                    # Just repulsive
                    F_d = -K^2 / d
                end
                force_vec_x += F_d * d_x
                force_vec_y += F_d * d_y
            end
            force_x[i] = force_vec_x
            force_y[i] = force_vec_y
        end
        # Cool down
        TEMP = initialtemperature / iter
        # apply but limit to temperature
        for i = 1:N
            force_mag  = sqrt(force_x[i] ^ 2 + force_y[i] ^ 2)
            scale      = min(force_mag, TEMP)/force_mag
            locs_x[i] += force_x[i] * scale
            locs_y[i] += force_y[i] * scale
        end
    end

    # fit to bounding box argument
    min_x, max_x = boxbottomleft(boundingbox).x, boxbottomright(boundingbox).x
    min_y, max_y = boxtopleft(boundingbox).y, boxbottomleft(boundingbox).y

    return (
        rescale.(locs_x, extrema(locs_x)..., min_x, max_x),
        rescale.(locs_y, extrema(locs_y)..., min_y, max_y)
        )
end

# layout algorithm for a simple tree structure

struct TreeLayout{A<:AbstractVector, F}
   nodes::A # []
   mod::F # []
   thread::Vector{Int} # []
   ancestor::Vector{Int} # []
   prelim::F # []
   shift::F # []
   change::F # []
   positions # Point[]
   nodesize::F # []
end

function TreeLayout(tree::AbstractVector, nodesize)
   len = length(tree)
   mod = zeros(len)
   thread = zeros(Int, len)
   prelim = zeros(len)
   shift = zeros(len)
   change = zeros(len)
   ancestor = collect(1:len)
   nodes = copy(tree)
   positions = [Point(0, 0) for i in 1:len]
   t = TreeLayout(nodes, mod, thread, ancestor, prelim, shift, change, positions, nodesize)
   return t
end

"""
   layout_tree(adjlist;
       nodesize=ones(size(am, 1)))
"""
function layout_tree(adjlist::AbstractVector;
   nodesize=ones(size(adjlist, 1)))
   # can check for stuff here
   _layout_tree(adjlist, nodesize)
end

# initialize modifiers, threads, and ancestors, then start the bottom-up
# and top-down traversals of the tree

function _layout_tree(t::AbstractVector, nodesize)
   tree = TreeLayout(t, nodesize)
   _layout_first_walk(1, tree)
   _layout_second_walk(1, -tree.prelim[1], 0.0, tree)
   return tree.positions
end

function _layout_parent(v, t::TreeLayout)
   tree = t.nodes
   for i in 1:length(tree)
       y = findfirst(x -> (x==v), tree[i])
       if !isnothing(y)
           return i
       end
   end
   return nothing
end

# computes a preliminary x-coordinate for v. Before that, FIRSTWALK is applied recursively
# to the children of v, as well as the function APPORTION. After spacing out the children by
# calling EXECUTESHIFTS, the node v is placed to the midpoint of its outermost children.

function _layout_first_walk(v, t::TreeLayout)
   prelim = t.prelim
   mod = t.mod
   tree = t.nodes
   nodesize = t.nodesize
   p = _layout_parent(v,t)
   if p != nothing
       index = findfirst(x -> (x==v), tree[p])
   else
       index = 1
   end
   if length(tree[v]) == 0
       if v != tree[p][1]
         prelim[v] = prelim[tree[p][index-1]] + (nodesize[tree[p][index-1]])
       else
         prelim[v] = 0
       end
   else
       defaultAncestor = tree[v][1]
       for w in tree[v]
           _layout_first_walk(w, t)
           defaultAncestor = _layout_apportion(w, defaultAncestor, t)
       end
       _layout_execute_shifts(v, t)
       midpoint = (prelim[tree[v][1]] + prelim[tree[v][end]]) / 2
       if index > 1
           w = tree[p][index-1]
           prelim[v] = prelim[w] + (nodesize[w] + 1.0)
           mod[v] = prelim[v] - midpoint
       else
           prelim[v] = midpoint
       end
   end
end

# a new subtree is combined with the previous subtrees
function _layout_apportion(v::T, defaultAncestor::T, t::TreeLayout) where {T}
   tree = t.nodes
   ancestor = t.ancestor
   prelim = t.prelim
   mod = t.mod
   thread = t.thread
   p = _layout_parent(v, t)
   nodesize = t.nodesize
   if p != nothing
       index = findfirst(x-> (x==v), tree[p])
   else
       index = 1
   end
   # "threads" are used to traverse the inside and outside contours of the left and right subtree up to the highest common level
   if index > 1
       w = tree[p][index-1]
       v_in_right = v_out_right = v
       v_in_left = w
       v_out_left = tree[_layout_parent(v_in_right, t)][1]
       s_in_right = mod[v_in_right]
       s_out_right = mod[v_out_right]
       s_in_left = mod[v_in_left]
       s_out_left = mod[v_out_left]
       while _layout_next_right(v_in_left, t) !=0 && _layout_next_left(v_in_right, t) !=0
           v_in_left = _layout_next_right(v_in_left, t)
           v_in_right = _layout_next_left(v_in_right, t)
           v_out_left = _layout_next_left(v_out_left, t)
           v_out_right = _layout_next_right(v_out_right, t)
           ancestor[v_out_right] = v
           shift = (prelim[v_in_left] + s_in_left) - (prelim[v_in_right] + s_in_right) + (nodesize[v_in_left])
           if shift > 0
               _layout_move_subtree(_layout_find_ancestor(v_in_left, v, defaultAncestor, t), v, shift, t)
               s_in_right += shift
               s_out_right += shift
           end
           s_in_left +=   mod[v_in_left]
           s_in_right +=  mod[v_in_right]
           s_out_left +=  mod[v_out_left]
           s_out_right += mod[v_out_right]
       end
       if _layout_next_right(v_in_left, t) != 0 && _layout_next_right(v_out_right, t) == 0
           thread[v_out_right] = _layout_next_right(v_in_left, t)
           mod[v_out_right] += s_in_left - s_out_right
       else
           if _layout_next_left(v_in_right,t) != 0 && _layout_next_left(v_out_left,t) == 0
               thread[v_out_left] = _layout_next_left(v_in_right, t)
               mod[v_out_left] += s_in_right - s_out_left
               defaultAncestor = v
           end
       end
   end
   return defaultAncestor
end

function _layout_number(v, t::TreeLayout)
   p = _layout_parent(v, t)
   index = findfirst(x -> (x==v), t.nodes[p])
   return index
end

# Shifting a subtree can be done in linear time..
# Calling MOVESUBTREE(w_ ,w+, shift) first shifts the current subtree, rooted at w+. This
# is done by increasing prelim(w+) and mod(w+) by shift.

function _layout_move_subtree(w_left::T, w_right::T, shift::Float64, t::TreeLayout) where {T}
   change = t.change
   prelim = t.prelim
   tree = t.nodes
   shifttree = t.shift
   mod = t.mod
   n_wl = _layout_number(w_left, t)
   n_wr = _layout_number(w_right, t)
   subtrees = n_wr - n_wl
   change[w_right] -= shift / subtrees
   shifttree[w_right] += shift
   change[w_left] += shift / subtrees
   prelim[w_right] += shift
   mod[w_right] += shift
end

# SECONDWALK is used to compute all real x-coordinates by
# summing up the modifiers recursively

function _layout_second_walk(v, m::Float64, depth::Float64, t::TreeLayout)
   prelim = t.prelim
   mod = t.mod
   positions = t.positions
   nodesize = t.nodesize
   positions[v] = Point(prelim[v] + m, depth)
   if length(t.nodes[v]) != 0
       maxdist = maximum([nodesize[i] for i in t.nodes[v]])
   else
       maxdist = 0
   end
   for w in t.nodes[v]
       _layout_second_walk(w, m + mod[v], Float64(depth + 1 + maxdist), t)
   end
end

# ANCESTOR returns the left one of the greatest uncommon ancestors of
# vi- and its right neighbour

function _layout_find_ancestor(w::T, v::T, defaultAncestor::T, tree::TreeLayout) where {T}
   ancestor = tree.ancestor
   if ancestor[w] in tree.nodes[_layout_parent(v, tree)]
       return ancestor[w]
   else
       return defaultAncestor
   end
end

# EXECUTESHIFTS(v) only needs one traversal of the children of v to execute
# all shifts computed and memorized in MOVESUBTREE

function _layout_execute_shifts(v, t::TreeLayout)
   tree = t.nodes
   shift = t.shift
   change = t.change
   prelim = t.prelim
   mod = t.mod
   shiftnode = 0
   changenode = 0
   for w in reverse(tree[v])
       prelim[w] += shiftnode
       mod[w] += shiftnode
       changenode += change[w]
       shiftnode += shift[w] + changenode
   end
end

function _layout_next_left(v, t::TreeLayout)
   tree = t.nodes
   thread = t.thread
   if length(tree[v]) != 0
       return tree[v][1]
   else
       return thread[v]
   end
end

function _layout_next_right(v, t::TreeLayout)
   tree = t.nodes
   thread = t.thread
   if length(tree[v]) != 0
       return tree[v][end]
   else
       return thread[v]
   end
end


text:
"""
    abstract type ImageQualityIndex <: Any

Root of `ImageQualityIndexes` type hierarchy.

Classified by whether a reference image is needed, there are three more abstract types defined: [`FullReferenceIQI`](@ref), [`ReducedReferenceIQI`](@ref) and [`NoReferenceIQI`](@ref)

!!! info

    Image quality index is also known as image quality assessment, image quality metrics and image quality methods.

# References

[1] Wang, Z., & Bovik, A. C. (2006). _Modern Image Quality Assessment_. Morgan & Claypool Publishers.
"""
abstract type ImageQualityIndex end

"""
    abstract type FullReferenceIQI <: ImageQualityIndex
    assess(::FullReferenceIQI, img, ref_img, args...)

Image quality index that requires its reference to be complete/full image or image series.

See also: [`ImageQualityIndex`](@ref), [`ReducedReferenceIQI`](@ref), [`NoReferenceIQI`](@ref)
"""
abstract type FullReferenceIQI <: ImageQualityIndex end

"""
    abstract type ReducedReferenceIQI <: ImageQualityIndex
    assess(::ReducedReferenceIQI, img, ref_info, args...)

Image quality index that requires some information from reference image.

Different from [`FullReferenceIQI`](@ref), `ReducedReferenceIQI` doesn't need the complete reference image.

See also: [`ImageQualityIndex`](@ref), [`FullReferenceIQI`](@ref), [`NoReferenceIQI`](@ref)
"""
abstract type ReducedReferenceIQI <: ImageQualityIndex end

"""
    abstract type NoReferenceIQI <: ImageQualityIndex
    assess(::NoReferenceIQI, img, args...)

Image quality index that doesn't require any information.

See also: [`ImageQualityIndex`](@ref), [`FullReferenceIQI`](@ref), [`ReducedReferenceIQI`](@ref)
"""
abstract type NoReferenceIQI <: ImageQualityIndex end

@doc (@doc FullReferenceIQI)
assess(iqi::FullReferenceIQI, x, ref, args...) = iqi(x, ref, args...)

@doc (@doc ReducedReferenceIQI)
assess(iqi::ReducedReferenceIQI, x, ref, args...) = iqi(x, ref, args...)

@doc (@doc NoReferenceIQI)
assess(iqi::NoReferenceIQI, x, args...) = iqi(x, args...)


text:
using Plots, LaTeXStrings, Statistics, StatsPlots, JLD, ProgressBars, StatsBase, Ranges

function FindEnergy(System::Matrix{Int8}, dim::Integer)
    E = 0.0
    for i ∈ 1:dim
        for j ∈ 1:dim
            Neighbors = [[1 + i % dim, j], [i, 1 + j % dim], [(-2 + i + dim) % dim + 1, j], [i, (-2 + j + dim) % dim + 1]]
            E += -sum(System[i, j] * System[Neighbors]) / 2
        end
    end
    return E / dim^2
end

function DE(System::Matrix{Int8}, i::Integer, j::Integer)
    Neighbors = [[1 + i % dim, j], [i, 1 + j % dim], [(-2 + i + dim) % dim + 1, j], [i, (-2 + j + dim) % dim + 1]]
    return system[i, j] * sum(System[CartesianIndex.(Tuple.(Neighbors))])
end

function MC(System::Matrix{Int8}, dim::Integer, N::Integer, β::Real)
    ΔE_sum = 0
    ΔM_sum = 0
    for _ = 1:N
        i = rand(1:dim)
        j = rand(1:dim)
        System[i, j] = spin
        ΔE = DE(System, i, j)
        if ΔE <= 0 || rand() < exp(-ΔE * β)
            System[i, j] = -spin
            ΔE_sum += ΔE
            ΔM_sum += 2 * (-spin)
        end
    end
    return System, ΔE_sum, ΔM_sum
end

function ising2d(dim::Integer, β::Real, MCS::Real, N::Integer)
    System = rand([Int8(-1), Int8(1)], dim, dim)

    Eᵢ = FindEnergy(system, dim)
    Mᵢ = mean(System)
    M_list = [Eᵢ]
    E_list = [Mᵢ]
    for ـ = 1:MCS
        System, ΔE, ΔM = MC(System, dim, N, β)
        E = E_list[end] + ΔE
        M = M_list[end] + ΔM
        push!(E_list, E)
        push!(M_list, M)
    end
    return mean(E_list), mean(M_list), β^2 * var(E_list[ceil(Int64, MCS / 5):end]) * dim^2, β * var(E_list[ceil(Int64, MCS / 5):end])
end

#Initials
dim = 4
β_list = range(0.1, 0.9, 80)

E_list = []
M_list = []
C_list = []
χ_list = []

for β in β_list
    E, M, C, χ = ising2d(dim, β, 1000, 100)
    push!(E_list, E)
    push!(M_list, M)
    push!(C_list, C)
    push!(χ_list, χ)
end

plot(β_list, E_list)

plot(β_list, M_list)

plot(β_list, C_list)

plot(β_list, χ_list)

E_Data = []
M_Data = []
C_Data = []
χ_Data = []

for l in 3:12
    E_list = []
    M_list = []
    C_list = []
    χ_list = []

    for β in β_list
        E, M, C, χ = ising2d(dim, β, 1000, 100)
        push!(E_list, E)
        push!(M_list, M)
        push!(C_list, C)
        push!(χ_list, χ)
    end
    push!(E_Data, E_list)
    push!(M_Data, M_list)
    push!(C_Data, C_list)
    push!(χ_Data, χ_list)
end


text:
using Statistics: quantile
using Interpolations: LinearInterpolation

normal_pdf(Δ, σ) = exp(-0.5*Δ^2 / σ^2) / √(2π) / σ

"""
    constant_R_LSF(flux, wls, R)

Applies a gaussian line spread function the the spectrum with flux vector `flux` and wavelength
vector `wls` with constant spectral resolution, ``R = \\lambda/\\Delta\\lambda.``

!!! warning
    This is a naive, slow implementation.  Do not use it when performance matters.

    `constant_R_LSF` will have weird behavior if your wavelength grid is not locally linearly-spaced.
    It is intended to be run on a fine wavelength grid, then downsampled to the observational (or
    otherwise desired) grid.
"""
function constant_R_LSF(flux::AbstractVector{F}, wls, R) where F <: Real
    #ideas - require wls to be a range object? Use erf to account for grid edges?
    convF = zeros(F, length(flux))
    for i in 1:length(wls)
        λ0 = wls[i]
        σ = λ0 / R / 2
        mask = λ0 - 4σ .< wls .< λ0 + 4σ
        ϕ = normal_pdf.(wls[mask] .- λ0, λ0 / R / 2)
        ϕ ./= sum(ϕ)
        convF[mask] += flux[i]*ϕ
    end
    convF
end

"""
    rectify(flux, wls; bandwidth=50, q=0.95, wl_step=1.0)

Rectify the spectrum with flux vector `flux` and wavelengths `wls` by dividing out a moving
`q`-quantile with window size `bandwidth`.  `wl_step` controls the size of the grid that the moving
quantile is calculated on and interpolated from.  Setting `wl_step` to 0 results in the exact
calculation with no interpolation, but note that this is very slow.

Experiments on real spectra show an agreement between the interpolated rectified spectrum and the
"exact" one (with default values) at the 3 × 10^-4 level.
"""
function rectify(flux::AbstractVector{F}, wls; bandwidth=50, q=0.95, wl_step=1.0) where F <: Real
    inds = 1 : max(1, Int(floor(wl_step/step(wls)))) : length(wls)
    lb = 1
    ub = 1
    moving_quantile = map(wls[inds]) do λ
        #move_bounds is defined in line_opacity.jl
        lb, ub = move_bounds(wls, lb, ub, λ, bandwidth)
        quantile(flux[lb:ub], q)
    end
    itp = LinearInterpolation(wls[inds], moving_quantile)
    flux ./ itp.(wls)
end

"""
    air_to_vacuum(λ; cgs=λ<1)

Convert λ from an air wavelength to a vacuum wavelength.
λ is assumed to be in Å if it is ⩾ 1, in cm otherwise.
Formula from Birch and Downs (1994) via the VALD website.
"""
function air_to_vacuum(λ; cgs=λ<1)
    if cgs
        λ *= 1e8
    end
    s = 1e4/λ
    n = 1 + 0.00008336624212083 + 0.02408926869968 / (130.1065924522 - s^2) + 0.0001599740894897 / (38.92568793293 - s^2)
    if cgs
        λ *= 1e-8
    end
    λ * n
end

"""
    vacuum_to_air(λ; cgs=λ<1)

convert λ from a vacuum wavelength to an air wavelength
λ is assumed to be in Å if it is ⩾ 1, in cm otherwise.
Formula from Birch and Downs (1994) via the VALD website.
"""
function vacuum_to_air(λ; cgs=λ<1)
    if cgs
        λ *= 1e8
    end
    s = 1e4/λ
    n = 1 + 0.0000834254 + 0.02406147 / (130 - s^2) + 0.00015998 / (38.9 - s^2)
    if cgs
        λ *= 1e-8
    end
    λ / n
end

"""
    parse_fwf(source, rowspec; datarow, lastrow)

Parse fixed-width files, returning a vector of `NamedTuples`.
`source` can be either an `AbstractVector` of `Strings` or a filename.
`rowspec` should be a vector of tuples each containing, in order
- a range object corresponding to the column postion
- the type.  Strings will be pulled directly & other types will be parsed with `parse`.
- the column name, a `Symbol`.
- optionally, a function to apply to the parsed value.
"""
function parse_fwf(filename::String, rowspec; datarow=1, lastrow=0)
    open(filename, "r") do file
        for _ in 1:(datarow-1)
            readline(file)
        end
        _parse_fwf_core(eachline(file), rowspec, lastrow-datarow+1)
    end
end
function parse_fwf(lines::AbstractVector{String}, rowspec; datarow=1, lastrow=0)
    lastrow = lastrow == 0 ? length(lines) : lastrow
    _parse_fwf_core(lines[datarow:lastrow], rowspec)
end

function _parse_fwf_core(lines, rowspec, n=0)
    #use n
    map(lines) do line
        #broadcasts over each itemspec and passes to NamedTuple constuctor as kwargs
        (; _parse_item.(line, rowspec)...)
    end
end

function _parse_item(line, itemspec)
    if length(itemspec) == 4
        r, t, n, f = itemspec
    else
        r, t, n = itemspec
        f = identity
    end
    if t == String
        (n, f(line[r]))
    elseif r[2] > length(line) || strip(line[r]) == ""
        (n, zero(t))
    else
        (n, f(parse(t, line[r])))
    end
end

# take some care to avoid subnormal values (they can slow things down a lot!)
_nextfloat_skipsubnorm(v::F) where {F<:AbstractFloat} =
    ifelse(-floatmin(F) ≤ v < 0, F(0), ifelse(0 ≤ v < floatmin(F), floatmin(F), nextfloat(v)))
_prevfloat_skipsubnorm(v::F) where {F<:AbstractFloat} =
    ifelse(-floatmin(F) < v ≤ 0, -floatmin(F), ifelse(0 < v ≤ floatmin(F), F(0), prevfloat(v)))

struct Interval # represents an exclusive interval
    lower::Float64
    upper::Float64

    function Interval(lower::Real, upper::Real; exclusive_lower = true, exclusive_upper = true)
        @assert lower < upper "the upper bound must exceed the lower bound"
        lower, upper = Float64(lower), Float64(upper)
        new((exclusive_lower || isinf(lower)) ? lower : _prevfloat_skipsubnorm(lower),
            (exclusive_upper || isinf(upper)) ? upper : _nextfloat_skipsubnorm(upper))
    end
end

# convenience function for defining interval where both bounds are inclusive
closed_interval(lo, up) = Interval(lo, up; exclusive_lower = false, exclusive_upper = false)

"""
    contained(value, interval)

Returns whether `value` is contained by `interval`.

# Examples
```julia-repl
julia> contained(0.5, Interval(1.0,10.0))
false
julia> contained(5.0, Interval(1.0,10.0))
true
```
"""
contained(value::Real, interval::Interval) = interval.lower < value < interval.upper

"""
    contained_slice(vals, interval)

Returns a range of indices denoting the elements of `vals` (which are assumed to be sorted in
 increasing order) that are contained by `interval`. When no entries are contained by interval,
this returns `(1,0)` (which is a valid empty slice).
"""
contained_slice(vals::AbstractVector, interval::Interval) =
    searchsortedfirst(vals, interval.lower):searchsortedlast(vals, interval.upper)

function _convert_λ_endpoint(λ_endpoint::AbstractFloat, λ_lower_bound::Bool)
    # determine the functions that:
    # - retrieve the neighboring ν val in the in-bounds and out-of-bounds directions
    # - specify the desired relationship between an in-bounds λ and λ_lower_bound
    inbound_ν_neighbor, oobound_ν_neighbor, inbound_λ_to_endpoint_relation =
        λ_lower_bound ? (prevfloat, nextfloat, >) : (nextfloat, prevfloat, <)

    ν_endpoint = (λ_endpoint == 0) ? Inf : c_cgs/λ_endpoint
    if isfinite(ν_endpoint) && (ν_endpoint != 0)
        # Adjust ν_endpoint in 2 cases:
        # Case 1: The neighboring ν to ν_endpoint that should be in-bounds, corresponds to a λ
        # value that is out-of-bounds. Nudge ν_endpoint in the "in-bounds" direction until resolved
        while !inbound_λ_to_endpoint_relation(c_cgs/inbound_ν_neighbor(ν_endpoint), λ_endpoint)
            ν_endpoint = inbound_ν_neighbor(ν_endpoint)
        end
        # Case 2: The current value of ν_endpoint corresponds to a λ that is in-bounds. Nudge
        # ν_endpoint in the "out-of-bounds" direction until resolved.
        while inbound_λ_to_endpoint_relation(c_cgs/ν_endpoint, λ_endpoint)
            ν_endpoint = oobound_ν_neighbor(ν_endpoint)
        end
    end
    ν_endpoint
end

"""
    λ_to_ν_bound(λ_bound)

Converts a λ `Inverval` (in cm) to an equivalent ν `Interval` (in Hz), correctly accounting for
tricky floating point details at the bounds.
"""
λ_to_ν_bound(λ_bound::Interval) =
    Interval(_convert_λ_endpoint(λ_bound.upper, false), _convert_λ_endpoint(λ_bound.lower, true))


text:
module Docker

using AWS
using ..ECR

export login, pull, push, build

"""
    login(registry_id::Union{AbstractString, Integer}="")

Gets the AWS ECR authorization token and runs the corresponding docker login command.
Takes an optional AWS `config` keyword argument or uses the default.
"""
function login(
    registry_id::Union{AbstractString,Integer}=""; config::AWSConfig=global_aws_config()
)
    return login(ECR.get_login(registry_id; config=config))
end

function login(docker_login::Cmd)
    success(pipeline(docker_login, stdout=stdout, stderr=stderr))
end

"""
    pull(image::AbstractString, tags::AbstractVector{<:AbstractString}=String[])

Pulls a docker image and tags it if `tags` is specified.
"""
function pull(image::AbstractString, tags::AbstractVector{<:AbstractString}=String[])
    run(`docker pull $image`)
    for tag in tags
        run(`docker tag $image $tag`)
    end
end

"""
    push(image::AbstractString)

Pushes a docker image.
"""
function push(image::AbstractString)
    run(`docker push $image`)
end

"""
    build(dir::AbstractString, tag::AbstractString="")

Builds the docker image.
"""
function build(dir::AbstractString, tag::AbstractString="")
    opts = isempty(tag) ? `` : `-t $tag`
    run(`docker build $opts $dir`)
end

end  # Docker


text:
module OperatorLearning

using Base: Integer, ident_cmp, Float32
using CUDA
using Flux
using FFTW
using FFTW: assert_applicable, unsafe_execute!, FORWARD, BACKWARD, rFFTWPlan
using Random
using Random: AbstractRNG
using Flux: nfan, glorot_uniform, batch
using OMEinsum

export FourierLayer, DeepONet

include("FourierLayer.jl")
include("DeepONet.jl")
include("ComplexWeights.jl")
include("batched.jl")
include("subnets.jl")

end # module


text:
### TODO: rewrite this type similarly to FinancialTimeSeries

type TickData{T<:Float64,N}  <: AbstractTimeSeries

    timestamp::Vector{Date}
    values::Array{T,N}
    colnames::Vector{ASCIIString}
    instrument::AbstractInstrument

    function TickData(timestamp::Vector{Date},
                      values::Array{T,N},
                      colnames::Vector{ASCIIString},
                      instrument::AbstractInstrument)

                      v = findfirst(colnames .== "Volume") # first sort the values in ascending order
                      #vals = sortrows(values, by=x->x[v])
                      if v > 0
                          vals = sortrows(values, by=x->x[v])
                      else
                          vals = values # don't sort if no Valume column available
                      end
                      nrow, ncol = size(values, 1), size(values, 2)
                      nrow != size(timestamp, 1) ? error("values must match length of timestamp"):
                      ncol != size(colnames,1) ? error("column names must match width of array"):
                      ~(flipud(timestamp) == sort(timestamp) || timestamp == sort(timestamp)) ? error("dates are mangled"):
                      flipud(timestamp) == sort(timestamp) ?
                      new(flipud(timestamp), vals, colnames, instrument): # then the dates if necessary
                      new(timestamp, vals, colnames, instrument)
    end
end

TickData{T<:Float64,N}(d::Vector{Date}, v::Array{T,N}, c::Vector{ASCIIString}, t::AbstractInstrument) = TickData{T,N}(d,v,c,t)
TickData{T<:Float64,N}(d::Date, v::Array{T,N}, c::Array{ASCIIString,1}, t::AbstractInstrument) = TickData([d],v,c,t)

# ########################### TickData
#
# function show(io::IO, ft::TickData)
#   # variables
#   nrow          = size(ft.values, 1)
#   ncol          = size(ft.values, 2)
#   intcatcher    = falses(ncol)
#   for c in 1:ncol
#       rowcheck =  trunc(ft.values[:,c]) - ft.values[:,c] .== 0
#       if sum(rowcheck) == length(rowcheck)
#           intcatcher[c] = true
#       end
#   end
#   spacetime     = strwidth(string(ft.timestamp[1])) + 3
#   firstcolwidth = strwidth(ft.colnames[1])
#   colwidth      = Int[]
#       for m in 1:ncol
#           push!(colwidth, max(strwidth(ft.colnames[m]), strwidth(@sprintf("%.2f", maximum(ft.values[:,m])))))
#       end
#
#   # summary line
#   print(io,@sprintf("%dx%d %s %s to %s", nrow, ncol, typeof(ft), string(ft.timestamp[1]), string(ft.timestamp[nrow])))
#   println(io,"")
#   println(io,"")
#
#   # row label line
#    print(io, ^(" ", spacetime), ft.colnames[1], ^(" ", colwidth[1] + 2 -firstcolwidth))
#
#    for p in 2:length(colwidth)
#      print(io, ft.colnames[p], ^(" ", colwidth[p] - strwidth(ft.colnames[p]) + 2))
#    end
#    println(io,"")
#
#   # timestamp and values line
#     if nrow > 7
#         for i in 1:4
#             print(io, ft.timestamp[i], " | ")
#         for j in 1:ncol
#             intcatcher[j] ?
#             print(io, rpad(iround(ft.values[i,j]), colwidth[j] + 2, " ")) :
#             print(io, rpad(round(ft.values[i,j], 2), colwidth[j] + 2, " "))
#         end
#         println(io,"")
#         end
#         println(io,'\u22EE')
#         for i in nrow-3:nrow
#             print(io, ft.timestamp[i], " | ")
#         for j in 1:ncol
#             intcatcher[j] ?
#             print(io, rpad(iround(ft.values[i,j]), colwidth[j] + 2, " ")) :
#             print(io, rpad(round(ft.values[i,j], 2), colwidth[j] + 2, " "))
#         end
#         println(io,"")
#         end
#     else
#         for i in 1:nrow
#             print(io, ft.timestamp[i], " | ")
#         for j in 1:ncol
#             intcatcher[j] ?
#             print(io, rpad(iround(ft.values[i,j]), colwidth[j] + 2, " ")) :
#             print(io, rpad(round(ft.values[i,j], 2), colwidth[j] + 2, " "))
#         end
#         println(io,"")
#         end
#     end
# end


text:
module VimBindings

using REPL: history_next, history_prev
using Base: AnyDict
using REPL
using REPL.LineEdit
import REPL.LineEdit: KeyAlias, edit_splice!, buffer, refresh_line
import Base: AnyDict, show_unquoted
using Sockets

const LE = LineEdit

include("util.jl")
include("textutils.jl")
include("types.jl")
include("command.jl")
include("textobject.jl")
include("execute.jl")
include("motion.jl")
include("operator.jl")
include("keys.jl")
include("parse.jl")

using .Parse
using .Commands
mutable struct VimState
    registers :: Dict{Char, String}
    register :: Char
end

VimState() = VimState(Dict{Char, String}(), '"')

const vim = VimState()
"""
dispatch_key method which takes in the char, VB.mode, s.
Same dispatch_key method for all keys.
It will do the eval of the symbol at runtime.

e.g. for normal mode (union of many modes, possible), it will dispatch_key
to the specific keybinds

but for "find" mode, there is a dispatch_key method which uses the character
to find, instead of calling that character's function.

for normal mode, take the @eval code below in order to call
the correct function at runtime.

This way, we do not need function specific bindings for every single
key, only the needed (and/or implemented) ones.
"""
function dispatch_key(c, mode, s::LE.MIState)
    log("dispatching key: ", string(c))
    fn_name = if c in keys(special_keys)
        Symbol(special_keys[c])
    elseif alphabetic(c)
        if is_lowercase(c)
            Symbol(c)
        else
            Symbol(string(lowercase(c), "_uppercase"))
        end
    else
        log("not dispatching key")
        @log c
        return
    end
    if !isdefined(VimBindings, fn_name)
        log("function $fn_name does not exist")
        return
    end
    eval(Expr(:call, fn_name, vim.mode, s))
end

# function dispatch_key(query_c :: Char, mode :: FindChar{}, s::LE.MIState)
#     buf = buffer(s)
#     motion = find_c(buf, query_c)
#     execute(mode, s, motion)
#     LE.refresh_line(s)
#     vim_reset()
# end

# function dispatch_key(c, mode :: ToChar, s::LE.MIState)
#     # TODO
# end

# function dispatch_key(c :: Char, mode :: SelectRegister, s::LE.MIState)
#     if alphanumeric(c) || c == '_'
#         @log vim.register = c
#     end
#     vim_reset()
# end

global key_stack = Char[]

function strike_key(c, s::LE.MIState)
    # if(c == '`')
    #     empty!(key_stack)
    #     @log key_stack
    # else
    append!(key_stack, c)
    s_cmd = String(key_stack)
    if well_formed(s_cmd)
        log("well formed command: $s_cmd")
        empty!(key_stack)
        cmd = parse_command(s_cmd)
        if cmd !== nothing
            refresh :: Bool = execute(s, cmd)
            if refresh
                LE.refresh_line(s)
            end
        end
    else
        @log key_stack
    end
    # end
end

function init()
    repl = Base.active_repl
    global juliamode = repl.interface.modes[1]
    historymode = repl.interface.modes[4]
    prefixhistorymode = repl.interface.modes[5]
    juliamode.prompt = "julia[i]> "
    # juliamode.keymap_dict['`'] = trigger_normal_mode
    LE.add_nested_key!(juliamode.keymap_dict, "\e\e", trigger_normal_mode)
    # LE.add_nested_key!(historymode.keymap_dict, "\e\e", trigger_normal_mode)
    LE.add_nested_key!(prefixhistorymode.keymap_dict, "\e\e", trigger_normal_mode)
    # julia_mode_new_keys = AnyDict(
    #     "\e\e" =>  trigger_normal_mode,
    #     '`' =>  trigger_normal_mode
    # ) |> LE.normalize_key
    # juliamode.keymap_dict = merge(juliamode.keymap_dict, julia_mode_new_keys)
    # juliamode.keymap_dict["\e\e"] = trigger_normal_mode

    # remove normal mode if it's already added
    normalindex = 0
    for (i, m) in enumerate(repl.interface.modes)
        if hasproperty(m, :prompt) && m.prompt == "julia[n]> "
            normalindex = i
        end
    end
    if normalindex != 0
        deleteat!(repl.interface.modes, normalindex)
    end

    binds = AnyDict()
    for c in all_keys
        bind = (s::LE.MIState, o...)->begin
            strike_key(c, s)
            # dispatch_key(c, vim.mode, s)
        end
        binds[c] = bind
    end

    keymap = AnyDict(
        '*' => (s::LE.MIState, o...)->begin
            log("keymap fallthrough: *")
            # @log o
        end,
        # "\e" => (s, o...)->begin
        #     log("keymap: \\e, escape")
        # end,

        "\e\e" => (s, o...)->begin
            empty!(key_stack)
        end,
        # "\e[A" => (s::LE.MIState, o...)->begin
            # log("Up Arrow")
            # @log o
        # end,
        # backspace
        # '\b' => (s::LE.MIState, o...)->LE.edit_move_left(s),
             # '`' => (s::LE.MIState, o...)->vim_reset(),
    )
    keymap = merge(keymap,
                   AnyDict(binds))

    # keys to copy from `juliamode`
    copy_keys = [
        # enter
        '\r',
        # tab
        # '\t',
        # newline
        '\n',
        # home
        "\e[H",
        # end
        "\e[F",
        # clear
        "^L",
        # right arrow
        "\e[C",
        # left arrow
        "\e[D",
        # up arrow
        # "\e[A",
        # down arrow
        # "\e[B",
        # delete
        "\e[3~",
    ]

    for c in copy_keys
        keymap[c] = LE.default_keymap[c]
    end
    global normalmode =
        REPL.Prompt("julia[n]> ",
                    keymap_dict = LE.keymap([keymap]),
                    # on_done = REPL.respond(split, repl, juliamode),
                    on_done = (s, o...) -> begin
                        log("on_done")
                        @log s.last_action
                        @log s.current_action
                    end
                    # on_enter = juliamode.on_enter,
                    )
    normalmode.on_done = juliamode.on_done
    normalmode.on_enter = juliamode.on_enter
    normalmode.hist = juliamode.hist
    normalmode.hist.mode_mapping[:julia] = normalmode

    push!(repl.interface.modes, normalmode)
    return
end

import REPL.LineEdit: TextTerminal, ModalInterface, MIState, activate, keymap, match_input, keymap_data, transition, mode, terminal
import REPL.Terminals: raw!, enable_bracketed_paste, disable_bracketed_paste
function LE.prompt!(term::TextTerminal, prompt::ModalInterface, s::MIState = init_state(term, prompt))
    # println("initializing prompt from VimBindings.jl")
    Base.reseteof(term)
    raw!(term, true)
    enable_bracketed_paste(term)
    try
        activate(prompt, s, term, term)
        old_state = mode(s)
        while true
            kmap = keymap(s, prompt)
            fcn = match_input(kmap, s)
            kdata = keymap_data(s, prompt)
            s.current_action = :unknown # if the to-be-run action doesn't update this field,
                                        # :unknown will be recorded in the last_action field
            local status
            # errors in keymaps shouldn't cause the REPL to fail, so wrap in a
            # try/catch block
            try
                status = fcn(s, kdata)
            catch e
                @error "Error in the keymap" exception=e,catch_backtrace()
                # try to cleanup and get `s` back to its original state before returning
                transition(s, :reset)
                transition(s, old_state)
                status = :done
            end
            status !== :ignore && (s.last_action = s.current_action)
            if status === :abort
                s.aborted = true
                return buffer(s), false, false
            elseif status === :done
                return buffer(s), true, false
            elseif status === :suspend
                if Sys.isunix()
                    return buffer(s), true, true
                end
            else
                @assert status ∈ (:ok, :ignore)
            end
        end
    finally
        raw!(term, false) && disable_bracketed_paste(term)
    end
    # unreachable
end

# function LE.match_input(f::Function, s::Union{Nothing,LE.MIState}, term, cs::Vector{Char}, keymap)
#     log("match function")
#     LE.update_key_repeats(s, cs)
#     c = String(cs)
#     return function (s, p)  # s::Union{Nothing,MIState}; p can be (at least) a LineEditREPL, PrefixSearchState, Nothing
#         r = Base.invokelatest(f, s, p, c)
#         if isa(r, Symbol)
#             return r
#         else
#             return :ok
#         end
#     end
# end

# function LE.match_input(k::Nothing, s, term, cs, keymap)
#     log("match nothing")
#     # @log cs
#     return (s,p) -> begin
#         log("nothing")
#         return :ok
#     end
# end
# LE.match_input(k::KeyAlias, s::Union{Nothing,LE.MIState}, term, cs, keymap::Dict{Char}) = LE.match_input(keymap, s, IOBuffer(k.seq), Char[], keymap)

function LE.match_input(k::Dict{Char}, s::Union{Nothing,LE.MIState}, term::Union{LE.AbstractTerminal,IOBuffer}=terminal(s), cs::Vector{Char}=Char[], keymap::Dict{Char} = k)
    # if we run out of characters to match before resolving an action,
    # return an empty keymap function
    eof(term) && return (s, p) -> :abort
    c = read(term, Char)

    if isempty(cs) && c == '\e'
        is_escape_task = @async begin
            sleep(0.03)
            avail = bytesavailable(term)
            if avail > 0
                log("bytes available to read: suspected encoded sequence")
                return false
            else
                log("no bytes available to read: suspected Escape key")
                return true
            end
        end

        if fetch(is_escape_task)
            @log keys(k['\e'])
            # short-circuit completion here for Escape key
            if @log haskey(k['\e'], '\e')
                return LE.match_input(k['\e']['\e'], s, term, cs, nothing)
            end
        end
    end
    # Ignore any `wildcard` as this is used as a
    # placeholder for the wildcard (see normalize_key("*"))
    c == LE.wildcard && return (s, p) -> :ok
    push!(cs, c)
    key = haskey(k, c) ? c : LE.wildcard
    # if we don't match on the key, look for a default action then fallback on 'nothing' to ignore
    return LE.match_input(get(k, key, nothing), s, term, cs, keymap)
end

function vim_reset()
    vim.mode = NormalMode()
    return true
end

function edit_move_end(s::LE.MIState)
    buf = LE.buffer(s)
    @show typeof(buf)
    while !eof(buf)
        if linebreak(
            LE.char_move_right(buf))
            break
        end
    end
    pos = max(position(buf) - 1, # correcting adjustment
              0)
    seek(buf,pos)
    LE.refresh_line(s)
    return true
end

function edit_move_start(s::LE.MIState)
    buf = LE.buffer(s)
    while position(buf) > 0
        if linebreak(LE.char_move_left(buf))
            LE.char_move_right(buf)
            break
        end
    end
    pos = position(buf)
    seek(buf,pos)
    LE.refresh_line(s)
    return true
end

function edit_move_phrase_right(s::LE.MIState)
    buf = LE.buffer(s)
    if !eof(buf)
        LE.char_move_word_right(buf, whitespace)
        return LE.refresh_line(s)
    end
    return true
end

function edit_move_phrase_left(s::LE.MIState)
    buf = LE.buffer(s)
    if position(buf) > 0
        LE.char_move_word_left(buf, whitespace)
        return LE.refresh_line(s)
    end
    return true
end

function trigger_insert_mode(state::LineEdit.MIState)
    iobuffer = LineEdit.buffer(state)
    LineEdit.transition(state, juliamode) do
        prompt_state = LineEdit.state(state, juliamode)
        prompt_state.input_buffer = copy(iobuffer)
    end
end

function trigger_normal_mode(state::LineEdit.MIState, o...)
    iobuffer = LineEdit.buffer(state)
    # vim.mode = NormalMode()
    LineEdit.transition(state, normalmode) do
        prompt_state = LineEdit.state(state, normalmode)
        prompt_state.input_buffer = copy(iobuffer)
    end
end

function key_press(state::REPL.LineEdit.MIState, repl::LineEditREPL, char::String)
end

function getsocket()
    if !isdefined(VimBindings, :socket) || isa(socket, Base.DevNull)
        try
            global socket = connect(1234)
        catch e
            global socket = devnull
        end
    end
    return socket
end

function log(any::Any)
    socket = getsocket()
    println(socket, any)
end

function debug_mode(state::REPL.LineEdit.MIState, repl::LineEditREPL, char::String)
    socket = getsocket()
    # iobuffer is not displayed to the user
    iobuffer = LineEdit.buffer(state)

    # write character typed into line buffer
    #= LineEdit.edit_insert(iobuffer, char) =#

    # write character typed into repl
    #= LineEdit.edit_insert(state, char) =#
    #
    # move to start of iostream to read what the user
    # has typed
    seekstart(iobuffer)
    line = read(iobuffer, String)
    println(socket, "line: ", line)

    # Show what character user typed
    println(socket, "character: ", char)
end

#= juliamode.keymap_dict['9'] = debug_mode =#

function funcdump(args...)
    socket = getsocket()
    for (i, arg) in enumerate(args)
        println(socket, "$i $(typeof(arg))")
        println(socket, string("\t", propertynames(arg)))
        println(socket)
    end
end

end


text:
# NOTE: You shouldn't really have to run this
# This makes the Project.toml and Manifest.toml files

# Activate this environment
using Pkg
Pkg.activate(@__DIR__)

# Change package directory
ENV["JULIA_PKG_DEVDIR"] = joinpath(@__DIR__, ".julia")

# Add unregistered packages
# Do this once you've stopped developing
# NOTE: You need access to the VTT GitLab for this
# Pkg.add(PackageSpec(
#         url="git@gitlab.vtt.fi:spine/SpineInterface.jl.git",
#         rev="8c0414860e01222c8c30b7b8e0cfddb0b677f069"
#     )
# )
# Pkg.add(PackageSpec(
#         url="git@gitlab.vtt.fi:spine/model.git",
#         rev="4224cbe0ba3c4f0e9807323c761abe1fa2c8352a"
#     )
# )
# Pkg.add(PackageSpec(
#         url="https://gitlab.kuleuven.be/u0128861/GEPPR.jl.git",
#         rev="d7ef4a9cbb6c024df333984ae761f69ef04ec9a8"
#     )
# )

# For development, uncomment and run the lines below
# NOTE: You'll need to manually change branch after
Pkg.develop(PackageSpec(url="git@gitlab.vtt.fi:spine/SpineInterface.jl.git"))
Pkg.develop(abspath(@__DIR__, ".julia", "SpineOpt"))
Pkg.develop(PackageSpec(url="https://gitlab.kuleuven.be/u0128861/GEPPR.jl.git"))

# Instantiate the other packages (Project.toml file will make sure
# versions are correct)
Pkg.instantiate()


text:
using ModelingToolkit
using Test

@parameters t
@variables a b(t) c d e(t)

b = ParentScope(b)
c = ParentScope(ParentScope(c))
d = GlobalScope(d)

# ensure it works on Term too
LocalScope(e.val)
ParentScope(e.val)
GlobalScope(e.val)

eqs = [0 ~ a
       0 ~ b
       0 ~ c
       0 ~ d]
@named sub4 = NonlinearSystem(eqs, [a, b, c, d], [])
@named sub3 = NonlinearSystem(eqs, [a, b, c, d], [])
@named sub2 = NonlinearSystem([], [], [], systems = [sub3, sub4])
@named sub1 = NonlinearSystem([], [], [], systems = [sub2])
@named sys = NonlinearSystem([], [], [], systems = [sub1])

names = ModelingToolkit.getname.(states(sys))
@test :d in names
@test Symbol("sub1₊c") in names
@test Symbol("sub1₊sub2₊b") in names
@test Symbol("sub1₊sub2₊sub3₊a") in names
@test Symbol("sub1₊sub2₊sub4₊a") in names

@named foo = NonlinearSystem(eqs, [a, b, c, d], [])
@named bar = NonlinearSystem(eqs, [a, b, c, d], [])
@test ModelingToolkit.getname(ModelingToolkit.namespace_expr(ModelingToolkit.namespace_expr(b,
                                                                                foo),
                                                             bar)) == Symbol("bar₊b")

function renamed(nss, sym)
    ModelingToolkit.getname(foldr(ModelingToolkit.renamespace, nss, init = sym))
end

@test renamed([:foo :bar :baz], a) == Symbol("foo₊bar₊baz₊a")
@test renamed([:foo :bar :baz], b) == Symbol("foo₊bar₊b")
@test renamed([:foo :bar :baz], c) == Symbol("foo₊c")
@test renamed([:foo :bar :baz], d) == :d


text:
using BayesianDataFusion
using Base.Test

using DataFrames

A = randn(20,2);
B = randn(30,2);

X = Float64[ sum(A[i,:].*B[j,:]) for i in 1:size(A,1), j in 1:size(B,1) ]

df = DataFrame(A=Int64[], B=Int64[], v=Float64[])
for i=1:size(A,1), j=1:size(B,1)
  push!(df, Any[i, j, X[i,j]])
end

rd = RelationData(df)
rd.entities[1].F = randn(size(A,1), 3)

assignToTest!(rd.relations[1], 10)

tmpdir = mktempdir()
try
  result = macau(rd, burnin=5, psamples=10, num_latent=2, verbose=false, output_beta=true, output="$tmpdir/macau-betasaving", output_type="binary")
  @test isfile("$tmpdir/macau-betasaving-A-01.binary")
  @test isfile("$tmpdir/macau-betasaving-A-01.beta.binary")
  @test isfile("$tmpdir/macau-betasaving-A-02.beta.binary")
  beta_sample1 = read_binary_float32("$tmpdir/macau-betasaving-A-01.beta.binary")
  @test 2 == size(beta_sample1, 2) ## number of latents
  @test 3 == size(beta_sample1, 1) ## number of features
  beta_sample2 = read_binary_float32("$tmpdir/macau-betasaving-A-10.beta.binary")
  beta_last    = convert(Array{Float32}, rd.entities[1].model.beta)
  @test_approx_eq   beta_sample2 beta_last
finally
  rm(tmpdir, recursive = true)
end


text:
"""
    complex_powers!(zpowers, z)

Compute integer powers of `z` from `z^0` through `z^m`, recursively, where `m` is
one less than the length of the input `zpowers` vector.

Note that `z` is assumed to be normalized, with complex amplitude approximately 1.

See also: [`complex_powers`](@ref)
"""
function complex_powers!(zpowers::AbstractVector{Complex{T}}, z::Complex{T}) where {T<:Real}
    Base.require_one_based_indexing(zpowers)
    @fastmath @inbounds begin
        M = length(zpowers)
        if M == 0
            return zpowers
        end
        M -= 1
        θ = one(z)
        zpowers[1] = θ
        if M == 0
            return zpowers
        end
        if M == 1
            zpowers[2] = z
            return zpowers
        end
        while z.re<0 || z.im<0
            θ *= 1im
            z /= 1im
        end
        zpowers[2] = z
        clock = θ
        dc = -2 * sqrt(z).im^2
        t = 2 * dc
        dz = dc * (1 + 2 * z) + 1im * sqrt(-dc * (2 + dc))
        for m in 2:M
            zpowers[m+1] = zpowers[m] + dz
            zpowers[m] *= clock
            clock *= θ
            dz += t * zpowers[m+1]
        end
        zpowers[M+1] *= clock
    end
    zpowers
end

"""
    complex_powers(z, m)

Compute integer powers of `z` from `z^0` through `z^m`, recursively.

Note that `z` is assumed to be normalized, with complex amplitude approximately 1.

This algorithm is mostly due to Stoer and Bulirsch in "Introduction to
Numerical Analysis" (page 24) — with a little help from de Moivre's formula,
which is essentially exp(iθ)ⁿ = exp(inθ), as well as my own alterations to deal
with different behaviors in different quadrants.

There isn't usually a huge advantage to using this specialized function.  If
you just need a particular power, it will generally be far more efficient and
just as accurate to compute either exp(iθ)ⁿ or exp(inθ) explicitly.  However,
if you need all powers from 0 to m, this function is about 10 or 5 times faster
than those options, respectively, for large m.  Like those options, this function
is numerically stable, in the sense that its errors are usually smaller than `m`
times the error from machine-precision errors in the input argument — or at worst
about 50% larger, which occurs as the phase approaches multiples of π/2.

See also: [`complex_powers!`](@ref)
"""
function complex_powers(z::Complex{T}, m::Int) where {T<:Real}
    if abs(z) ≉ one(z)
        throw(DomainError("z = $z",
            "This function is only valid for `z` with complex amplitude approximately 1; "
            * "abs(z) = $(abs(z))"
        ))
    end
    zpowers = zeros(typeof(z), m+1)
    complex_powers!(zpowers, z)
end