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

using Markdown
using InteractiveUtils

# ╔═╡ 3bf40bbf-539c-45f5-956c-b9650b08c1d5
using Unitful, Statistics

# ╔═╡ 22f27dce-3ed6-41df-92e1-712eb08d5c88
using DataFrames, DataFramesMeta

# ╔═╡ e303bf14-1aa4-4d4f-884a-cd4f43f82186
begin
	using CommonMark,
		LaTeXStrings,
		PlutoUI,
		HypertextLiteral,
		Lazy,
		MarkdownLiteral,
		Hyperscript

	import MarkdownLiteral: @markdown

	md"Load Julia \"infrastructure\" pkgs"
end

# ╔═╡ 4d965bb6-567e-4f05-8cb0-ef49349c5f12
cm"""

Given
-----

A piston–cylinder device contains 0.85 kg of refrigerant- 134a at –10°C. The piston that is free to move has a mass of 12 kg and a diameter of 25 cm. The local atmospheric pressure is 88 kPa. Now, heat is transferred to refrigerant-134a until the temperature is 15°C.
"""

# ╔═╡ 9804b746-3ad2-4b9e-af82-d3adec5eeda3
cm"""

{style="width:35%; margin-inline:auto;"}
![head](https://imgur.com/cXHhjIw.png)

"""

# ╔═╡ 34570e7a-60ad-4851-ac9a-356e866c0146
cm"""

Find
----

1. the final pressure,
1. the change in the volume of the cylinder, and
1. the change in the enthalpy of the refrigerant-134a.
"""

# ╔═╡ a4b02708-524c-40af-8d94-ab3272c3b81d
cm"""

Get organized
-------------

These things you must notice right away:

{style="list-style-type: lower-greek;"}
1. Volume change ``\Rightarrow`` _boundary work_{.killerorange}.  No other forms of work
1. "... _piston that is free to move..._" ``\Rightarrow`` _frictionless piston_{.killerorange}
1. No variable forces on piston + frictionless piston ``\Rightarrow`` _isobaric process_{.killerorange}
1. Properties/characteristics away from the norm
	1. Fluid is R-134a, not water.
	1. Atmospheric pressure = 88 kPa

This is a classic before/after scenario. The information we're tasked to find describes what happened to the system in its move from one state to another. All we need to do is resolve both states, then calculate the property changes requested.

"""

# ╔═╡ de2ef376-c065-4244-b455-d91ff78ffcd8
cm"""
We're not overtly given enough information to resolve either state, but we are given what we need to calculate ``P_1``. Moreover, we know that our states are connected by an isobaric process.

Moreover, we are given enough information Let's help ourselves by clearly expressing our states with respect to the properties we know and the properties we need:

```math
\left.
\begin{aligned}
T_1 &= \pu{-10 °C}  \\
P_1 &= ?
\end{aligned}
\right\} \ce{①->[\quad isobaric \quad]②}
\left\{
\begin{aligned}
T_2 &= \pu{15 °C}  \\
P_2 &= P_1
\end{aligned}
\right.
```

Checking inventory, we need
1. ``P_1`` : calculate it
1. ``\Delta V = m (v_2 - v_1)`` : ``v`` in tables
1. ``\Delta H = m (h_2 - h_1)`` : ``h`` in tables

"""

# ╔═╡ e01ff0d4-103b-4ca6-aa3f-7bce7a690083
begin
  # given
	m_R134a  = 0.85u"kg"
	m_piston = 12u"kg"
	P_atm    = 88u"kPa"
	D_piston = 25u"cm"

  # calc
	A_piston = π/4 * D_piston^2
end;

# ╔═╡ 794342fc-96f6-49fa-bfa3-3c3229c49e74
P1 = P_atm + m_piston * 1u"ge" / A_piston |> u"kPa"

# ╔═╡ cd80963f-bdab-406f-bb5c-7fc3b4e849df
cm"""

Now that we have the pressure, we can go to the property tables to resolve our states.

![head](https://imgur.com/VoARkJf.png)
![body](https://imgur.com/KV702G7.png)

As usual, we head to the saturation tables to find our state. There's no break point at our temperature, but we don't need it. Our pressure is lower than the saturation pressures of the temperatures bracketing our temperature; we can safely say that our pressure would be lower than the saturation pressure belonging to our temperature. That tells us to head to the Superheated Vapor tables.

"""

# ╔═╡ 11e873c6-f27c-4251-a9df-a09e59a9267c
cm"""
![body](https://imgur.com/gfy1P9I.png)

In the superheated vapor tables we find the neighborhood in which out states exist, but only one of four properties lands on a breakpoint. We need to interpolate in two directions, affectionately known as _bilinear interpolation._

In bilinear interpolation, we interpolate in one direction, then interpolate that result in the other direction. There are techniques for interpolating both directions simultaneously, but they're rarely needed in undergraduate thermodynamics.
"""

# ╔═╡ b1d33668-3ed0-4ee0-b14a-78d501929cbc
cm"""

Attack
------

We will interpolate pressure first, mainly because of personal preference. We are essentially creating a subtable at our pressure. Then we interpolate the new subtable for temperature. Liberal subscripting makes this manageable. Julia helps with generous latitude in naming variables.

The only thing that goes wrong is losing track of which terms go into what positions in the equations. Unfortunately, it is painfully easy to lose track of the terms.

I **strongly**{.killerorange} recommend you lay everything out algebraically first, don't combine any steps, and use a spreadsheet, dataframes (lean and mean database), or other computational tools to do the actual arithmetic.

I used subscripting for readability. My tool of choice would be a dataframe. At the end of this problem, I'll show a dataframe inplementation.
"""

# ╔═╡ f81bb289-fd95-442f-8b86-dd59a9df685c
# from R-134a superheated vapor table

begin
	v_60kPa_⁻10°C  = 0.34992u"m^3/kg"
	v_60kPa_10°C   = 0.37861u"m^3/kg"
	v_60kPa_20°C   = 0.39279u"m^3/kg"

	v_100kPa_⁻10°C = 0.20686u"m^3/kg"
	v_100kPa_10°C  = 0.22473u"m^3/kg"
	v_100kPa_20°C  = 0.23349u"m^3/kg"

	h_60kPa_⁻10°C  = 245.96u"kJ/kg"
	h_60kPa_10°C   = 262.41u"kJ/kg"
	h_60kPa_20°C   = 270.89u"kJ/kg"

	h_100kPa_⁻10°C = 244.70u"kJ/kg"
	h_100kPa_10°C  = 261.43u"kJ/kg"
	h_100kPa_20°C  = 270.02u"kJ/kg"
end;

# ╔═╡ 751e9afd-6938-45cd-bd99-2a8aae5f1075

# ╔═╡ a510380b-8791-4339-857f-136e89eaa931
# interpolate across P for all states
# P1 = P2, but labelling them separately helps readability

begin
	𝑓p = (P1 - 60u"kPa")/(100u"kPa" - 60u"kPa")

	v_P1_⁻10°C = v_60kPa_⁻10°C + 𝑓p*(v_100kPa_⁻10°C - v_60kPa_⁻10°C)
	h_P1_⁻10°C = h_60kPa_⁻10°C + 𝑓p*(h_100kPa_⁻10°C - h_60kPa_⁻10°C)

	v_P2_10°C  = v_60kPa_10°C  + 𝑓p*(v_100kPa_10°C  - v_60kPa_10°C)
	v_P2_20°C  = v_60kPa_20°C  + 𝑓p*(v_100kPa_20°C  - v_60kPa_20°C)

	h_P2_10°C  = h_60kPa_10°C  + 𝑓p*(h_100kPa_10°C  - h_60kPa_10°C)
	h_P2_20°C  = h_60kPa_20°C  + 𝑓p*(h_100kPa_20°C  - h_60kPa_20°C)
end;

# ╔═╡ c11e59c2-0c7d-4bde-9d02-c2ae9a2cc80e
# interpolate across T for State 2

begin
	v_P2_15°C  = [v_P2_10°C, v_P2_20°C] |> mean
	h_P2_15°C  = [h_P2_10°C, h_P2_20°C] |> mean
end;

# ╔═╡ b8fd73e6-1408-41d4-9ca2-0406d3a4dc1d
# final property calculations

begin
	v1 = v_P1_⁻10°C
	v2 = v_P2_15°C

	h1 = h_P1_⁻10°C
	h2 = h_P2_15°C
end;

# ╔═╡ 6a60a492-812b-4273-94d9-e9cc744b148c
ΔH = m_R134a * (h2 - h1)

# ╔═╡ dbcdc200-3044-4376-bf0a-4adc103f4681
ΔV = m_R134a * (v2 - v1)

# ╔═╡ 32c816f3-ed9d-485a-8d78-b6f8c2735bf9
cm"""
----

{.message .success}
> _finis_{.success}
"""

# ╔═╡ 04496cd3-ed1a-4a77-bd6b-b3e80780fe16
begin

	P_  = round(P1 |> typeof, P1, digits = 5)
	ΔV_ = round(ΔV |> typeof, ΔV, digits = 5)
	ΔH_ = round(ΔH |> typeof, ΔH, digits = 5)

	L"""
	\begin{align}
		P &= \pu{%$P_} \\
		\Delta V &= \pu{%$ΔV_} \\
		\Delta H &= \pu{%$ΔH_}
	\end{align}
	"""
end

# ╔═╡ c53969f6-5ec7-4448-bdf7-28b7819b6d33
cm"----
**Show/tell:**{.killerorange} Dataframes
"

# ╔═╡ 39d64047-5eba-42a1-a78c-dfc1fb4200a4
R134a = DataFrame(
	P = [     60,      60,      60,     100,     100,     100]u"kPa"    ,
	T = [ 	 -10,      10,      20,     -10,      10,      20]u"°C"     ,
	v = [0.34992, 0.37861, 0.39279, 0.20686, 0.22473, 0.23349]u"m^3/kg" ,
	h = [ 245.96,  262.41,  270.89,  244.70,  261.43,  270.02]u"kJ/kg"
)

# ╔═╡ 92ac28e6-43a6-4e4a-9e52-a065bbb31cde
@subset(R134a, :P .== 60u"kPa")

# ╔═╡ 552fb87a-b1e0-4b74-8be6-25dcb40a1294
@subset(R134a, :P .== 60u"kPa").h

# ╔═╡ 7db6ecdb-c6d3-49ae-85e6-569271c23285
@subset(R134a, :P .== 60u"kPa", :T .== 20u"°C")

# ╔═╡ aadbbc39-6702-4ea2-9167-c35efade6740
@subset(R134a, :P .== 60u"kPa", :T .== 20u"°C").v

# ╔═╡ d9120f9a-d084-471b-8bf7-250674f63b0a
cm"## Appendix"

# ╔═╡ 5dc45506-e5c8-4992-902b-5f8d93f31dcd
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
					""",
		"ms8" 	=> 	"""
					Moran, MJ Ed (2014)
					Fundamentals of engineering thermodynamics 8e,
					978-1118412930
					""",
		"ms9" 	=> 	"""
					Moran, MJ, _et al_ (2018)
					Fundamentals of engineering thermodynamics 9e,
					978-1119391388
					""",
		"cb7" 	=> 	"""
					Çengel, YA & MA Boles (2011)
					Thermodynamics: An engineering approach 7e,
					978-0073529325
					""",
		"cb8" 	=> 	"""
					Çengel, YA & MA Boles (2015)
					Thermodynamics: An engineering approach 8e,
					978-0073398174
					""",
		"cb9" 	=> 	"""
					Çengel, YA, Boles, MA, & M Kanoglu (2019)
					Thermodynamics: An engineering approach 9e,
					978-1259822674
					""",
	)

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

	the_bibliography =
		@as xx bib[bookID] begin
			split(xx, '\n', keepempty=false)
			strip.(xx)
			join(xx, " _", "_ ")
		end

	# return the_bibliography

	return @markdown("""

			   # $title

			   {.source}
			   $the_bibliography

			""")

end

# ╔═╡ 806d1314-9909-11ec-05f5-c194fa8b2722
title"Problem 3.30"cb9

# ╔═╡ c2a438b9-8be7-42f7-9871-293f8fd19790
begin
	mcs_css_url = "https://gitlab.com/marcucius/mcsBibliotheca/-/raw/main/styles/mcs-pluto.css"

	mcs_css = open(download(mcs_css_url)) do io
	    read(io, String)
	end

	@htl """

Load stylesheets

	"""

end

# ╔═╡ 592a917a-7b3b-4ae8-85dd-14726efdc373
# TableOfContents()

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CommonMark = "a80b9123-70ca-4bc0-993e-6e3bcb318db6"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
DataFramesMeta = "1313f7d8-7da2-5740-9ea0-a2ca25f37964"
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"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"

[compat]
CommonMark = "~0.8.6"
DataFrames = "~1.3.2"
DataFramesMeta = "~0.10.0"
Hyperscript = "~0.0.4"
HypertextLiteral = "~0.9.3"
LaTeXStrings = "~1.3.0"
Lazy = "~0.15.1"
MarkdownLiteral = "~0.1.1"
PlutoUI = "~0.7.35"
Unitful = "~1.11.0"
"""

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

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

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

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

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

[[deps.Chain]]
git-tree-sha1 = "339237319ef4712e6e5df7758d0bccddf5c237d9"
uuid = "8be319e6-bccf-4806-a6f7-6fae938471bc"
version = "0.4.10"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

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

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

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

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

[[deps.Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

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

[[deps.DataFramesMeta]]
deps = ["Chain", "DataFrames", "MacroTools", "OrderedCollections", "Reexport"]
git-tree-sha1 = "ab4768d2cc6ab000cd0cec78e8e1ea6b03c7c3e2"
uuid = "1313f7d8-7da2-5740-9ea0-a2ca25f37964"
version = "0.10.0"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

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

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

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

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

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

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

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

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

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

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

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

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

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

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

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "13468f237353112a01b2d6b32f3d0f80219944aa"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.2"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# ╔═╡ Cell order:
# ╟─806d1314-9909-11ec-05f5-c194fa8b2722
# ╟─4d965bb6-567e-4f05-8cb0-ef49349c5f12
# ╠═9804b746-3ad2-4b9e-af82-d3adec5eeda3
# ╟─34570e7a-60ad-4851-ac9a-356e866c0146
# ╟─a4b02708-524c-40af-8d94-ab3272c3b81d
# ╟─de2ef376-c065-4244-b455-d91ff78ffcd8
# ╠═3bf40bbf-539c-45f5-956c-b9650b08c1d5
# ╠═e01ff0d4-103b-4ca6-aa3f-7bce7a690083
# ╠═794342fc-96f6-49fa-bfa3-3c3229c49e74
# ╟─cd80963f-bdab-406f-bb5c-7fc3b4e849df
# ╟─11e873c6-f27c-4251-a9df-a09e59a9267c
# ╟─b1d33668-3ed0-4ee0-b14a-78d501929cbc
# ╠═f81bb289-fd95-442f-8b86-dd59a9df685c
# ╠═751e9afd-6938-45cd-bd99-2a8aae5f1075
# ╠═a510380b-8791-4339-857f-136e89eaa931
# ╠═c11e59c2-0c7d-4bde-9d02-c2ae9a2cc80e
# ╠═b8fd73e6-1408-41d4-9ca2-0406d3a4dc1d
# ╠═6a60a492-812b-4273-94d9-e9cc744b148c
# ╠═dbcdc200-3044-4376-bf0a-4adc103f4681
# ╟─32c816f3-ed9d-485a-8d78-b6f8c2735bf9
# ╟─04496cd3-ed1a-4a77-bd6b-b3e80780fe16
# ╟─c53969f6-5ec7-4448-bdf7-28b7819b6d33
# ╠═22f27dce-3ed6-41df-92e1-712eb08d5c88
# ╠═39d64047-5eba-42a1-a78c-dfc1fb4200a4
# ╠═92ac28e6-43a6-4e4a-9e52-a065bbb31cde
# ╠═552fb87a-b1e0-4b74-8be6-25dcb40a1294
# ╠═7db6ecdb-c6d3-49ae-85e6-569271c23285
# ╠═aadbbc39-6702-4ea2-9167-c35efade6740
# ╟─d9120f9a-d084-471b-8bf7-250674f63b0a
# ╟─e303bf14-1aa4-4d4f-884a-cd4f43f82186
# ╟─5dc45506-e5c8-4992-902b-5f8d93f31dcd
# ╟─c2a438b9-8be7-42f7-9871-293f8fd19790
# ╠═592a917a-7b3b-4ae8-85dd-14726efdc373
# ╟─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:
#=
ExtenderNSmoother.jl
Jacob Child
May 20th, 2022
=#

ext: jl


text:
#### PRREQUISITES ####
 include("objectives.jl")

#### for module multiobjectiveGraphs ####
export  multiobj # , getObj, setObj!, summarize
# export descr, defaultValue, mimic, combine # already in objectives.jl
# export lexicoBetter, dominates, dominatesStrictly, objLength

import Base.size, Base.length # Julia 1+ , allow redefinition of theses functions

##############################################################
#################### generic type for multi-objective values     #########
##############################################################

#### note : all objectives do not need to share the same type of values ##################

mutable struct multiobj # <: genericMultiobj # Tobjval is union of all the Tobjvals of the objectives
	nb::UInt8;  ######### number of objectives is (severely :-) limited to 255
	objectives::Vector@{ genericWeightCategory}
	############## constructors
	function multiobj()
			## create an monoobj with initialised array of  1 summable objective
		ar= Vector{ genericWeightCategory }(undef, 1)    # Vector{ weightMinSum{Tobjval} }( 1)  for Julia 0.63
		ar[1]=weightMinSum{Float64}()
		return new(1,ar)
	end													##############

	function multiobj( mo::multiobj )
	 	return new( mo.nb, mo.objectives ) ## clone with same types and same values
	end 												##############
	function multiobj(n::Number; typeofobjectives=weightMinSum{Float64})
			## create an multiobj with initialised array of summable objectives
		nnn::Int8=convert(Int8,n)
		ar= Vector{ genericWeightCategory }(undef,  nnn )  # suppress undef first parameter for Julia 0.63
		for i::Int64 in 1:nnn 	ar[i]=typeofobjectives() 	end
		return new(nnn,ar)
	end													##############
	function multiobj(p::Number,q::Number,r::Number ; typeofvalues=Float64)
			## create an multiobj with initialised array of
			## p summable objectives, q multiplicative and r bottleneck
		ar=Vector{genericWeightCategory}(undef, p+q+r) # suppress undef first parameter for Julia 0.63
		for i::Int64 in 1:p 	ar[i]=weightMinSum{typeofvalues}() 	end
		for i::Int64 in (p+1):(p+q) 	ar[i]=weightMinProduct{typeofvalues}() 	end
		for i::Int64 in (p+q+1):(p+q+r) 	ar[i]=weightMinMax{typeofvalues}() 	end
		return new(p+q+r,ar)
	end													##############
	function multiobj(mobjTypes::Vector{genericWeightCategory } )
			## create an multiobj with initialised types of weights (and matching default values)
		nn = length(mobjTypes)
		ar = Vector{genericWeightCategory}(undef,  nn ) # suppress undef first parameter for Julia 0.63
		   # with the same types in the same indices
		for i::Int64 in 1:nn 	ar[i]=mimic(mobjTypes[i], defaultValue( mobjTypes[i] ) ) 	end
		return new{Tobjval}( nn, ar )
	end
	function multiobj(mobjTypes::Array{genericWeightCategory } )
			## create an multiobj with initialised types of weights (and matching default values)
		nn = length(mobjTypes)
		ar = Vector{genericWeightCategory}(undef,  nn ) # suppress undef first parameter for Julia 0.63
		   # with the same types in the same indices
		for i::Int64 in 1:nn 	ar[i]=mimic(mobjTypes[i], defaultValue( mobjTypes[i] ) ) 	end
		return new{Tobjval}( nn, ar )
	end													##############
	function multiobj(mobjTypes::Vector{ genericWeightCategory },
								mobjVals::Vector{Tvals} )  where  Tvals <: Number
			## create an multiobj with array of objectives types and initialised array of values of weights
		nn = length(mobjVals)
		ar = Vector{genericWeightCategory}(undef, nn )  # suppress undef first parameter for Julia 0.63
		   # with the same types in the same indices
		for i::Int64 in 1:nn 	ar[i]=mimic(mobjTypes[i], mobjVals[i] ) 	end
		return new{Tobjval}( nn, ar )
	end
	function multiobj(a::multiobj,
								mobjVals::Vector{Tvals} )  where Tvals <: Number
			## create an multiobj with model of objectives types and initialised array of values of weights
		nn = length(mobjVals)
		ar = Vector{genericWeightCategory}(undef, nn )  # suppress undef first parameter for Julia 0.63
		   # with the same types in the same indices
		for i::Int64 in 1:nn 	ar[i]=mimic(a.objectives[i], mobjVals[i] ) 	end
		return new( nn, ar )
	end													##############
	function multiobj(a::multiobj,
								mobjVals::Array{Tvals} )  where Tvals <: Number
			## create an multiobj with model of objectives types and initialised array of values of weights
		nn = length(mobjVals)
		ar = Vector{genericWeightCategory}(undef, nn )  # suppress undef first parameter for Julia 0.63
		   # with the same types in the same indices
		for i::Int64 in 1:nn 	ar[i]=mimic(a.objectives[i], mobjVals[i] ) 	end
		return new( nn, ar )
	end													##############
	function multiobj(a::multiobj , b::multiobj )
			## create an multiobj with combinations of weights in a && b
		nn::Int8 = objLength(a)
		ar= Vector{ genericWeightCategory }(undef,  nn )  # a && b MUST be the same size  # suppress undef first parameter for Julia 0.63
		   # with the same types in the same indices
		for i::Int64 in 1:nn 	ar[i]=combine(a.objectives[i],b.objectives[i]) 	end
		return new(nn,ar)
	end
end

######### accessor ##############
 function getObj( mo::multiobj, numobj::Integer )
   return mo.objectives[numobj]
 end
 function setObj!( mo::multiobj, numobj::Integer, objvalue)
   old = mo.objectives[numobj] ; mo.objectives[numobj] = objvalue ; return old
 end
######### displaying a multi-objective ############
 function summarize( mo::multiobj ; log=stdout)
 	print(log,"[ ")
 	 for obj in mo.objectives
 	 	print(log,obj.value," ")
 	 end
 	 print(log,"]");  #### be carefull using it : no line break here
 end

 ######### description of a multi-objective ############
 function descr( mo::multiobj ; withName=true )
  txt::String="[ "
  if withName
 	 txt*=descr(mo.objectives[1])
 	 for i in 2:length(mo.objectives) txt*=","*descr(mo.objectives[i])  end
 	 txt*="]"
  else
 	 txt*=string(mo.objectives[1].value)
 	 for i in 2:length(mo.objectives)  txt*=","*string(mo.objectives[i].value)  end
 	 txt*="]"
  end
  return txt
 end

############## initialisation values of given types
	function defaultValue( mo::multiobj )
		return multiobj( mo.objectives )
	end
############## initialisation values of given types
	function defaultValue( objTypes::Vector{ genericWeightCategory } )
		return multiobj( objTypes )
	end
############## same types with default values ( same than preceding )
	function mimic( objTypes::Vector{ genericWeightCategory} )
	 	return multiobj( objTypes )
	end
############## same types with other values
	function mimic( objTypes::Vector{ genericWeightCategory },
								objVals::Vector{Tvals} ) where   Tvals <: Number
	 	return multiobj( objTypes, objVals )
	end

	function mimic( objTypes::Vector{ genericWeightCategory },
								objVals::Array{Tvals} ) where  Tvals <: Number
	 	return multiobj( objTypes, objVals )
	end
############## same types with other values
	function mimic( objTypes::multiobj,
									objVals::Vector{Tvals} ) where  Tvals <: Number
		 return multiobj( objTypes, objVals )
	end
############## clone
	function mimic( mo::multiobj )
	 	return multiobj( mo )
	end

############# combination ##########
function combine(a::multiobj, b::multiobj )::multiobj
 return multiobj(a , b )
end

############# lexicographic order ##########
function lexicoBetter(a::multiobj, b::multiobj )::Bool
	num::Int8 = 1 ; maxNum::Int8 = objLength(a)
	while num < maxNum && compare(a.objectives[num], b.objectives[num])==(true,true) ## sequential evaluation not needed here
		num += 1
	end
	return compare(a.objectives[num], b.objectives[num])[1]
end
############# lexicographic order cycling, beginning with k-th  ##########
function lexicoBetter(a::multiobj, b::multiobj, k::Int8 )::Bool
	cntr::Int8 = 1 ; maxNum::Int8 = objLength(a) ; num::Int8 =k
	while cntr < maxNum && compare(a.objectives[num], b.objectives[num])==(true,true) ## sequential evaluation not needed here
		cntr += 1 ; num = (num % maxNum) + 1 ;
	end
	return compare(a.objectives[num], b.objectives[num])[1]
end
############# dominance ##########
function compare(a::multiobj, b::multiobj )::Tuple{Bool,Bool}
              ## returns (a<=b, b<=a) ; hence returns  (true,true) if  a == b
 			i::Int8=1; aisbetter::Bool=true  # better and worse must be seen here as large ( _ or equal)
 			while i <= a.nb &&  aisbetter                                # a.nb is assumed == b.nb
 					aisbetter = aisbetter && compare(a.objectives[i],  b.objectives[i])[1]
 				    i += 1
 			end
 			i=1; aisworse::Bool=true
 			while i <= a.nb &&  aisworse                               # a.nb is assumed == b.nb
 					aisworse = aisworse && compare(a.objectives[i], b.objectives[i])[2]
 				    i += 1
 			end
 			return  aisbetter, aisworse
end

function dominates(a::multiobj, b::multiobj )::Bool
  cp = compare(a,b)
  return cp[1]
end

function dominatesStrictly(a::multiobj, b::multiobj )::Bool
  cp = compare(a,b)
  return cp[1] && ! cp[2]
end

######### specific dominance (one objective only -> Dijsktra) ###########
function dominates(a::multiobj, b::multiobj, objnum::Int8 )::Bool
  return compare(a.objectives[objnum],  b.objectives[objnum])[1]
end

############## some well known functions redefined     ##########
function size(mo::multiobj)::Int8  where Tobjval<:Number  return mo.nb end
function length(mo::multiobj)::Int8  where Tobjval<:Number  return mo.nb end
function objLength(mo::multiobj)::Int8  where Tobjval<:Number return mo.nb end

println("multiObjectives loaded")

max_stars_repo_name: DRobbes/multiobjectiveGraphs

max_issues_repo_name: DRobbes/multiobjectiveGraphs

max_forks_repo_name: DRobbes/multiobjectiveGraphs


text:
using DifferentialEquations
using IterableTables, DataFrames
using NamedTuples
using Gadfly

function F(du,u,p,t)
    S_H, E_H, I_H, R_H, S_V, E_V, I_V = u

    # host dynamics
    host_infection = (p.β*S_H*I_V)/p.N_H
    host_mortality = p.μ_H .* u[1:4] # include S_H, so easier to remove mortality
    host_births = sum(host_mortality)
    host_progression = p.σ_H*E_H
    recovery = p.λ*I_H

    du[1] = -host_infection + host_births
    du[2] = host_infection - host_progression
    du[3] = host_progression - recovery
    du[4] = recovery
    du[1:4] -= host_mortality

    # vector dynamics
    vec_infection = (p.β*S_V*I_H)/p.N_H
    vec_mortality = p.μ_V .* u[5:7] # include S_V, so easier to remove mortality
    vec_births = sum(vec_mortality)
    vec_progression = p.σ_V*E_V

    du[5] = -vec_infection + vec_births
    du[6] = vec_infection - vec_progression
    du[7] = vec_progression
    du[5:7] -= vec_mortality

end

# nb: in >= Julia v0.7, can eliminate this import
#  and the @NT syntax
u0 = [
    S_H=100.0,   E_H=0.0, I_H=1.0, R_H=0.0,
    S_V=10000.0, E_V=0.0, I_V=0.0
]
p = @NT(
  μ_H=1/365, μ_V=1/30, σ_H=1/3, σ_V=1/7, λ=1/14,
  β=0.05, N_H = sum(u0[1:4])
)
tspan = (0.0, 365.0)
prob = ODEProblem(F, u0, tspan, p)
sol = @time solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8,saveat=linspace(0,365,365*10+1))

df = DataFrame(sol)
rename!(df, Dict(:timestamp => :t,
  :value1 => :S_H, :value2 => :E_H, :value3 => :I_H, :value4 => :R_H,
  :value5 => :S_V, :value6 => :E_V, :value7 => :I_V
))
mlt = melt(df,:t) # convert results into long format for plotting
mlt[:host] = contains.(string.(mlt[:variable]),"H"); # tag which entries are host vs vector
df

# Plot
fig1a = plot(mlt[mlt[:host] .== true,:], x=:t, y=:value, color=:variable, Geom.line)
fig1b = plot(mlt[mlt[:host] .!= true,:], x=:t, y=:value, color=:variable, Geom.line)
vstack(fig1a,fig1b)

max_stars_repo_name: epimodels/epicookbook

max_issues_repo_name: epimodels/epicookbook

max_forks_repo_name: epimodels/epicookbook


text:
module Rocketland
#using Mosek
using MathOptInterface
using ECOS
using MosekTools
using Rotations
using StaticArrays
using LinearAlgebra
using ForwardDiff
using ..RocketlandDefns
import ..Dynamics
import ..FirstRound

#state indexing
const state_dim = 14
const control_dim = 3
const r_idx_it = 2:4
const v_idx_it = 5:7
const qbi_idx_it = 8:11
const omb_idx_it = 12:14
const acc_width = state_dim+control_dim*2+3
const acc_height = state_dim
const mass_idx = 1

#=
how to use
julia> Dynamics.make_dynamics_module(RocketlandDefns.ProbInfo(SampleProblems.base_prob_aero_scaled))
julia> cache=Dynamics.IntegratorCache(SampleProblems.base_prob_aero_scaled, RocketlandDefns.ProbInfo(SampleProblems.base_prob_aero_scaled), Linearizer);
julia> pi = Rocketland.create_initial(SampleProblems.base_prob_aero_scaled, cache);

=#

function create_initial(problem::DescentProblem, linear_cache::IntegratorCache)
    K = problem.K
    initial_points,linpoints = FirstRound.linear_initial(problem, linear_cache)
    model = build_model(problem, K, linpoints, initial_points, problem.tf_guess)
    return ProblemIteration(problem, linear_cache, problem.tf_guess, initial_points, linpoints, model, 0, 100.0, Inf)
end

const MOI=MathOptInterface
const SAF=MOI.ScalarAffineFunction
const VAF=MOI.VectorAffineFunction
const SA=MOI.ScalarAffineTerm
const VA=MOI.VectorAffineTerm
const VoV=MOI.VectorOfVariables
const SOC=MOI.SecondOrderCone

# SCvx defns
# Lk = -1 x[1, K+1] + 1e4 * norm(nuv)
# Jk = -1 x[1, K+1] + 1e4 * norm(x[:, i+1] - act[i+1])

function build_model(prob, K, iterDynam, iterAbout, sigHat)
<<<<<<< HEAD
	#model =
=======
>>>>>>> 6d18d726748568198fe41d32ebcf2420ce96670c
	model = MOI.instantiate(() -> Mosek.Optimizer(#=MSK_IPAR_LOG=0,=#MSK_IPAR_INFEAS_REPORT_AUTO=1,MSK_IPAR_BI_IGNORE_MAX_ITER=1,
				MSK_IPAR_INTPNT_MAX_ITERATIONS=10000); with_bridge_type=Float64)
	dcs = MOI.ConstraintIndex[]
	state_nuc = MOI.ConstraintIndex[]

	tggs = tand(prob.gammaGs)
	sqcm = sqrt((1-cosd(prob.thetaMax))/2)
	delMax = cosd(prob.deltaMax)
	cosma = cosd(45.0)

	all_state_dim = state_dim*(K+1)
	all_control_dim = control_dim*(K+1)
	#main state variables
	xv = reshape(MOI.add_variables(model, all_state_dim), state_dim, K+1)
	uv = reshape(MOI.add_variables(model, all_control_dim), control_dim, K+1)
	dxv = reshape(MOI.add_variables(model, all_state_dim), state_dim, K+1)
	duv = reshape(MOI.add_variables(model, all_control_dim), control_dim, K+1)
	dsig = MOI.add_variable(model)
	nuv = reshape(MOI.add_variables(model, all_state_dim), state_dim, K+1)

	# trust regions
	Jvnu = MOI.add_variable(model)
	Jtr = MOI.add_variable(model)
	Jsig = MOI.add_variable(model)

	#objective
	MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
		SAF([SA(-1.0, xv[1, K+1]), SA(prob.wNu, Jvnu), SA(0.5, Jtr), SA(1.0, Jsig)],0.0))
	MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

	# couple the xs and us to the dxs and the dus
	# about + dx = x; about + du = u; etc
    # vcat(iterAbout[:].state...) + reshape(xv, all_state_dim) - reshape(xv, all_state_dim) = 0
    # vcat(iterAbout[:].control...) + reshape(duv, all_control_dim) - reshape(uv, all_control_dim) = 0
    state_base = MOI.add_constraint(model, VAF(VA.([1:all_state_dim; 1:all_state_dim],
    	[SA.(1.0, reshape(dxv, all_state_dim)); SA.(-1.0, reshape(xv, all_state_dim))]),
    	Array(vcat(getfield.(iterAbout[:], :state)...))), MOI.Zeros(all_state_dim))
    control_base = MOI.add_constraint(model, VAF(VA.([1:all_control_dim; 1:all_control_dim],
    	[SA.(1.0, reshape(duv, all_control_dim)); SA.(-1.0, reshape(uv, all_control_dim))]),
    	Array(vcat(getfield.(iterAbout[:], :control)...))), MOI.Zeros(all_control_dim))

    # build the trust regions
    MOI.add_constraint(model, VoV([Jvnu; reshape(nuv, all_state_dim)]), SOC(all_state_dim + 1))
    MOI.add_constraint(model, VoV([Jtr; reshape(dxv, all_state_dim); reshape(duv, all_control_dim)]), SOC(1+all_state_dim+all_control_dim))
    MOI.add_constraint(model, VoV([Jsig; dsig]), SOC(2))

	#initial and final state constraints
	# r_1 = r_i, v_1 = v_i, omb_1 = omb_i; same for final
	# we assume that the base solution satisfies the required state constraints. Therefore, it is sufficient
	# to lock the differences to 0
	# done as one big vector equality constraint; 1 * vars - vals = 0 -> vars = vals
	vars = map(t -> VA(t...), enumerate(
		SA.(1.0, [xv[mass_idx, 1]; xv[r_idx_it,1]; xv[v_idx_it,1]; xv[omb_idx_it,1];
				  xv[r_idx_it,K+1]; xv[v_idx_it,K+1]; xv[qbi_idx_it,K+1]; xv[omb_idx_it,K+1];
				  uv[2:3, K+1]])))
	vals = -[prob.mwet; prob.rIi; prob.vIi; prob.wBi;
			            prob.rIf; prob.vIf; prob.qBIf; prob.wBf; 0.0; 0.0]
	MOI.add_constraint(model, VAF(vars, vals), MOI.Zeros(length(vars)))

	# linearized dynamics
	# where inp_n = [state_n; control_n; control_n+1; sigma]
	# at each step, ostate_{n+1} + dstate_{n+1} = linearized*(lininp_n + dinp_n) + nu
	# at each step, ostate_{n+1} + dstate_{n+1} = endpoint + linearized*dinp_n + nu
	# at each step, dstate_{n+1} = linearized*dinp_n + nu + (endpoint - ostate_{n+1})
	# so linearized*[dstate_n; dcontrol_n; dcontrol_n+1; dsigma] + nu + (endpoint - ostate_{n+1}) - dstate_{n+1} = 0
	dcstrs = MOI.ConstraintIndex[]
	for n=1:K
		matmul = vcat(map((col,var) -> vcat(map(x->VA(x...), enumerate(map(x -> SA(x,var), col)))...),
					      eachcol(iterDynam[n].derivative), [dxv[:,n]; duv[:,n]; duv[:,n+1]; dsig])...)
		nu_eqn = map(invar -> VA(invar[1], SA(1.0, invar[2])), enumerate(nuv[:,n+1]))
		eq_eqn = map(invar -> VA(invar[1], SA(-1.0, invar[2])), enumerate(dxv[:,n+1]))
		lin_err = Array(iterDynam[n].endpoint - iterAbout[n+1].state)
		func = VAF(vcat(matmul, nu_eqn, eq_eqn), lin_err)
		cstr = MOI.add_constraint(model, func, MOI.Zeros(state_dim))
		push!(dcstrs, cstr)
	end

	# state constraints wooo; equation nums from Szmuk and Acikmese 2018
	# eq. 6; mdry <= mk
	MOI.add_constraint(model, VAF(VA.(1:K, SA.(1.0, xv[1,2:K+1])), fill(-prob.mdry, K)), MOI.Nonnegatives(K))

	# eq. 7; tan (gamma_gs) ||rIk[2:3]||_2 <= rIk[1]
	# we model as ||rIk[2:3]||_2 <= (gshelp: rIk[1]/tan(gamma_gs))
	# initalize gshelp_k - 1/tan(gamma_ga) * rIk[1] = 0
	gshelp = MOI.add_variables(model, K)
	MOI.add_constraint(model,
		VAF([VA.(1:K, SA.(1.0, gshelp)); VA.(1:K, SA.(-1.0/tggs, xv[r_idx_it[1], 1:K]))], zeros(K)), MOI.Zeros(K))
	# SECOND ORDER CONE WOOOO; ||rIk[2:3]||_2 <= gshelp_k
	for n=1:K
		MOI.add_constraint(model, VoV([gshelp[n]; xv[r_idx_it[2:3], n]]), SOC(3))
	end

	# eq. 28; cos thmax <= 1 - 2 ||qBIk[3:4]||^2_2
	# or ||qBIk[3:4]||^2_2  <= (1 - cos thmax)/2
	# or ||qBIk[3:4]||_2  <= (aoa_help: sqrt((1 - cos thmax)/2))
	# WE DON'T NEED NO STINKIN' ROTATED CONES
	# initialize aoa_help; need one for each cone :(
	aoa_help = MOI.add_variables(model, K)
	MOI.add_constraint(model, VAF(VA.(1:K, SA.(1.0, aoa_help)), fill(-sqcm, K)), MOI.Zeros(K))
	# ||qBIk[3:4]||_2 <= aoa_help
	for n=1:K
		MOI.add_constraint(model, VoV([aoa_help[n]; xv[qbi_idx_it[3:4], n]]), SOC(3))
	end

	# eq. 10; ||ombk||_2 <= (ang_sp_help: ommax)
	ang_sp_help = MOI.add_variables(model, K)
	MOI.add_constraint(model, VAF(VA.(1:K, SA.(1.0, ang_sp_help)), fill(-prob.omMax, K)), MOI.Zeros(K))
	for n=1:K
		MOI.add_constraint(model, VoV([ang_sp_help[n]; xv[omb_idx_it, n]]), SOC(4))
	end

	# constraints now from Szmuk, Reynolds, and Acikmese 2018 https://arxiv.org/pdf/1811.10803.pdf
	# max thrust constraint ||uk|| <= tMax (8)
	# max_thrust_vk = max_thrust
	#max_thrust = MOI.add_variables(model, 	K+1)
	#MOI.add_constraint(model, VAF(VA.(1:K+1, SA.(1.0, max_thrust)), fill(-prob.Tmax, K+1)), MOI.Zeros(K+1))
	# thrust angle constraint (9)
	# ||uk|| <= uk[3]/(cos delMax)
	# thrust_anglek - uk[3]/(cos delMax) = 0
	#thrust_angle = MOI.add_variables(model, K+1)
	#MOI.add_constraint(model, VAF([VA.(1:K+1, SA.(1.0, thrust_angle)); VA.(1:K+1, SA.(-1/delMax, uv[3,:]))], zeros(K+1)), MOI.Zeros(K+1))

	# combined thrust constraint
	# ||uk|| <= min(max_thrust_vk, thrust_anglek)
	# mtk <= max_thrust_vk /\ mtk <= thrust_anglek
	# ||uk|| <= mtk
	mtk = MOI.add_variables(model, K+1)
	# mtk <= max_thrust_vk <=> mtk - max_thrust_vk <= 0 <=> mtk - prob.Tmax <= 0
	MOI.add_constraint(model, VAF(VA.(1:K+1, SA.(1.0,mtk)), fill(-prob.Tmax, K+1)), MOI.Nonpositives(K+1))
	# mtk <= thrust_anglek <=> mtk - thrust_anglek <= 0 <=> mtk - uk[3]/(cos delMax) <= 0
	MOI.add_constraint(model, VAF([VA.(1:K+1, SA.(1.0,mtk)); VA.(1:K+1, SA.(-1/delMax, uv[1,:]))], zeros(K+1)), MOI.Nonpositives(K+1))

	for n=1:K+1
		MOI.add_constraint(model, VoV([mtk[n]; uv[1:3, n]]), SOC(4))
	end

	#linearized thrust lower bound constraint (8 & 35)
	# Tmin - ||uk|| <= 0
	# h_tlb,k + H_tlb,k duk <= 0
	# h_tlb,k = Tmin - ||uk_lin||
	# H_tlb,k = d(Tmin-||uk_lin||)/duk_lin = -d(uk_lin)/duk_lin = -uk_lin.../norm(uk_lin)
	tlbvars = vcat((VA.(fill(n,3), SA.(-(iterAbout[n].control[1:3] ./ norm(iterAbout[n].control[1:3])), duv[1:3,n])) for n=1:K+1)...)
	tlbconsts = [(prob.Tmin - norm(iterAbout[n].control)) for n=1:K+1]
	thrust_lb_constraint = MOI.add_constraint(model, VAF(tlbvars, tlbconsts), MOI.Nonpositives(K+1))

	# simple fin constraints
	#finmxf = MOI.add_variables(model, K+1)
	#MOI.add_constraint(model, VAF(VA.(1:K+1, SA.(1.0, finmxf)), fill(-0.01,K+1)), MOI.Zeros(K+1))
	for n=1:K+1
		#MOI.add_constraint(model, VAF(VA.(1:2, SA.(1.0, uv[4:5, n])), fill(-0.000001,2)), MOI.Nonpositives(2))
		#MOI.add_constraint(model, VoV([finmxf[n]; uv[4:5, n]]), SOC(3))
	end

	# state triggered constraints
	# todo

	# SCvx trust region
	rK = MOI.add_variable(model)
	rKc = MOI.add_constraint(model, VAF([VA(1, SA(1.0, Jtr))], [-1.0]), MOI.Nonpositives(1))

	return ProblemModel(model, xv, uv, dxv, duv, dsig, nuv, rK, state_base, control_base, dcstrs, thrust_lb_constraint, rKc, (Jtr, ))
end

const MOI=MathOptInterface
function vsq_sub_dotsq(inp::Vector{T}) where T
	return dot(inp[1:3], Dynamics.DCM(inp[4:7])*[1.0,0,0])^2
end

function solve_step(iteration::ProblemIteration, linear_cache::IntegratorCache)
	prob = iteration.problem
	itermodel = iteration.model
	model = itermodel.socp_model
	K = prob.K
	xv = itermodel.xv
	uv = itermodel.uv
	dxv = itermodel.dxv
	duv = itermodel.duv
	dsig = itermodel.dsv
	nuv = itermodel.nuv
	rk = itermodel.rk
	dbg = itermodel.debug
	all_state_dim = state_dim*(K+1)
	all_control_dim = control_dim*(K+1)
	# fix the base constraints

	about = iteration.about
	dynam = iteration.dynam
	MOI.modify(model, itermodel.state_base,
		MOI.VectorConstantChange(Array(vcat(getfield.(about[:], :state)...))))
	MOI.modify(model, itermodel.control_base,
		MOI.VectorConstantChange(Array(vcat(getfield.(about[:], :control)...))))

	# fix the dynamic constraints
	for n=1:K
		vars = [dxv[:,n]; duv[:,n]; duv[:,n+1]; dsig]
		derivs = map(x -> collect(enumerate(x)), eachcol(dynam[n].derivative))
		changes = MOI.MultirowChange.(vars, derivs)
		map(x->MOI.modify(model, itermodel.dynamic_constraints[n], x), changes)
		lin_err = Array(dynam[n].endpoint - about[n+1].state)
		MOI.modify(model, itermodel.dynamic_constraints[n], MOI.VectorConstantChange(lin_err))
	end

	# fix the thrust lower bound
	mrcs = vcat((MOI.MultirowChange.(duv[1:3, n],
				 map(x->[(n, x)], -about[n].control[1:3] ./ norm(about[n].control[1:3]))) for n=1:K+1)...)
	ncsts = [(prob.Tmin - norm(about[n].control[1:3])) for n=1:K+1]
	map(cstr -> MOI.modify(model, itermodel.thrust_lb_constraint, cstr), mrcs)
	MOI.modify(model, itermodel.thrust_lb_constraint, MOI.VectorConstantChange(ncsts))

	#update the trust region
	println("rhk $(iteration.rk)")
	MOI.modify(model, itermodel.rkc, MOI.VectorConstantChange([-iteration.rk]))

	MOI.optimize!(model);

	status = MOI.get(model, MOI.TerminationStatus())
	if status != MOI.OPTIMAL
	    error("Non-optimal result $status exiting")
	end

	xr = reshape(MOI.get(model, MOI.VariablePrimal(), reshape(xv, all_state_dim)), state_dim, K+1)
	ur = reshape(MOI.get(model, MOI.VariablePrimal(), reshape(uv, all_control_dim)), control_dim, K+1)
	dsr = MOI.get(model, MOI.VariablePrimal(), dsig)
	nur = (reshape(MOI.get(model, MOI.VariablePrimal(), reshape(nuv, all_state_dim)), state_dim, K+1))
	dxvr = (reshape(MOI.get(model, MOI.VariablePrimal(), reshape(dxv, all_state_dim)), state_dim, K+1))
	Jtrr = MOI.get(model, MOI.VariablePrimal(), dbg[1])
	println("Jtr $Jtrr")

	info = ProbInfo(prob)
	# Lk = -1 x[1, K+1] + prob.wNu * norm(nuv)
	# Jk = -1 x[1, K+1] + prob.wNu * norm(x[:, i+1] - act[i+1])
	jK = -1 * xr[1, K+1] + prob.wNu * norm(xr[:,k+1] - Dynamics.predict_state(xr[:,k], ur[:,k], ur[:,k+1], iteration.sigma + dsr, 1.0/(K+1), info, linear_cache) for k=1:K)
	lK = -1 * xr[1, K+1] + prob.wNu * norm(nur)

	if iteration.rk == Inf
		next_rk = prob.ri
	else
		jKm = iteration.cost
		djk = jKm - jK
		dlk = jKm - lK
		rhk = djk/dlk
		if rhk < prob.rh0
			println("reject $rhk")
			return ProblemIteration(prob, iteration.cache, iteration.sigma, about, dynam, iteration.model, iteration.iter+1, iteration.rk/prob.alph, iteration.cost), norm(nur), Inf
		elseif rhk < prob.rh1
			next_rk = iteration.rk/prob.alph
			case = 1
		elseif prob.rh1 <= rhk && rhk < prob.rh2
			next_rk = iteration.rk
			case = 2
		else
			next_rk = prob.bet*iteration.rk
			case = 3
		end
		println("nrhk $rhk case $case")
	end

	traj_points = [LinPoint(xr[:,n], ur[:,n]) for n=1:K+1]
	#return traj_points
	nsig = iteration.sigma + dsr
	linpoints = Dynamics.linearize_dynamics(traj_points, iteration.sigma + dsr, 1/(prob.K+1), linear_cache)
	return ProblemIteration(prob, iteration.cache, iteration.sigma + dsr, traj_points, linpoints,
		iteration.model, iteration.iter+1, next_rk, jK), norm(nur), djk

	#=
	MOI.optimize!(model)
	status = MOI.get(model, MOI.TerminationStatus())
	if status != MOI.OPTIMAL
		println(status)
		if status != MOI.ITERATION_LIMIT
			return
		end
	else
		println("optimal found")
	end
	result_status = MOI.get(model, MOI.PrimalStatus())
	if result_status != MOI.FEASIBLE_POINT
		println(result_status)
	    println("Solver ran successfully did not return a feasible point. The problem may be infeasible. Trying to continue.")
	end

	println(MOI.get(model, MOI.VariablePrimal(), rkv))
	xs = reshape(MOI.get(model, MOI.VariablePrimal(), reshape(dxv, length(xv))), state_dim, K+1)
	us = reshape(MOI.get(model, MOI.VariablePrimal(), reshape(duv, length(uv))), control_dim, K+1)
	dxs = reshape(MOI.get(model, MOI.VariablePrimal(), reshape(dxv, length(xv))), state_dim, K+1)
	dus = reshape(MOI.get(model, MOI.VariablePrimal(), reshape(duv, length(uv))), control_dim, K+1)
	nus = reshape(MOI.get(model, MOI.VariablePrimal(), reshape(nuv, length(nuv))), state_dim, K+1)
	dus = reshape(MOI.get(model, MOI.VariablePrimal(), reshape(duv, length(duv))), control_dim, K+1)
	ds = MOI.get(model, MOI.VariablePrimal(), dsig)
	viol = MOI.get(model, MOI.VariablePrimal(), ctrl_viol)
	aviol = MOI.get(model, MOI.VariablePrimal(), aoa_viol)
	traj_points = Array{LinPoint,1}(undef, K+1)
	for k=1:K+1
		traj_points[k] = LinPoint(xs[:,k], us[:,k])
	end
	infos = ProbInfo(prob)
	pv = 0.0
	for k=1:K
		ps = Dynamics.predict_state(xs[:,k], us[:,k], us[:,k+1], sigHat + ds, 1.0/(K+1), infos)
		pv += norm(xs[:,k+1] - ps,1)
	end
	ic = norm(collect(min(norm(us[:,k])-prob.Tmin,0.0) for k=1:K+1),1)
	daoa = 0.0
	aoas = []
	for k=1:K+1
		quat = xs[qbi_idx_it,k]
		q0 = quat[1]
		q1 = quat[2]
		q2 = quat[3]
		q3 = quat[4]
		qv1 = -(1-2*(q2^2 + q3^2))
		qv2 = -(2*(q1 * q2 + q0 * q3))
		qv3 = -(2*(q1 * q3 - q0 * q2))
		cosaoa = min(dot([qv1,qv2,qv3], xs[v_idx_it,k])/(norm(quat)*norm(xs[v_idx_it,k])),1.0)
		push!(aoas, acosd(cosaoa))
		daoa += max(cosd(45.0) - cosaoa, 0)
	end
	println("costs ic:$(norm(viol,1)) daoa:$daoa rp:$(norm(aviol,1)) $(maximum(aoas)) $(iteration.rk)")
	cost = prob.wNu*pv + prob.wTviol*ic + prob.wTviol*daoa

	#linpoints = Dynamics.linearize_dynamics_rk4(traj_points, sigHat + ds, 1.0/(K+1), ProbInfo(prob))
	#linpoints = Dynamics.linearize_dynamics(traj_points, sigHat + ds, 1.0/(K+1), ProbInfo(prob))

	#=
	dstate = sum(norm(linpoints[i].endpoint - linpoints_rk4[i].endpoint) for i=1:K-1)
	djacob = sum(norm(linpoints[i].derivative .- linpoints_rk4[i].derivative) for i=1:K-1)
	println((linpoints[1].derivative .- linpoints_rk4[1].derivative) ./ (linpoints[1].derivative))
	println("dstate:$dstate djacob:$djacob")
	=#

	if iteration.rk == Inf
		next_rk = prob.ri
	else
		jK = iteration.cost
		lk = prob.wNu*sum(norm(nus[:,k],1) for k=1:K+1) + prob.wTviol*norm(viol,1) + prob.wTviol*norm(aviol,1)
		djk = jK - cost
		dlk = jK - lk
		rhk = djk/dlk
		if rhk < prob.rh0
			println("reject $rhk $djk $jK $(sum(norm(dus))) $(iteration.rk)")
			next_rk = 1.0
			return ProblemIteration(prob,  iteration.cache, iteration.sigma, iteration.about, iteration.dynam, iteration.model, iteration.iter+1, iteration.rk/prob.alph, iteration.cost)
		else
			case = 0
			if rhk < prob.rh1
				next_rk = iteration.rk/prob.alph
				case = 1
			elseif prob.rh1 <= rhk && rhk < prob.rh2
				next_rk = iteration.rk
				case = 2
			else
				next_rk = prob.bet*iteration.rk
				case = 3
			end
			println("accept $rhk $cost pc:$dlk $rhk $case")
		end
	end
	linpoints = Dynamics.linearize_dynamics_symb(traj_points, sigHat + ds, iteration.cache)
	return ProblemIteration(prob, iteration.cache, sigHat + ds, traj_points, linpoints, iteration.model, iteration.iter+1, next_rk, cost)
	=#
end

function run_iters(iprob::DescentProblem, niters::Int)
	ip = Rocketland.create_initial(iprob)
	trjs = []
	tfs = []
	for i = 1:niters
		ip = solve_step(ip)
		push!(tfs, ip.sigma)
		push!(trjs, hcat([pt.state[2:4] for pt in ip.about]...))
	end
	return trjs, tfs
end

function solve_problem(iprob::DescentProblem, cache::LinearCache)
    prob = create_initial(iprob, cache)
    cnu = Inf
    cdel = Inf
    iter = 1
    while (iprob.nuTol < cnu || iprob.delTol < cdel) && iter < iprob.imax
        println(cnu, "|", cdel, "|", iprob.nuTol < cnu, "|", iprob.delTol < cdel)
        prob,cnu,cdel = solve_step(prob, cache)
        iter = iter+1
    end
    return prob,cnu,cdel
end

function solve_problem(iprob::DescentProblem)
	return eval(quote
			lc = Dynamics.initalize_cache($iprob);
			Rocketland.solve_problem($iprob, lc)
		end)
end

using Plots
using StaticArrays
function plot_solution(ip::ProblemIteration)
	ys = hcat(([pt.state[2],pt.state[2]] for pt in ip.about)...)'
	xs = hcat(([pt.state[3],pt.state[4]] for pt in ip.about)...)'
	tmax = max(maximum(pt.state[2] for pt in ip.about),
			maximum(pt.state[3] for pt in ip.about))
	pmax = max(maximum(pt.state[2] for pt in ip.about),
			maximum(pt.state[4] for pt in ip.about))
	tmin = min(minimum(pt.state[2] for pt in ip.about),
			minimum(pt.state[3] for pt in ip.about))
	pmin = min(minimum(pt.state[2] for pt in ip.about),
			minimum(pt.state[4] for pt in ip.about))

	thr = [norm(pt.control)/ip.problem.Tmax for pt in ip.about]
	println(thr)
	p=plot(xs,ys, xlims = (tmin,tmax), layout=2, legend=false)
	for pt in ip.about
		dp = dot(Dynamics.DCM(pt.state[qbi_idx_it])*[1.0,0,0], pt.state[v_idx_it])/norm(pt.state[v_idx_it])
		println(dp)
		dv = Dynamics.DCM(pt.state[qbi_idx_it])	*[1.0,0,0]
		xls = [pt.state[3] pt.state[4]; pt.state[3]+dv[2]/3 pt.state[4]+dv[3]/3]
		yls = [pt.state[2] pt.state[2]; pt.state[2]+dv[1]/3 pt.state[2]+dv[1]/3]
		plot!(p,xls,yls, layout=2)
	end
	display(p)
end
export solve_problem, DescentProblem
end

max_stars_repo_name: BenChung/SuccessiveConvexification

max_issues_repo_name: BenChung/SuccessiveConvexification

max_forks_repo_name: BenChung/SuccessiveConvexification


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

polynomial optimization problem
see:
D. Papp and S. Yildiz
Sum-of-squares optimization without semidefinite programming

real polynomials and real-valued complex polynomials

TODO
- implement PSD formulation for complex case
=#

include(joinpath(@__DIR__, "data_real.jl"))
include(joinpath(@__DIR__, "data_complex.jl"))

struct PolyMinNative{T <: Real} <: ExampleInstanceNative{T}
    is_complex::Bool
    interp_vals::Vector
    Ps::Vector{<:Matrix}
    true_min::Real
    use_primal::Bool # solve primal, else solve dual
    use_wsos::Bool # use wsosinterpnonnegative cone, else PSD formulation
end

function PolyMinNative{T}(
    is_complex::Bool,
    poly_name::Symbol,
    halfdeg::Int,
    args...) where {T <: Real}
    R = (is_complex ? Complex{T} : T)
    interp = get_interp_data(R, poly_name, halfdeg)
    return PolyMinNative{T}(is_complex, interp..., args...)
end

function PolyMinNative{T}(
    is_complex::Bool,
    n::Int,
    halfdeg::Int,
    args...) where {T <: Real}
    interp = random_interp_data(T, n, halfdeg)
    return PolyMinNative{T}(is_complex, interp..., args...)
end

function build(inst::PolyMinNative{T}) where {T <: Real}
    if inst.is_complex
        if !inst.use_wsos
            error("PSD formulation is not implemented yet")
        end
        return build_complex(inst)
    else
        if inst.use_primal && !inst.use_wsos
            error("primal psd formulation is not implemented yet")
        end
        return build_real(inst)
    end
end

function build_real(inst::PolyMinNative{T}) where {T <: Real}
    U = length(inst.interp_vals)

    cones = Cones.Cone{T}[]
    if inst.use_wsos
        push!(cones, Cones.WSOSInterpNonnegative{T, T}(U, inst.Ps, use_dual = !inst.use_primal))
    end

    if inst.use_primal
        c = T[-1]
        A = zeros(T, 0, 1)
        G = ones(T, U, 1)
        b = T[]
        h = inst.interp_vals
    else
        c = inst.interp_vals
        A = ones(T, 1, U) # can eliminate constraint and a variable
        b = T[1]
        if inst.use_wsos
            G = -one(T) * I
            h = zeros(T, U)
        else
            svec_lengths = [Cones.svec_length(size(Pk, 2)) for Pk in inst.Ps]
            G = zeros(T, sum(svec_lengths), U)
            offset = 0
            nonneg_cone_size = 0
            for (Pk, dk) in zip(inst.Ps, svec_lengths)
                Lk = size(Pk, 2)
                if dk == 1
                    nonneg_cone_size += 1
                else
                    if nonneg_cone_size > 0
                        push!(cones, Cones.Nonnegative{T}(nonneg_cone_size))
                    end
                    push!(cones, Cones.PosSemidefTri{T, T}(dk))
                end
                l = 1
                for i in 1:Lk, j in 1:i
                    @. @views G[offset + l, :] = -Pk[:, i] * Pk[:, j]
                    l += 1
                end
                @views Cones.scale_svec!(G[offset .+ (1:dk), :], sqrt(T(2)))
                offset += dk
            end
            if nonneg_cone_size > 0
                push!(cones, Cones.Nonnegative{T}(nonneg_cone_size))
            end
            h = zeros(T, size(G, 1))
        end
    end

    model = Models.Model{T}(c, A, b, G, h, cones)
    return model
end

function build_complex(inst::PolyMinNative{T}) where {T <: Real}
    U = length(inst.interp_vals)

    if inst.use_primal
        c = T[-1]
        A = zeros(T, 0, 1)
        b = T[]
        G = ones(T, U, 1)
        h = inst.interp_vals
    else
        c = inst.interp_vals
        A = ones(T, 1, U)
        b = T[1]
        G = -one(T) * I
        h = zeros(T, U)
    end

    cones = Cones.Cone{T}[Cones.WSOSInterpNonnegative{T, Complex{T}}(U, inst.Ps, use_dual = !inst.use_primal)]

    model = Models.Model{T}(c, A, b, G, h, cones)
    return model
end

function test_extra(inst::PolyMinNative{T}, solve_stats::NamedTuple, ::NamedTuple) where T
    @test solve_stats.status == Solvers.Optimal
    if solve_stats.status == Solvers.Optimal && !isnan(inst.true_min)
        # check objective value is correct
        tol = eps(T)^0.1
        true_min = (inst.use_primal ? -1 : 1) * inst.true_min
        @test solve_stats.primal_obj ≈ true_min atol=tol rtol=tol
    end
end


text:
module FloatingTableView

using Blink: Window, body!, active
using TableView

export browse, showtable

# Method for keeping current window active borrowed
# from Plots.jl
mutable struct CurrentWin
    nullablewin::Union{Window, Nothing}
end
const CURRENT_WIN = CurrentWin(nothing)

iswinnull() = CURRENT_WIN.nullablewin === nothing

"""
    current()
Returns the Blink Window object for the current plot
"""
function current()
    if iswinnull()
        error("No current window")
    end
    CURRENT_WIN.nullablewin
end

current(win::Window) = (CURRENT_WIN.nullablewin = win)

"""
	browse(t; kwargs)

Browse data source in a new Blink window using
Tables.jl's `showtable` function.

# Arguments

- `t`, a Tables.jl-compatible data source,
  for example a `DataFrame` or a `NamedTuple`
  of `Vector`s
- `kwargs`, keyword arguments passed to the
  `showtable` command from TableView.jl to
  control the details of the new data viewer.
  See [`showtable`](@ref) for details.

# Examples

```julia
julia> t = (a = rand(10), b = rand(10))

julia> browse(t)
```

```julia
julia> t = rand(100, 100)

julia> browse(t, dark = true)
```

See also: [`showtable`](@ref)
"""
function browse(df; kwargs...)
	if iswinnull() || !active(current())
		w = Window()
		current(w)
		body!(current(), showtable(df; kwargs...))
	else
		body!(current(), showtable(df; kwargs...))
	end

    return nothing
end

end # module


text:
@testset "Embed" begin
    for f ∈ readdir("./embed/")
        include("./embed/$f")
    end
end


text:
# indicator of a generic box

"""
  IndBox(lb, ub)

Returns the function `g = ind{x : lb ⩽ x ⩽ ub}`. Parameters `lb` and `ub` can be
either scalars or arrays of the same dimension as `x`, and must satisfy `lb <= ub`.
Bounds are allowed to take values `-Inf` and `+Inf`.
"""

immutable IndBox{T <: Union{Real, AbstractArray}, S <: Union{Real, AbstractArray}} <: IndicatorConvex
  lb::T
  ub::S
  function IndBox(lb::T, ub::S)
    if !(eltype(lb) <: Real && eltype(ub) <: Real)
      error("lb and ub must be real")
    end
    if any(lb .> ub)
      error("lb and ub must satisfy lb <= ub")
    else
      new(lb, ub)
    end
  end
end

IndBox{T <: Real}(lb::T, ub::T) = IndBox{T, T}(lb, ub)

IndBox{T <: AbstractArray, S <: Real}(lb::T, ub::S) = IndBox{T, S}(lb, ub)

IndBox{T <: Real, S <: AbstractArray}(lb::T, ub::S) = IndBox{T, S}(lb, ub)

IndBox{T <: AbstractArray, S <: AbstractArray}(lb::T, ub::S) =
  size(lb) != size(ub) ? error("bounds must have the same dimensions, or at least one of them be scalar") :
  IndBox{T, S}(lb, ub)

IndBox_lb{T <: Real, S}(f::IndBox{T, S}, i) = f.lb
IndBox_lb{T <: AbstractArray, S}(f::IndBox{T, S}, i) = f.lb[i]
IndBox_ub{T, S <: Real}(f::IndBox{T, S}, i) = f.ub
IndBox_ub{T, S <: AbstractArray}(f::IndBox{T, S}, i) = f.ub[i]

function (f::IndBox){R <: Real}(x::AbstractArray{R})
  for k in eachindex(x)
    if x[k] < IndBox_lb(f,k) || x[k] > IndBox_ub(f,k)
      return +Inf
    end
  end
  return 0.0
end

function prox!{R <: Real}(f::IndBox, x::AbstractArray{R}, y::AbstractArray{R}, gamma::Real=1.0)
  for k in eachindex(x)
    if x[k] < IndBox_lb(f,k)
      y[k] = IndBox_lb(f,k)
    elseif x[k] > IndBox_ub(f,k)
      y[k] = IndBox_ub(f,k)
    else
      y[k] = x[k]
    end
  end
  return 0.0
end

"""
  IndBallLinf(r::Real=1.0)

Returns the indicator function of an infinity-norm ball, that is function
`g(x) = ind{maximum(abs(x)) ⩽ r}` for `r ⩾ 0`.
"""

IndBallLinf{R <: Real}(r::R=1.0) = IndBox(-r, r)

fun_name(f::IndBox) = "indicator of a box"
fun_dom(f::IndBox) = "AbstractArray{Real}"
fun_expr(f::IndBox) = "x ↦ 0 if all(lb ⩽ x ⩽ ub), +∞ otherwise"
fun_params(f::IndBox) =
  string( "lb = ", typeof(f.lb) <: AbstractArray ? string(typeof(f.lb), " of size ", size(f.lb)) : f.lb, ", ",
          "ub = ", typeof(f.ub) <: AbstractArray ? string(typeof(f.ub), " of size ", size(f.ub)) : f.ub)

function prox_naive{R <: Real}(f::IndBox, x::AbstractArray{R}, gamma::Real=1.0)
  y = min.(f.ub, max.(f.lb, x))
  return y, 0.0
end


text:
abstract type AbstractDataContainer end

nSamples(d::AbstractDataContainer) = d.nSamples
nDim(d::AbstractDataContainer) = d.nDim
nOutput(d::AbstractDataContainer) = 1

input(d::AbstractDataContainer) = d.X
output(d::AbstractDataContainer) = d.y

struct DataContainer{Tx<:Real,TX<:AbstractVector,Ty<:Real,TY<:AbstractVector} <:
       AbstractDataContainer
    X::TX # Feature vectors
    y::TY # Output (-1,1 for classification, real for regression, vector{vector} for multiclass)
    nSamples::Int # Number of samples
    nDim::Int # Number of features per sample
end

function wrap_data(X::TX, y::TY) where {TX,TY<:AbstractVector{<:Real}}
    length(y) == length(X) || error(
        "There is not the same number of samples in X ($(length(X))) and y ($(length(y)))",
    )
    Tx = eltype(first(X))
    Ty = eltype(first(y))
    return DataContainer{Tx,TX,Ty,TY}(X, y, length(X), length(first(X)))
end

function wrap_data(X::TX, y::TY) where {TX,TY<:AbstractVector}
    all(length.(y) .== length(X)) || error(
        "There is not the same number of samples in X ($(length(X))) and y ($(length.(y))))",
    )
    Tx = eltype(first(X))
    Ty = eltype(first(y))
    return DataContainer{Tx,TX,Ty,TY}(X, y, length(X), length(first(X)))
end

struct MODataContainer{Tx<:Real,TX<:AbstractVector,TY<:AbstractVector} <:
       AbstractDataContainer
    X::TX # Feature vectors
    y::TY # Output (-1,1 for classification, real for regression, matrix for multiclass)
    nSamples::Int # Number of samples
    nDim::Int # Number of features per sample
    nOutput::Int # Number of outputs
end

function wrap_modata(X::TX, y::TY) where {TX,TY<:AbstractVector}
    all(length.(y) .== length(X)) || error(
        "There is not the same number of samples in X ($(length(X))) and y ($(length.(y))))",
    )
    Tx = eltype(first(X))
    return MODataContainer{Tx,TX,TY}(X, y, length(X), length(first(X)), length(y))
end

nOutput(d::MODataContainer) = d.nOutput

function wrap_X(X::AbstractMatrix{T}, obsdim::Int=1) where {T<:Real}
    return KernelFunctions.vec_of_vecs(X; obsdim=obsdim), T
end

function wrap_X(X::AbstractVector{T}, ::Int=1) where {T<:Real}
    return X, T
end

function wrap_X(X::AbstractVector{<:AbstractVector{T}}, ::Int=1) where {T<:Real}
    return X, T
end

mutable struct OnlineDataContainer <: AbstractDataContainer
    X::AbstractVector # Feature vectors
    y::AbstractVector # Output (-1,1 for classification, real for regression, matrix for multiclass)
    nSamples::Int # Number of samples
    nDim::Int # Number of features per sample
    function OnlineDataContainer()
        return new()
    end
end

function wrap_data!(data::OnlineDataContainer, X::AbstractVector, y::AbstractVector)
    data.X = X
    data.y = y
    data.nSamples = size(X, 1)
    return data.nDim = size(first(X), 1)
end


text:
module RecoverPose
export GEEV4x4Cache
export five_point, five_point_ransac, essential_ransac
export recover_pose, triangulate, ransac
export p3p_select_model, p3p, p3p_ransac, pre_divide_normalize

using Statistics
using Random
using LinearAlgebra
using StaticArrays
using TypedPolynomials
using MultivariatePolynomials
using RowEchelon
import Polynomials

struct GEEV4x4Cache
    A::Matrix{Float64}
    VL::Matrix{Float64}
    VR::Matrix{Float64}
    WR::Vector{Float64}
    WI::Vector{Float64}
    work::Vector{Float64}
end

function GEEV4x4Cache()
    A = Matrix{Float64}(undef, 4, 4)
    VL = Matrix{Float64}(undef, 4, 4)
    VR = Matrix{Float64}(undef, 0, 0)
    WR = Vector{Float64}(undef, 4)
    WI = Vector{Float64}(undef, 4)
    work = Vector{Float64}(undef, 520)
    GEEV4x4Cache(A, VL, VR, WR, WI, work)
end

@inline function __geev_4x4!(cache::GEEV4x4Cache)
    n, A_stride = 4, 4
    jobvl, jobvr = 'V', 'N'
    lwork = LAPACK.BlasInt(520)
    info = Ref{LAPACK.BlasInt}()
    ccall(
        (:dgeev_64_, LAPACK.liblapack), Cvoid,
        (Ref{UInt8}, Ref{UInt8}, Ref{LAPACK.BlasInt}, Ptr{Float64},
         Ref{LAPACK.BlasInt}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64},
         Ref{LAPACK.BlasInt}, Ptr{Float64}, Ref{LAPACK.BlasInt}, Ptr{Float64},
         Ref{LAPACK.BlasInt}, Ptr{LAPACK.BlasInt}, Clong, Clong),
        jobvl, jobvr, n, cache.A, A_stride,
        cache.WR, cache.WI, cache.VL, n, cache.VR, n,
        cache.work, lwork, info, 1, 1)
end

const PPolynomial = Polynomials.Polynomial

include("ransac.jl")
include("five_point/utils.jl")
include("five_point/chirality.jl")
include("five_point/five_point.jl")
include("pnp/p3p.jl")

end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 39a6b3e2-4c50-11eb-16c8-a53147f5e269
md" ## _If Statements_
"

# ╔═╡ 4269dcbe-4c50-11eb-388d-1b32a4f02a1b
md"""Enter the first number:

`num1 = ` $(@bind num1 html"<input type=number >")"""

# ╔═╡ bc59aeca-4c50-11eb-2f6d-d100da19b862
md"""Enter the second number:

`num2 = ` $(@bind num2 html"<input type=number >")"""

# ╔═╡ c70a7e6a-4c50-11eb-1e95-917470deca87
begin
	if num1 <= num2
			lowest = num1
		else
			lowest = num2
		end
		md"The lowest number is $lowest"
end

# ╔═╡ 0a14f9b2-4c51-11eb-10e6-03e0a1047d07
md" ## Question 2
"

# ╔═╡ 17e389b4-4c51-11eb-1c46-c9f83ecb7c17
md"""Enter a value for x:

`x = ` $(@bind x html"<input type=number >")"""

# ╔═╡ 2377fd6c-4c51-11eb-253b-451a7b8ea60a
if x < 4 | x>7
    md"$x is NOT between 4 and 7"
else
	md"$x is between 4 and 7"
end

# ╔═╡ e0d5642e-4c57-11eb-3498-2b4f485f95e6
md" ## Question 3
"

# ╔═╡ f5f19696-4c57-11eb-3f8b-3d1ccb46cc3f
md" Using a function so that the code can be used more universally"

# ╔═╡ 221b0e70-4c58-11eb-1211-97a2e61f5e60
function lower_num(num3, num4)
	if num3 <= num4
			lowest = num3
		else
			lowest = num4
		end
		md"The lowest number is $lowest"
end

# ╔═╡ 318d97bc-4c5e-11eb-18db-118d5869c534
md"""Enter the first number:

`num3 = ` $(@bind num3 html"<input type=number >")"""

# ╔═╡ 3d0a0940-4c5e-11eb-1e57-55302ba763fa
md"""Enter the second number:

`num4 = ` $(@bind num4 html"<input type=number >")"""

# ╔═╡ 63c9b65c-4c5e-11eb-2fea-ad3ca9b2dde3
lower_num(num3,num4)

# ╔═╡ 16d4569c-4dd8-11eb-0048-b16a121fe624
md" ## Challenge 1 & 2: Highest Number + Calling Highest Number"

# ╔═╡ 242fda3c-4dd8-11eb-136d-4d6df172db52
function highest_number(numa, numb, numc)
	if numa >= numb && numa>= numc
		highest = numa
	elseif numb >= numa && numb >= numc
		higest = numb
	else
		highest = numbc
	end
	md"The highest number is $highest"
end

# ╔═╡ 8ba743e6-4dd8-11eb-00d9-257f6b27bfb5
md"""Enter the first number:

`numa = ` $(@bind numa html"<input type=number >")"""

# ╔═╡ 930aba6c-4dd8-11eb-1f9a-85185b0e7033
md"""Enter the second number:

`numb = ` $(@bind numb html"<input type=number >")"""

# ╔═╡ 9a44b4c4-4dd8-11eb-1f72-c36a2a812343
md"""Enter the third number:

`numc = ` $(@bind numc html"<input type=number >")"""

# ╔═╡ a78a9368-4dd8-11eb-0ae1-1185aaa95fcc
highest_number(numa, numb, numc)

# ╔═╡ e2a13c06-4dd8-11eb-1550-55dae34ebc87
md" ## Challenge 3 + Challenge 4: Options and Calling Options"

# ╔═╡ f25fadce-4dd8-11eb-324f-d300053e5be1
function options(value)
	if value == 1
		md"Computer Science"
	elseif value == 2
		md"Music"
	elseif value == 3
		md" Dance"
	elseif value == 4
		md"PE"
	else
		md"Error"
	end
end

# ╔═╡ 6140a984-4dd9-11eb-1e1b-9fa512972905
md"""Enter a value:

`value = ` $(@bind value html"<input type=number >")"""

# ╔═╡ 790a7368-4dd9-11eb-16c7-c5ec34be601c
options(value)

# ╔═╡ a6741d86-4dd9-11eb-29a4-9549137cfa74
md" ## Challenge 6: Refining"

# ╔═╡ 28b74464-4dda-11eb-2453-39a77ae08ecd
function mystery_number(value1)
	if value1 < 3
		md"1"
	elseif value1 == 3
		md"3"
	elseif value1 < 5
		md"8"
	else
		md"$value1"
	end
end

# ╔═╡ 73797900-4dda-11eb-19b1-3592df24113e
md"""Enter a value:

`value1 = ` $(@bind value1 html"<input type=number >")"""

# ╔═╡ 7ef61c52-4dda-11eb-0401-ad5e3405c596
mystery_number(value1)

# ╔═╡ 9448f5cc-4dda-11eb-04d4-c5f044823a51

# ╔═╡ Cell order:
# ╟─39a6b3e2-4c50-11eb-16c8-a53147f5e269
# ╟─4269dcbe-4c50-11eb-388d-1b32a4f02a1b
# ╟─bc59aeca-4c50-11eb-2f6d-d100da19b862
# ╟─c70a7e6a-4c50-11eb-1e95-917470deca87
# ╟─0a14f9b2-4c51-11eb-10e6-03e0a1047d07
# ╟─17e389b4-4c51-11eb-1c46-c9f83ecb7c17
# ╟─2377fd6c-4c51-11eb-253b-451a7b8ea60a
# ╟─e0d5642e-4c57-11eb-3498-2b4f485f95e6
# ╟─f5f19696-4c57-11eb-3f8b-3d1ccb46cc3f
# ╟─221b0e70-4c58-11eb-1211-97a2e61f5e60
# ╟─318d97bc-4c5e-11eb-18db-118d5869c534
# ╟─3d0a0940-4c5e-11eb-1e57-55302ba763fa
# ╠═63c9b65c-4c5e-11eb-2fea-ad3ca9b2dde3
# ╟─16d4569c-4dd8-11eb-0048-b16a121fe624
# ╟─242fda3c-4dd8-11eb-136d-4d6df172db52
# ╟─8ba743e6-4dd8-11eb-00d9-257f6b27bfb5
# ╟─930aba6c-4dd8-11eb-1f9a-85185b0e7033
# ╟─9a44b4c4-4dd8-11eb-1f72-c36a2a812343
# ╠═a78a9368-4dd8-11eb-0ae1-1185aaa95fcc
# ╟─e2a13c06-4dd8-11eb-1550-55dae34ebc87
# ╟─f25fadce-4dd8-11eb-324f-d300053e5be1
# ╟─6140a984-4dd9-11eb-1e1b-9fa512972905
# ╟─790a7368-4dd9-11eb-16c7-c5ec34be601c
# ╟─a6741d86-4dd9-11eb-29a4-9549137cfa74
# ╟─28b74464-4dda-11eb-2453-39a77ae08ecd
# ╟─73797900-4dda-11eb-19b1-3592df24113e
# ╠═7ef61c52-4dda-11eb-0401-ad5e3405c596
# ╠═9448f5cc-4dda-11eb-04d4-c5f044823a51


text:
using FEMSparse, SparseArrays, Test

K = sparse(Float64[1 0 1 1;
                   0 1 0 1;
                   1 0 1 0;
                   1 1 0 1];)

f = zeros(4)
fill!(K, 0.0)

dofs1 = [1, 3]
dofs2 = [2, 4]
dofs3 = [1, 4]
Ke = ones(2, 2)
fe = ones(2)
assembler = FEMSparse.start_assemble(K, f)
for dofs in (dofs1, dofs2, dofs3)
    FEMSparse.assemble_local!(assembler, dofs, Ke, fe)
end

@test Matrix(K) ≈  [2.0  0.0  1.0  1.0;
                    0.0  1.0  0.0  1.0;
                    1.0  0.0  1.0  0.0;
                    1.0  1.0  0.0  2.0]

@test f == [2.0, 1.0, 1.0, 2.0]


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

@doc raw"""Scheduling specifies the scheduling constraints for nodes supporting a RuntimeClass.

    IoK8sApiNodeV1alpha1Scheduling(;
        nodeSelector=nothing,
        tolerations=nothing,
    )

    - nodeSelector::Dict{String, String} : nodeSelector lists labels that must be present on nodes that support this RuntimeClass. Pods using this RuntimeClass can only be scheduled to a node matched by this selector. The RuntimeClass nodeSelector is merged with a pod's existing nodeSelector. Any conflicts will cause the pod to be rejected in admission.
    - tolerations::Vector{IoK8sApiCoreV1Toleration} : tolerations are appended (excluding duplicates) to pods running with this RuntimeClass during admission, effectively unioning the set of nodes tolerated by the pod and the RuntimeClass.
"""
mutable struct IoK8sApiNodeV1alpha1Scheduling <: SwaggerModel
    nodeSelector::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: nodeSelector
    tolerations::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1Toleration} } # spec name: tolerations

    function IoK8sApiNodeV1alpha1Scheduling(;nodeSelector=nothing, tolerations=nothing)
        o = new()
        validate_property(IoK8sApiNodeV1alpha1Scheduling, Symbol("nodeSelector"), nodeSelector)
        setfield!(o, Symbol("nodeSelector"), nodeSelector)
        validate_property(IoK8sApiNodeV1alpha1Scheduling, Symbol("tolerations"), tolerations)
        setfield!(o, Symbol("tolerations"), tolerations)
        o
    end
end # type IoK8sApiNodeV1alpha1Scheduling

const _property_map_IoK8sApiNodeV1alpha1Scheduling = Dict{Symbol,Symbol}(Symbol("nodeSelector")=>Symbol("nodeSelector"), Symbol("tolerations")=>Symbol("tolerations"))
const _property_types_IoK8sApiNodeV1alpha1Scheduling = Dict{Symbol,String}(Symbol("nodeSelector")=>"Dict{String, String}", Symbol("tolerations")=>"Vector{IoK8sApiCoreV1Toleration}")
Base.propertynames(::Type{ IoK8sApiNodeV1alpha1Scheduling }) = collect(keys(_property_map_IoK8sApiNodeV1alpha1Scheduling))
Swagger.property_type(::Type{ IoK8sApiNodeV1alpha1Scheduling }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiNodeV1alpha1Scheduling[name]))}
Swagger.field_name(::Type{ IoK8sApiNodeV1alpha1Scheduling }, property_name::Symbol) =  _property_map_IoK8sApiNodeV1alpha1Scheduling[property_name]

function check_required(o::IoK8sApiNodeV1alpha1Scheduling)
    true
end

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


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

@doc raw"""SupplementalGroupsStrategyOptions defines the strategy type and options used to create the strategy.

    IoK8sApiPolicyV1beta1SupplementalGroupsStrategyOptions(;
        ranges=nothing,
        rule=nothing,
    )

    - ranges::Vector{IoK8sApiPolicyV1beta1IDRange} : ranges are the allowed ranges of supplemental groups.  If you would like to force a single supplemental group then supply a single range with the same start and end. Required for MustRunAs.
    - rule::String : rule is the strategy that will dictate what supplemental groups is used in the SecurityContext.
"""
mutable struct IoK8sApiPolicyV1beta1SupplementalGroupsStrategyOptions <: SwaggerModel
    ranges::Any # spec type: Union{ Nothing, Vector{IoK8sApiPolicyV1beta1IDRange} } # spec name: ranges
    rule::Any # spec type: Union{ Nothing, String } # spec name: rule

    function IoK8sApiPolicyV1beta1SupplementalGroupsStrategyOptions(;ranges=nothing, rule=nothing)
        o = new()
        validate_property(IoK8sApiPolicyV1beta1SupplementalGroupsStrategyOptions, Symbol("ranges"), ranges)
        setfield!(o, Symbol("ranges"), ranges)
        validate_property(IoK8sApiPolicyV1beta1SupplementalGroupsStrategyOptions, Symbol("rule"), rule)
        setfield!(o, Symbol("rule"), rule)
        o
    end
end # type IoK8sApiPolicyV1beta1SupplementalGroupsStrategyOptions

const _property_map_IoK8sApiPolicyV1beta1SupplementalGroupsStrategyOptions = Dict{Symbol,Symbol}(Symbol("ranges")=>Symbol("ranges"), Symbol("rule")=>Symbol("rule"))
const _property_types_IoK8sApiPolicyV1beta1SupplementalGroupsStrategyOptions = Dict{Symbol,String}(Symbol("ranges")=>"Vector{IoK8sApiPolicyV1beta1IDRange}", Symbol("rule")=>"String")
Base.propertynames(::Type{ IoK8sApiPolicyV1beta1SupplementalGroupsStrategyOptions }) = collect(keys(_property_map_IoK8sApiPolicyV1beta1SupplementalGroupsStrategyOptions))
Swagger.property_type(::Type{ IoK8sApiPolicyV1beta1SupplementalGroupsStrategyOptions }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiPolicyV1beta1SupplementalGroupsStrategyOptions[name]))}
Swagger.field_name(::Type{ IoK8sApiPolicyV1beta1SupplementalGroupsStrategyOptions }, property_name::Symbol) =  _property_map_IoK8sApiPolicyV1beta1SupplementalGroupsStrategyOptions[property_name]

function check_required(o::IoK8sApiPolicyV1beta1SupplementalGroupsStrategyOptions)
    true
end

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


text:
module SCP

include("SCP_Exact.jl")
include("SCP_EM.jl")
include("SCP_VB.jl")
include("SCP_MC.jl")
include("BackwardKernel.jl")

using .SCP_Exact, .SCP_EM, .SCP_VB, .SCP_MC, .BackwardKernel

export generate, log_marginal
export standard_EM, dual_EM, standard_VB
export Particle, smc_weight, particle_mcmc
export EventQueue, sum, length, eltype, iterate

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
using ROC
using Statistics

auc_(scores,labels) = AUC(roc(scores, labels))

function R2_(scores, labels)
    dY = scores -labels
    Ysc = labels - mean(labels)*ones(size(labels))
    return 1 - (sum(dY.*dY) / sum(Ysc.*Ysc))
end

"""
#R2_(scores, labels) = 1 - sum((scores-labels).*scores-labels)

function explained_variance_score(y_true, y_pred)
    return(sum((y_pred .- mean(y_true)) .^ 2))
end

function total_variance_score(y_true, y_pred)
    return(sum((y_true .- mean(y_true)) .^ 2))
end

function r2_score(y_true, y_pred)
    #println(y_true)
    #println(y_pred)
    println((explained_variance_score(y_true, y_pred), total_variance_score(y_true, y_pred)))
    return(explained_variance_score(y_true, y_pred) / total_variance_score(y_true, y_pred))
end
"""

max_stars_repo_name: padwulf/MyExample

max_issues_repo_name: padwulf/MyExample

max_forks_repo_name: padwulf/MyExample


text:
# convenience function to calculate the mean-squared error
function mse(arr1::AbstractArray, arr2::AbstractArray)
    @assert length(arr1) == length(arr2)
    N = length(arr1)
    err = 0.0
    for i in 1:N
        err += (arr2[i] - arr1[i])^2
    end
    err /= N
end

mse(X::AbstractArray, thresh=1e-8) = Y::AbstractArray -> begin
    if size(X) != size(Y)
        return false
    end

    return mse(X, Y) < thresh
end

issubtype(T::Type) = x -> typeof(x) <: T
lessthan(rhs) = lhs -> lhs < rhs
greaterthan(rhs) = lhs -> lhs > rhs


text:
import Base.find
#=
find : get the triplet of the input Assoc in three array.  Similar to findnz for sparse.
=#
function find(A::Assoc)
    row, col, val = findnz(A.A)
    n = nnz(A.A)
    #map if the Associative array isn't numerical
    if isa(A.row[1],AbstractString)
        row = [ A.row[row[i]] for i in 1:n]
    end

    if isa(A.col[1],AbstractString)
        col = [ A.col[col[i]] for i in 1:n]
    end

    if isa(A.val[1],AbstractString)
        val = [ A.val[val[i]] for i in 1:n]
    end

    return row, col, val
end

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
export add_ite, add_ite_const, add_eval_const
export add_leq
export add_comp

function add_ite(mgr::Ptr{Manager}, f::Ptr{Node}, g::Ptr{Node}, h::Ptr{Node})
    result = ccall((:Cudd_addIte, libcudd),
        Ptr{Node}, (Ptr{Manager}, Ptr{Node}, Ptr{Node}, Ptr{Node}), mgr, f, g, h)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end

function add_ite_const(mgr::Ptr{Manager}, f::Ptr{Node}, g::Ptr{Node}, h::Ptr{Node})
    result = ccall((:Cudd_addIteConstant, libcudd),
        Ptr{Node}, (Ptr{Manager}, Ptr{Node}, Ptr{Node}, Ptr{Node}), mgr, f, g, h)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end

function add_eval_const(mgr::Ptr{Manager}, f::Ptr{Node}, g::Ptr{Node})
    result = ccall((:Cudd_addEvalConst, libcudd),
        Ptr{Node}, (Ptr{Manager}, Ptr{Node}, Ptr{Node}), mgr, f, g)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end

function add_leq(mgr::Ptr{Manager}, f::Ptr{Node}, g::Ptr{Node})
    result = ccall((:Cudd_addLeq, libcudd),
        Cuint, (Ptr{Manager}, Ptr{Node}, Ptr{Node}), mgr, f, g)
    return result
end

function add_comp(mgr::Ptr{Manager}, f::Ptr{Node})
    result = ccall((:Cudd_addCmpl, libcudd),
        Ptr{Node}, (Ptr{Manager}, Ptr{Node}), mgr, f)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end


text:
module CRBM

include("conversion_functions.jl")
include("update_functions.jl")

using RBM
using ProgressMeter
#using PyPlot

export crbm_control_sample
export crbm_binary_train!

#export from RBM package, so that you don't have to import both
export RBM_t
export CRBM_cfg_t, crbm_create_config
export rbm_copy
export rbm_create,rbm_create_full
export rbm_init_weights_random!
export rbm_init_visible_bias!
export rbm_init_output_bias_random!, rbm_init_hidden_bias_random!
export rbm_write, rbm_read, rbm_read_old

export binarise_matrix, i2b, iv2b, binarise_vector
export unbinarise_matrix, b2i, b2iv, bv2dv
export bv, bin_matrix, bin_vector
export unbin_value, unbin_matrix, unbin_vector
export binary_draw
export up, down, binary_up, binary_down

type CRBM_cfg_t
  use_progress_meter::Bool
  use_pyplot::Bool
  plot_steps::Int64
  batchmode::String
  verbose::Bool
end

function crbm_create_config()
  return CRBM_cfg_t(true, true, 100, "random", false)
end

function crbm_learn_sampling(rbm::RBM_t, y::Array{Float64}, X::Array{Float64})
  Z = binary_up(rbm, y, X)
  for i=1:rbm.uditer-1
    X = binary_down(rbm, Z)
    Z = binary_up(rbm, y, X)
  end
  X = binary_down(rbm, Z)
  Z = binary_up(rbm, y, X)
  return X,Z
end

function crbm_control_sample(rbm::RBM_t, y::Array{Float64}, X::Array{Float64})
  Z = binary_up(rbm, y, X)
  for i=1:rbm.uditer-1
    X = binary_down(rbm, Z)
    Z = binary_up(rbm, y, X)
  end
  X = binary_down(rbm, Z)
  return X,Z
end

function crbm_binary_train!(cfg::CRBM_cfg_t, rbm::RBM_t, S::Matrix{Float64}, A::Matrix{Float64})
  @assert (rbm.dropout >= 0.0 && rbm.dropout <= 1.0) "Dropout must be in [0,1]"

  # TODO more assert statements needed
  N  = ceil(log2(rbm.bins))
  if cfg.verbose == true
    println("binarising sensor data")
  end
  binary_s_matrix = binarise_matrix(S, rbm.bins)
  if cfg.verbose == true
    println("binarising actuator data")
  end
  binary_a_matrix = binarise_matrix(A, rbm.bins)
  ns = size(S)[1]
  nb = rbm.batchsize

  if maximum(rbm.W) == 0.0 && minimum(rbm.W) == 0.0 &&
     maximum(rbm.V) == 0.0 && minimum(rbm.V) == 0.0
     if cfg.verbose == true
       println("Initialising W, V, and c.")
     end
    rbm_init_weights_random!(rbm)
    rbm.c = zeros(rbm.m)
  end

  if cfg.verbose == true
    println("Initialising visible bias")
  end
  rbm_init_visible_bias!(rbm, convert(Array{Int64},binary_s_matrix))

  if cfg.use_progress_meter == true
    pm = Progress(rbm.numepochs, 1, "Training progress:", 20)
  end

  m = ns - rbm.batchsize

  if cfg.verbose == true
    println("Starting learning")
  end
  for t=1:rbm.numepochs
    # extract data batch for current epoch
    start = int(1 + floor(rand() * m)) # 1 to m
    if cfg.batchmode == "sequential"
      r     = [start:start+rbm.batchsize-1]
    else #if cfg.batch_mode == "random"
      r     = rand(1:ns, rbm.batchsize)
    end

    s = binary_s_matrix[r,:]
    a = binary_a_matrix[r,:]

    # generate hidden states given the data
    #= z = binary_up(rbm, s, a) =#
    z = binary_up(rbm, s, a)

    # generate random outputs to start sampler
    A=zeros(size(s)[1], rbm.n)
    for i=1:size(s)[1]
      A[i,:] = transpose(i2b(int(floor(2^rbm.n * rand())), rbm.n ))
    end

    (A, Z) = crbm_learn_sampling(rbm, s, A)

    Eb  = transpose(mean(a,1) - mean(A,1))
    Ec  = transpose(mean(z,1) - mean(Z,1))
    EW  = (z' * a - Z' * A)/nb
    EV  = (z' * s - Z' * s)/nb

    Eb = squeeze(Eb,2)
    Ec = squeeze(Ec,2)

    rbm.b = rbm.b + rbm.alpha * Eb
    rbm.c = rbm.c + rbm.alpha * Ec
    rbm.W = rbm.W + rbm.alpha * EW
    rbm.V = rbm.V + rbm.alpha * EV

    if rbm.momentum > 0.0
      f = rbm.alpha * rbm.momentum
      rbm.b = rbm.b + f * rbm.vb
      rbm.c = rbm.c + f * rbm.vc
      rbm.W = rbm.W + f * rbm.vW
      rbm.V = rbm.V + f * rbm.vV
    end

    if rbm.weightcost > 0.0 # using L2
      f     = (1 - rbm.weightcost)
      rbm.W = rbm.W * f
      rbm.V = rbm.V * f
    end

    rbm.vb = Eb
    rbm.vc = Ec
    rbm.vW = EW
    rbm.vV = EV

    #= if cfg.use_pyplot == true =#
      #= if t % cfg.plot_steps == 0 || t == 0 =#
        #= clf() =#
        #= subplot(121) =#
        #= colorbar(imshow(rbm.W)) =#
        #= subplot(122) =#
        #= colorbar(imshow(rbm.V)) =#
      #= end =#
    #= end =#

    if cfg.use_progress_meter
      next!(pm)
    end
  end # training iteration
end

end


text:
using CSV
using DataFrames
using DecisionTree
using Distributions
using Gurobi
using JuMP
using MLBase
using Plots
using Random
using StatsBase
using StatsPlots
using GradientBoosting

data = CSV.read("test/forestfires.csv", header = false)
ys = convert(Matrix, data)[:,end]
Xs = convert(Matrix, data)[:,1:end-1]

shuffleray = hcat(ys,Xs) #combine into single matrix
shuffleray = shuffleray[shuffle(1:end), :]
y_shuffled = shuffleray[:,1] #seperate shuffled matrix into two again
X_shuffled = shuffleray[:,2:end]

y_trn = y_shuffled[1:388] #about 75%
y_tst = y_shuffled[389:end]

X_trn = X_shuffled[1:388,:]
X_tst = X_shuffled[389:end,:]

gb_is_prediction, gb_models = GradientBoosting.fit(y_trn, X_trn, 0.005, 1, 25) #in-sample
gb_os_prediction = GradientBoosting.predict(y_tst, X_tst, 0.005, gb_models) #out-of-sample

print(1/length(y_tst)*sum((y_tst[i]-gb_os_prediction[i])^2 for i = 1:length(y_tst)),)

max_stars_repo_name: ezuccarelli/GradientBoosting

max_issues_repo_name: ezuccarelli/GradientBoosting


text:
global t = 0

function braitenberg(n_::Ptr{Csize_t}, in_::Ptr{Cdouble},
                     m_::Ptr{Csize_t}, out_::Ptr{Cdouble},
                     l_::Ptr{Csize_t}, str_::Ptr{UInt8})

  n::Int64 = unsafe_load(n_)
  m::Int64 = unsafe_load(m_)
  l::Int64 = unsafe_load(l_)

  input::Vector{Float64}  = pointer_to_array(in_,  (n,))
  output::Vector{Float64} = pointer_to_array(out_, (m,))
  message::Vector{UInt8}  = pointer_to_array(str_, (l,))

  # START your code here

  global t

  t = t + 1
  status = 0

  if t % 1000 == 0
    status = 1
  end

  if t == 10000
    status = 2
  end

  o = [-tanh(1/3*(input[4] + input[5] + input[6])),
       -tanh(1/3*(input[1] + input[2] + input[3]))]

  str = string(input[1], " ", input[2], " ", input[3], "\n", input[4], " ", input[5], " ", input[6])

  # END your code here

  unsafe_copy!(pointer(output),  pointer(o),   m)
  unsafe_copy!(pointer(message), pointer(str), l)

  return convert(Cint, status)
end

function test_function()
  println("test function called")
end

cfunction(braitenberg, Cint, (Ptr{Csize_t}, Ptr{Cdouble}, Ptr{Csize_t}, Ptr{Cdouble}, Ptr{Csize_t}, Ptr{UInt8}))

max_stars_repo_name: kzahedi/YARS

max_issues_repo_name: kzahedi/YARS

max_forks_repo_name: kzahedi/YARS


text:
# ---
# title: Zoom
# cover: zoom.gif
# ---

# Scale without resize

using Augmentor
using ImageShow, ImageCore
using Random

# In the case that only a single Zoom factor is specified, the
# operation will assume that the intention is to Zoom all
# dimensions uniformly by that factor.

img_in = testpattern(RGB, ratio=0.5)

mosaicview(
    img_in,
    augment(img_in, Zoom(1.3)),
    augment(img_in, Zoom(1.3, 1));
    fillvalue=colorant"white", nrow=1, npad=10
)

# It is also possible to pass some abstract vector(s) to the
# constructor, in which case Augmentor will randomly sample one of
# its elements every time the operation is applied.

Random.seed!(1337)
img_out = [augment(img_in, Zoom(0.9:0.05:1.2)) for _ in 1:4]

mosaicview(img_out...; nrow=2)

# ## References

#md # ```@docs
#md # Zoom
#md # ```

## save covers #src
using ImageMagick #src
include(joinpath("..", "assets", "utilities.jl")) #src
cover = make_gif(testpattern(RGB, ratio=0.5), Zoom(0.9:0.1:1.5), 5) #src
ImageMagick.save("zoom.gif", cover; fps=1) #src


text:
using DataStructures
using Random

#############
# PlainVtree
#############

"Root of the plain vtree node hierarchy"
abstract type PlainVtreeNode <: VtreeNode end

mutable struct PlainVtreeLeafNode <: PlainVtreeNode
    var::Var
end

mutable struct PlainVtreeInnerNode <: PlainVtreeNode
    left::PlainVtreeNode
    right::PlainVtreeNode
    variables::Vector{Var}
end

const PlainVtree = AbstractVector{<:PlainVtreeNode}

#####################
# Constructor
#####################

function PlainVtreeInnerNode(left::PlainVtreeNode, right::PlainVtreeNode)
    @assert isempty(intersect(variables(left), variables(right)))
    PlainVtreeInnerNode(left, right, [variables(left); variables(right)])
end

function PlainVtreeLeafNode(vars::Vector{Var})
    @assert length(vars) == 1
    PlainVtreeLeafNode(vars[1])
end

PlainVtreeNode(v::Var) = PlainVtreeLeafNode(v)
PlainVtreeNode(left::PlainVtreeNode, right::PlainVtreeNode) = PlainVtreeInnerNode(left, right)

#####################
# Traits
#####################

@inline NodeType(::PlainVtreeLeafNode) = Leaf()
@inline NodeType(::PlainVtreeInnerNode) = Inner()

#####################
# Methods
#####################

@inline children(n::PlainVtreeInnerNode) = [n.left, n.right]

variables(n::PlainVtreeLeafNode) = [n.var]
variables(n::PlainVtreeInnerNode) = n.variables

import ..Utils.isequal_local
"""
Compare whether two vtree nodes are locally equal (enables `equals` and `equals_unordered` from Utils)
"""
isequal_local(leaf1::PlainVtreeNode, leaf2::PlainVtreeNode)::Bool = false #default
isequal_local(leaf1::PlainVtreeLeafNode, leaf2::PlainVtreeLeafNode)::Bool =
    (leaf1.var == leaf2.var)
isequal_local(inner1::PlainVtreeInnerNode, inner2::PlainVtreeInnerNode)::Bool =
    isequal(variables(inner1), variables(inner2))


text:
file_path = "../../../../../../externals/Core.Math.Data/data/HolisticWare.Core.Sample.Data/CSV/Pejcic_318.csv"

data = open(file_path, "r") do f
    filecontent_as_string_local = read(f,String)
    print(filecontent_as_string_local)
  end
# following line would crash
# error:
# ERROR: LoadError: UndefVarError: filecontent_as_string_local not defined
# print(filecontent_as_string_local)

data = open(file_path, "r") do f
  global filecontent_as_string_global = read(f,String)
end
# no error, because filecontent_as_string_global is global variable (`global` keyword)
print(filecontent_as_string_global)

# "r" for reading
data_read = open(file_path, "r") do file
    # attention that it can be used only once. The second time, without reopening the file,
    # read() would return an empty string
    filecontent_as_string_local = read(file ,String)
    println("filecontent_as_string_local = ")
    println(filecontent_as_string_local)
end

data_read = open(file_path, "r") do file
  # attention that it can be used only once. The second time, without reopening the file,
  # read() would return an empty string
  global filecontent_as_string_global = read(file ,String)
end

println("filecontent_as_string_global = ")
println(filecontent_as_string_global)


text:
## Estimation algorithm WGBS
using CpelNano

# Parameter vector
a = 0.25; b = -50.0; c = 5.0; ϕ = [a,b,c]; L = 200;

# Other parameters
x = CpelNano.RegStruct(); x.L = L; x.Nl = fill(1.0, x.L);
x.ρl = rand(0.001:0.001:0.1, x.L); x.dl = fill(10.0, x.L - 1);
α, β = CpelNano.get_αβ_from_ϕ(ϕ, x);
config = CpelNano.CpelNanoConfig(); config.max_em_iters = 50;
config.max_em_init = 10; config.verbose = true;

# Simulate data
M = 25; pobs = 0.25;
rs = CpelNano.cpel_samp_wgbs(M, α, β, pobs);
rs.L = x.L; rs.Nl = x.Nl; rs.ρl = x.ρl; rs.dl = x.dl;

# Run EM algorithm
CpelNano.get_ϕhat_wgbs!(rs,config); rs.ϕhat
α, β = CpelNano.get_αβ_from_ϕ(rs.ϕhat, rs)
CpelNano.gen_x_mc(α,β)

## EM instance
using CpelNano

# Parameter vector
a = 0.1; b = 0.2; c = 5.0; ϕ = [a,b,c];

# Other parameters
x = CpelNano.RegStruct(); x.L = 40; x.Nl = fill(5.0, x.L); x.ρl = rand(0.1:0.01:0.5, x.L); x.dl = fill(10.0, x.L - 1);
α, β = CpelNano.get_αβ_from_ϕ(ϕ, x); M = 100; pobs = 1.0; μ_m = 40.0; σ_m = 2.0; μ_u = 80.0; σ_u = 2.0;
config = CpelNano.CpelNanoConfig(); config.max_em_iters = 30; config.max_em_init = 1;

# Simulate data
rs = CpelNano.cpel_samp_ont(M, α, β, pobs, μ_m, σ_m, μ_u, σ_u); rs.L = x.L; rs.Nl = x.Nl; rs.ρl = x.ρl; rs.dl = x.dl;

# Run instance
ϕhat, Q = CpelNano.em_inst(rs, zeros(3), config)

## Estimation algorithm
using CpelNano

# Parameter vector
a = 1.0; b = -20.0; c = 12.0; ϕ = [a,b,c]; L = 60;

# Other parameters
x = CpelNano.RegStruct(); x.L = L; x.Nl = fill(1.0, x.L); x.ρl = rand(0.001:0.001:0.1, x.L); x.dl = fill(10.0, x.L - 1);
α, β = CpelNano.get_αβ_from_ϕ(ϕ, x); M = 50; pobs = 0.75; μ_m = 2.0; σ_m = 3.0; μ_u = -2.0; σ_u = 3.0;
config = CpelNano.CpelNanoConfig(); config.max_em_iters = 100; config.max_em_init = 1;config.verbose = true;

# Simulate data
rs = CpelNano.cpel_samp_ont(M, α, β, pobs, μ_m, σ_m, μ_u, σ_u);
rs.L = x.L; rs.Nl = x.Nl; rs.ρl = x.ρl; rs.dl = x.dl;

# Run EM algorithm
CpelNano.get_ϕhat!(rs,config); rs.ϕhat


text:
using PartialRejectionSampling

using Plots
using GraphPlot, Colors

using LightGraphs
const LG = LightGraphs

function plot(
    graph::LG.AbstractGraph,
    dims::Vector{Int}=zeros(Int, 2);
    kwargs...
)
    if any(dims .== 0)
        p = GraphPlot.gplot(graph;
                            nodelabel=LG.vertices(graph), kwargs...)
        return p
    else
        @assert prod(dims) == LG.nv(graph)
        pos = collect(Iterators.product(1:dims[2], 1:dims[1]))[:]
        locs_x, locs_y = map(x->x[1], pos), map(x->x[2], pos)

        p = GraphPlot.gplot(graph, locs_x, reverse(locs_y);
                            nodelabel=LG.vertices(graph), kwargs...)
        return p
    end
end

function plot(
    pp::PRS.Ising,
    state,
    dims=zeros(2);
    kwargs...
)
    c_nodes = ifelse.(
            state .== 1,
            Colors.colorant"gray",
            Colors.colorant"white")

    plot(pp.graph, dims; nodefillc=c_nodes, kwargs...)
end

function plot(
    pp::PRS.HardCoreGraph,
    state,
    dims=zeros(2);
    kwargs...
)
    c_nodes = [Colors.colorant"white" for _ in 1:LG.nv(pp.graph)]
    c_nodes[state] .= Colors.colorant"gray"

    plot(pp.graph, dims; nodefillc=c_nodes, kwargs...)
end

function plot(
    pp::PRS.SinkFreeGraph,
    sample,
    dims=zeros(2);
    kwargs...
)
    c_nodes, c_edges = color_sinks(sample)

    plot(sample, dims; nodefillc=c_nodes, edgestrokec=c_edges, kwargs...)
end

function plot(
    pp::PRS.RootedSpanningForest,
    sample,
    dims=zeros(2);
    kwargs...
)
    c_nodes, c_edges = color_cycles(sample)
    c_nodes[collect(pp.roots)] .= Colors.colorant"lightgreen"

    plot(sample, dims; nodefillc=c_nodes, edgestrokec=c_edges)
end

function color_cycles(
    graph::LG.SimpleDiGraph{T},
    c_node=Colors.colorant"turquoise",
    c_node_cycle=Colors.colorant"red",
    c_edge=Colors.colorant"lightgray",
    c_edge_cycle=Colors.colorant"red"
) where {T}

    edge_map = PRS.edgemap(graph)
    edge_idx = T[]
    nodes = T[]
    for cycle in LG.simplecycles(graph)
        union!(nodes, cycle)
        for (x, y) in zip(cycle, circshift(cycle, -1))
            push!(edge_idx, edge_map[LG.Edge(x, y)])
        end
    end

    c_nodes = fill(c_node, LG.nv(graph))
    c_nodes[nodes] .= c_node_cycle

    c_edges = fill(c_edge, LG.ne(graph))
    c_edges[edge_idx] .= c_edge_cycle

    return c_nodes, c_edges
end

function color_sinks(
    graph::LG.SimpleDiGraph{T},
    c_node=Colors.colorant"turquoise",
    c_node_sink=Colors.colorant"red",
    c_edge=Colors.colorant"lightgray",
    c_edge_sink=Colors.colorant"red"
) where {T}

    nodes = T[]
    edge_idx = T[]
    edge_map = PRS.edgemap(graph)
    for v in PRS.sink_nodes(graph)
        push!(nodes, v)
        for w in LG.inneighbors(graph, v)
            push!(edge_idx, edge_map[LG.Edge(w, v)])
        end
    end

    c_nodes = fill(c_node, LG.nv(graph))
    c_nodes[nodes] .= c_node_sink

    c_edges = fill(c_edge, LG.ne(graph))
    c_edges[edge_idx] .= c_edge_sink

    return c_nodes, c_edges
end

function plot(
    pp::PRS.AbstractSpatialPointProcess,
    points;
    show_center=true,
    radius=0,
    color="white",
    window::Union{Nothing, PRS.AbstractSpatialWindow}=nothing
)
    p = Plots.plot([0], [0],
            label="", legend=false,
            color="white",
            linewidth=0,
            aspect_ratio=:equal,
            grid=false,
            title="")

    for x in (points isa Matrix ? eachcol(points) : points)
        if radius > 0
            plot_disk!(p, x, radius, color)
        end
        if show_center
            Plots.scatter!(p, [x[1]], [x[2]], markersize=2, color=color, grid=false)
        end
    end

    win = isnothing(window) ? PRS.window(pp) : window
    Plots.xlims!(win.c[1], win.c[1] + win.w[1])
    Plots.ylims!(win.c[2], win.c[2] + win.w[win.w isa Real ? 1 : 2])

    return p
end

function plot_disk!(p, center, radius, color=:white)
    θ = LinRange(0, 2π, 15)
    x, y = center[1] .+ radius .* cos.(θ), center[2] .+ radius .* sin.(θ)

    Plots.plot!(p, x, y, seriestype=[:shape], c=color, linewidth=0.5, legend=false, fillapha=0.2)

    return p
end