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

using Markdown
using InteractiveUtils

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

# ╔═╡ 7b8a4360-aaa4-11ec-2efc-b1778f6d3a8c
begin
	using PlotlyBase
	using HypertextLiteral
end

# ╔═╡ 810fb486-10b5-460f-a25a-1a7c9d84e256
using LaTeXStrings

# ╔═╡ 676d318f-b4a4-4949-a5db-1c3a5fd9fa68
using AbstractPlutoDingetjes

# ╔═╡ fc52e423-1370-4ca9-95dc-090815278a4a
using PlutoUI

# ╔═╡ a9e7c04a-2d60-4727-9e24-a6dd9a0ee73c
using Dates

# ╔═╡ 7bd46437-8af0-4a15-87e9-1508869e1600
TableOfContents()

# ╔═╡ 4a18fa5d-7c73-468b-bed2-3acff51e3981
publish_to_js = if is_inside_pluto()
	PlutoRunner.publish_to_js
else
	# @warn "You loaded this package outside of Pluto, this is not the intended behavior and you should use either PlotlyBase or PlotlyJS directly"
	x -> x
end

# ╔═╡ 0ae3f943-4f9b-4cfb-aa76-3bcdc7dc9963
"""
	htl_js(x)
Simple convenience constructor for `HypertextLiteral.JavaScript` objects, renamed and re-exported from HypertextLiteral for convenience in case HypertextLiteral is not explicitly loaded alongisde PlutoPlotly.

See also: [`add_plotly_listeners!`](@ref)
"""
htl_js(x) = HypertextLiteral.JavaScript(x)

# ╔═╡ 90fd960f-65b3-4d8c-b8a8-42d3be8c770f
function Base.show(io::IO, mime::MIME"text/html", s::HypertextLiteral.JavaScript)
if is_inside_pluto()
	show(io, mime, Markdown.MD(Markdown.Code("js",s.content)))
else
	show(io, MIME"text/plain",s)
end
end

# ╔═╡ e9d43bc6-390e-43c3-becb-d1584202da41
# This is only used to simplify debugging (by setting this to false) the plotly internal functions using the developer console
const LOAD_MINIFIED = Ref(true)

# ╔═╡ 16f4b455-086b-4a8b-8767-26fb00a77aad
md"""
# Utility Functions
"""

# ╔═╡ 03bf1bc5-37a9-4b02-bff7-f8b42500c4fc
const JS = HypertextLiteral.JavaScript

# ╔═╡ c3dcd4d9-5e57-4189-a7f8-524afd6db1e6
md"""
## Hyperscript.content hack!
"""

# ╔═╡ 6a4a5cc2-dca5-4f5d-a7e2-9b1f2fbaa406
md"""
## ScriptContents struct
"""

# ╔═╡ 441b20b3-ef9a-4d8a-a6b0-6b6be151a3dd
struct ScriptContents
	vec::Vector{JS}
end

# ╔═╡ 271fd3a7-8347-407d-92d0-2d49758cb3f1
function HypertextLiteral.print_script(io::IO, value::ScriptContents)
	for el ∈ value.vec
		print(io, el.content, '\n')
	end
end

# ╔═╡ 907d51fd-9aaf-43d0-a83b-879cae330a0b
md"""
## Plotly Version
"""

# ╔═╡ 10da78b9-9a67-4cd8-9453-c01ea4baabeb
const PLOTLY_VERSION = Ref("2.11.1")

# ╔═╡ ea88edae-c1a1-4cd3-95da-fd6d5cf337ff
function change_plotly_version(ver::String)
	PLOTLY_VERSION[] = ver
end

# ╔═╡ 5e5cf476-609c-4669-bf5d-c3fc3b75b2fe
get_plotly_version() = PLOTLY_VERSION[]

# ╔═╡ 0c84069c-8362-43e7-8f0f-10f445ddc7fd
check_plotly_version() = @htl """
<script>
	let dv = document.createElement('div')
	let ver = window.Plotly?.version
	if (!ver) {
		dv.innerHTML = "Plotly not loaded!"
		return dv
	}
	if (ver === $(PLOTLY_VERSION[])) {
		dv.innerHTML = ver
	} else {
		dv.innerHTML = "The loaded Plotly version (" + ver + ") is different from the one specified in the package ($(HypertextLiteral.JavaScript(PLOTLY_VERSION[]))), reload the browser page to use the version from PlutoPlotly"
	}
	return dv

</script>
"""

# ╔═╡ 83f2fe95-a6c2-42ec-ac86-72a4b4ec95c3
md"""
## Mathjax hack
"""

# ╔═╡ 4b1688b2-677e-41be-9446-e395925a7311
const FORCE_MATHJAX_LOCAL = Ref(false)

# ╔═╡ 49cf85b1-bf09-49de-8468-4e240dc621fa
"""
	force_pluto_mathjax_local::Bool
	force_pluto_mathjax_local(flag::Bool)::Bool

Returns `true` if the `PlutoPlot` `show` method forces svgs produced by MathJax to be locally cached and `false` otherwise.

The flag can be set at package level by providing the intended boolean value as argument to the function

Local svg caching is used to make mathjax in recent plolty versions (>2.10) work as expected. The default `global` caching in Pluto creates problems with the math display.
"""
force_pluto_mathjax_local() = FORCE_MATHJAX_LOCAL[]

# ╔═╡ 997f1421-b2f7-40c2-bc5b-f8a21cb4b04a
force_pluto_mathjax_local(flag::Bool) = FORCE_MATHJAX_LOCAL[] = flag

# ╔═╡ e6b52b32-def4-4d71-80ca-e43530b1e704
md"""
## Preprocess data
"""

# ╔═╡ 77fe2c5d-f3dd-4779-92a4-e0ceadb639a9
md"""
This function is basically `_json_lower` from PlotlyBase, but we do it directly on the PlutoPlot to avoid the modifying the behavior of `_json_lower` for `Plot` objects (which is required to modify how matrices are passed to `publish_to_js`)
"""

# ╔═╡ 2380a265-700d-4fed-a52e-f6fa1ce41391
# Defaults to JSON.lower for generic non-overloaded types
_preprocess(x) = PlotlyBase.JSON.lower(x)

# ╔═╡ 65638528-ed6b-4e67-9cda-70390db81f64
_preprocess(x::TimeType) = sprint(print, x)

# ╔═╡ bc727ded-8675-420d-806e-0b49357118e5
begin
	_preprocess(x::Union{Bool,String,Number,Nothing,Missing}) = x
	_preprocess(x::Union{Tuple,AbstractArray}) = _preprocess.(x)
	_preprocess(m::AbstractMatrix{<:Number}) = [collect(r) for r ∈ eachcol(m)]
	_preprocess(d::Dict) = Dict{Any,Any}(k => _preprocess(v) for (k, v) in pairs(d))
	_preprocess(a::PlotlyBase.HasFields) = Dict{Any,Any}(k => _preprocess(v) for (k, v) in pairs(a.fields))
	_preprocess(c::PlotlyBase.Cycler) = c.vals
	function _preprocess(c::PlotlyBase.ColorScheme)::Vector{Tuple{Float64,String}}
	    N = length(c.colors)
	    map(ic -> ((ic[1] - 1) / (N - 1), _preprocess(ic[2])), enumerate(c.colors))
	end
end

# ╔═╡ f9c0a331-1f1c-4648-9c24-5e9e16d6be18
_preprocess(t::PlotlyBase.Template) = Dict(
    :data => _preprocess(t.data),
    :layout => _preprocess(t.layout)
)

# ╔═╡ b0d77b4f-da8f-4a0b-a244-043b2e3bdfae
function _preprocess(pc::PlotlyBase.PlotConfig)
    out = Dict{Symbol,Any}()
    for fn in fieldnames(PlotlyBase.PlotConfig)
        field = getfield(pc, fn)
        if !isnothing(field)
            out[fn] = field
        end
    end
    out
end

# ╔═╡ b8e1b177-6686-4b58-8c4c-991d9c148520
# Escape latexstrings
_preprocess(s::LaTeXString) = s.s

# ╔═╡ 4e296bdd-cbd4-4d43-a769-0b4a80d7dec9
md"""
## Unique Counter
"""

# ╔═╡ 628c6e1f-03eb-43a2-8092-a2f61cf6bcbd
function unique_io_counter(io::IO, identifier = "script_id")
	!get(io, :is_pluto, false) && return -1 # We simply return -1 if not inside pluto
	# By default pluto inserts a dict inside the IOContext under key :extra_items. See https://github.com/fonsp/Pluto.jl/blob/10747db7ed512c6b3a9881c5cdb2a4daadea766d/src/runner/PlutoRunner.jl#L786
	dict = io.dict[:extra_items]
	# The dict has key of type Tuple{ObjectID, Int64}, so we a custom key with our custom identifier where we will store the counter
	key = (objectid(identifier), 0)
	counter = get(dict, key, 0) + 1
	# Update the counter on the dict that is shared within this IOContext
	dict[key] = counter
end

# ╔═╡ ebcc9c42-9928-4a20-a307-02ee6ef726d0
# Using the unique_io_counter inside the show3 method allows to have unique counters for plots within a same cell.
# This does not ensure that the same plot object is always given the same unique script id if the plots are added to the cells with `if...end` blocks.
function plotly_script_id(io::IO)
	counter = unique_io_counter(io, "plotly-plot")
	return "plot_$counter"
end

# ╔═╡ fa975cb6-4ec1-419a-bcd6-527c0762a533
md"""
# Plot Wrapper
"""

# ╔═╡ 8b57581f-65b3-4edf-abe3-9dfa4ed82ed5
md"""
We define a wrapper around the PlotlyBase.Plot object
"""

# ╔═╡ f5491a94-5ea8-4459-b6ee-5d37f2ba6188
md"""
## Default script contents
"""

# ╔═╡ 92f5a728-6c57-47b7-9929-e2e19f91da2f
const _default_script_contents = htl_js.([
	"""
	// Flag to check if this cell was  manually ran or reactively ran
	const firstRun = this ? false : true
	const PLOT = this ?? document.createElement("div");
	const parent = currentScript.parentElement
	const isPlutoWrapper = parent.classList.contains('raw-html-wrapper')
	""",
	"""
	if (firstRun) {
		// It seem plot divs would not autosize themself inside flexbox containers without this
		parent.appendChild(PLOT)
	}
	""",
	"""
	// If width is not specified, set it to 100%
	PLOT.style.width = plot_obj.layout.width ? "" : "100%"

	// For the height we have to also put a fixed value in case the plot is put on a non-fixed-size container (like the default wrapper)
	PLOT.style.height = plot_obj.layout.height ? plot_obj.layout.height + "px" :
		(isPlutoWrapper || parent.clientHeight == 0) ? "400px" : "100%"
	""",
	"""

	PLOT.classList.forEach(cn => {
		if (cn !== 'js-plotly-plot' && !custom_classlist.includes(cn)) {
			PLOT.classList.toggle(cn, false)
		}
	})
	for (const className of custom_classlist) {
		PLOT.classList.toggle(className, true)
	}
	""",
	"""

	// Create the resizeObserver to make the plot even more responsive! :magic:
	const resizeObserver = new ResizeObserver(entries => {
		/*
		The addition of the invalid argument `plutoresize` seems to fix the problem with calling `relayout` simply with `{autosize: true}` as update breaking mouse relayout events tracking.
		See https://github.com/plotly/plotly.js/issues/6156 for details
		*/
		Plotly.relayout(PLOT, {autosize: true, plutoresize: true})
	})

	resizeObserver.observe(PLOT)
	""",
	"""

	Plotly.react(PLOT, plot_obj).then(() => {
		// Assign the Plotly event listeners
		for (const [key, listener_vec] of Object.entries(plotly_listeners)) {
			for (const listener of listener_vec) {
				PLOT.on(key, listener)
			}
		}
		// Assign the JS event listeners
		for (const [key, listener_vec] of Object.entries(js_listeners)) {
			for (const listener of listener_vec) {
				PLOT.addEventListener(key, listener)
			}
		}
	}
	)
	""",
	"""

	invalidation.then(() => {
		// Remove all plotly listeners
		PLOT.removeAllListeners()
		// Remove all JS listeners
		for (const [key, listener_vec] of Object.entries(js_listeners)) {
			for (const listener of listener_vec) {
				PLOT.removeEventListener(key, listener)
			}
		}
		// Remove the resizeObserver
		resizeObserver.disconnect()
	})
	""",
])

# ╔═╡ 656e982b-d805-4a75-b3e6-53a4444e5374
md"""
## Struct definition
"""

# ╔═╡ 0f088a21-7d5f-43f7-b99f-688338b61dc6
begin
"""
	PlutoPlot(p::Plot; kwargs...)

A wrapper around `PlotlyBase.Plot` to provide optimized visualization within Pluto notebooks exploiting `@htl` from HypertextLiteral.

# Fields
- `Plot::PlotlyBase.Plot`
- `plotly_listeners::Dict{String, Vector{HypertextLitera.JavaScript}}`
- `classList::Vector{String}`

Once the wrapper has been created, the underlying `Plot` object can be accessed from the `Plot` field of the `PlutoPlot` object.

Custom listeners to [plotly events](https://plotly.com/javascript/plotlyjs-events/) can be added to the `PlutoPlot` as *javascript* functions using the [`add_plotly_listener!`](@ref) function.

Multiple listeners can be associated to each event, and they are executed in the order they are added.

A list of custom CSS classes can be added to the PlutoPlot by using the [`add_class!`](@ref) and [`remove_class!`](@ref) functions.

# Examples
```julia
p = PlutoPlot(Plot(rand(10)))
add_plotly_listener!(p, "plotly_click", "e => console.log(e)")
add_class!(p, "custom_class")
```
"""
Base.@kwdef struct PlutoPlot
	Plot::PlotlyBase.Plot
	plotly_listeners::Dict{String, Vector{JS}} = Dict{String, Vector{JS}}()
	js_listeners::Dict{String, Vector{JS}} = Dict{String, Vector{JS}}()
	classList::Vector{String} = String[]
	script_contents::ScriptContents = ScriptContents(deepcopy(_default_script_contents))
end
PlutoPlot(p::PlotlyBase.Plot; kwargs...) = PlutoPlot(;kwargs..., Plot = p)
end

# ╔═╡ 214cae09-fb98-4ca8-8475-62563e31f665
# This is needed till next version of HypertextLiteral is out (current is 0.9.3). See https://github.com/JuliaPluto/HypertextLiteral.jl/issues/28 for details
HypertextLiteral.content(p::PlutoPlot) = HypertextLiteral.Render(p)

# ╔═╡ 1686debe-6d74-4ec5-bf25-12346c8045c2
function push_script!(p::PlutoPlot, items::Vararg{String,N}) where N
	@nospecialize
	push_script!(p, htl_js.(items)...)
end

# ╔═╡ 64ce91b4-aaa3-45ec-b4d6-f24457167667
function _preprocess(pp::PlutoPlot)
	p = pp.Plot
    out = Dict(
        :data => _preprocess(p.data),
        :layout => _preprocess(p.layout),
        :frames => _preprocess(p.frames),
        :config => _preprocess(p.config)
    )

    if templates.default !== "none" && PlotlyBase._isempty(get(out[:layout], :template, Dict()))
        out[:layout][:template] = _preprocess(templates[templates.default])
    end
    out
end

# ╔═╡ 4ebd0ae4-9f4f-42b2-980e-a25550d01b6b
md"""
## Add plotly listeners
"""

# ╔═╡ de2a547b-3ccd-4f56-96c0-81a7d9b2d272
"""
	add_plotly_listener!(p::PlutoPlot, event_name::String, listener::HypertextLiteral.JavaScript)
	add_plotly_listener!(p::PlutoPlot, event_name::String, listener::String)

Add a custom *javascript* `listener` (to be provided as `String` or directly as `HypertextLiteral.JavaScript`) to the `PlutoPlot` object `p`, and associated to the [plotly event](https://plotly.com/javascript/plotlyjs-events/) specified by `event_name`.

The listeners are added to the HTML plot div after rendering. The div where the plot is inserted can be accessed using the variable named `PLOT` inside the listener code.

# Differences with `add_js_listener!`
This function adds a listener using the plotly internal events via the `on` function. These events differ from the standard javascript ones and provide data specific to the plot.

See also: [`add_js_listener!`](@ref), [`htl_js`](@ref)

# Examples:
```julia
p = PlutoPlot(Plot(rand(10), Layout(uirevision = 1)))
add_plotly_listener!(p, "plotly_relayout", htl_js(\"\"\"
function(e) {

console.log(PLOT) // logs the plot div inside the developer console

}
\"\"\"
```
"""
function add_plotly_listener!(p::PlutoPlot, event_name::String, listener::JS)
	ldict = p.plotly_listeners
	listeners_array = get!(ldict, event_name, JS[])
	push!(listeners_array, listener)
	return p
end;

# ╔═╡ 35e643ab-e3ea-427b-85f2-685b6b6103b8
add_plotly_listener!(p::PlutoPlot, event_name, listener::String) = add_plotly_listener!(p, event_name, htl_js(listener))

# ╔═╡ 3cc48426-1c39-4afe-bf3e-7f4aa2197fca
md"""
## Add JS listeners
"""

# ╔═╡ 67137f84-a284-4615-b7a8-729f0a412939
"""
	add_js_listener!(p::PlutoPlot, event_name::String, listener::HypertextLiteral.JavaScript)
	add_js_listener!(p::PlutoPlot, event_name::String, listener::String)

Add a custom *javascript* `listener` (to be provided as `String` or directly as `HypertextLiteral.JavaScript`) to the `PlutoPlot` object `p`, and associated to the javascript event specified by `event_name`.

The listeners are added to the HTML plot div after rendering. The div where the plot is inserted can be accessed using the variable named `PLOT` inside the listener code.

# Differences with `add_plotly_listener!`
This function adds standard javascript events via the `addEventListener` function. These events differ from the plotly specific events.

See also: [`add_plotly_listener!`](@ref), [`htl_js`](@ref)

# Examples:
```julia
p = PlutoPlot(Plot(rand(10), Layout(uirevision = 1)))
add_js_listener!(p, "mousedown", htl_js(\"\"\"
function(e) {

console.log(PLOT) // logs the plot div inside the developer console when pressing down the mouse

}
\"\"\"
```
"""
function add_js_listener!(p::PlutoPlot, event_name::String, listener::JS)
	ldict = p.js_listeners
	listeners_array = get!(ldict, event_name, JS[])
	push!(listeners_array, listener)
	return p
end;

# ╔═╡ 83906aab-d4ac-4c2b-b7a4-718edb0c2a18
add_js_listener!(p::PlutoPlot, event_name, listener::String) = add_js_listener!(p, event_name, htl_js(listener))

# ╔═╡ 0215aea2-eb79-449e-8dee-a32ca3c5d5f9
md"""
## add_class!
"""

# ╔═╡ 4c6a1004-52ca-40c1-915a-081c0a3c5fbf
"""
	add_class!(p::PlutoPlot, className::String)

Add a CSS class with name `className` to the list of custom classes that are added to the PLOT div when displayed inside Pluto. This can be used to give custom CSS styles to certain plots.

See also: [`remove_class!`](@ref)
"""
function add_class!(p::PlutoPlot, className::String)
	cl = p.classList
	if className ∉ cl
		push!(cl, className)
	end
	return p
end

# ╔═╡ 87af8c6f-3d6d-44c6-87bb-588e01829339
md"""
## remove_class!
"""

# ╔═╡ 50077612-a858-48a0-a187-a9de1489f34f
"""
	remove_class!(p::PlutoPlot, className::String)

Remove a CSS class with name `className` (if present) from the list of custom classes that are added to the PLOT div when displayed inside Pluto. This can be used to give custom CSS styles to certain plots.

See also: [`add_class!`](@ref)
"""
function remove_class!(p::PlutoPlot, className::String)
	cl = p.classList
	idx = findfirst(x -> x === className, cl)
	if idx !== nothing
		deleteat!(cl, idx)
	end
	return p
end

# ╔═╡ 49fe75d5-844d-46d5-a251-7023706a7f92
let
	p = PlutoPlot(Plot())
	add_class!(p, "lol")
	remove_class!(p, "lola")
	p.classList
end

# ╔═╡ 73290c37-481a-4f7d-a92d-038766702890
md"""
# plot Method
"""

# ╔═╡ c3fcc72a-389c-456a-aba5-cbfb4a798c9e
function plot(args...;kwargs...)
	@nospecialize
	PlutoPlot(Plot(args...;kwargs...))
end

# ╔═╡ 18e74f8f-39b6-4c8f-a06f-214d4e9dc6fb
plot(scatter(x = 1:10, y = rand(10)), Layout(title = "TITLE", template = "none"))

# ╔═╡ 8a047414-cd5d-4491-a143-eb30578928ce
md"""
# Show Method
"""

# ╔═╡ f9d1e69f-7a07-486d-b43a-334c1c77790a
function _show(pp::PlutoPlot, script_id = "pluto-plotly-div")
ver = htl_js(PLOTLY_VERSION[])
suffix = htl_js(LOAD_MINIFIED[] ? "min.js" : "js")
@htl """
	<script id=$(script_id)>
		// We start by putting all the variable interpolation here at the beginning

		// Publish the plot object to JS
		let plot_obj = $(publish_to_js(_preprocess(pp)))
		// Get the plotly listeners
		const plotly_listeners = $(pp.plotly_listeners)
		// Get the JS listeners
		const js_listeners = $(pp.js_listeners)
		// Deal with eventual custom classes
		let custom_classlist = $(pp.classList)

		// Load the plotly library
		if (!window.Plotly) {
			const {plotly} = await import('https://cdn.plot.ly/plotly-$(ver).$(suffix)')
		}

		// Check if we have to force local mathjax font cache
		if ($(force_pluto_mathjax_local()) && window?.MathJax?.config?.svg?.fontCache === 'global') {
			window.MathJax.config.svg.fontCache = 'local'
		}

		$(pp.script_contents)

		return PLOT
	</script>
"""
end

# ╔═╡ d42d4694-e05d-4e0e-a198-79a3a5cb688a
function Base.show(io::IO, mime::MIME"text/html", plt::PlutoPlot)
	show(io, mime, _show(plt, plotly_script_id(io)))
	# show(io, mime, _show(plt))
end

# ╔═╡ acba5003-a456-4c1a-a53f-71a3bec30251
md"""
# Tests
"""

# ╔═╡ 0c30855c-6542-4b1a-9427-3a8427e75210
md"""
## Slider + UIRevision
"""

# ╔═╡ de0cb780-ff4e-4236-89c4-4c3163337cfc
@bind clk Clock()

# ╔═╡ dd23fe10-a8d5-461a-85a8-e03468cdcd97
@bind N Slider(50:50:250)
# N =let
# 	clk
# 	rand(50:100)
# end

# ╔═╡ 6da3c910-a350-4a7e-b481-88942c97686b
"""
	push_script!(p::PlutoPlot, items...)
Add script contents contained in collection `items` at the end of the plot show method script.
The `item` must either be a collection of `String` or `HypertextLiteral.JavaScript` elements
"""
function push_script!(p::PlutoPlot, items::Vararg{JS,N}) where N
	@nospecialize
	push!(p.script_contents.vec, items...)
	return p
end

# ╔═╡ fdc22972-b8aa-4202-bb4e-bfff92574814
let
	p = plot(rand(4))
	push_script!(p, "console.log('PUSHED')")
end

# ╔═╡ f6a63433-553c-4857-b767-33465eb22934
let
	p = plot(rand(10,4))
	push_script!(p,htl_js("console.log('MAGIC PUSH SCRIPT')"))
	p
end

# ╔═╡ 8bf75ceb-e4ae-4c6c-8ab0-a81350f19bc7
pp = Plot(scatter3d(x = rand(N), y = rand(N), z = rand(N), mode="markers"), Layout(
	uirevision = 1,
	scene = attr(
		xaxis_range = [-1,2],
		yaxis_range = [-1,2],
		zaxis_range = [-1,2],
		aspectmode = "cube",
	),
	height = 550
	# autosize = true,
));

# ╔═╡ ccf62e33-8fcf-45d9-83ed-c7de80800b76
let
	p = PlutoPlot(pp)
	add_plotly_listener!(p, "plotly_relayout", htl_js("""
	(e) => {

	console.log(e)
	//console.log(PLOT._fullLayout._preGUI)

	var eye = e['scene.camera']?.eye;

    if (eye) {
		console.log('update: ', eye);
	} else {
		console.log(e)
	}
	console.log('div: ',PLOT._fullLayout.scene.camera.eye)
   	console.log('plot_obj: ',plot_obj.layout.scene?.camera?.eye)

}
	"""))
	_show(p)
end

# ╔═╡ 1460ece1-7828-4e93-ac37-e979b874b492
md"""
## @bind click
"""

# ╔═╡ 18c80ea2-0df4-40ea-bd87-f8fee463161e
@bind asdasd let
	p = PlutoPlot(Plot(scatter(y = rand(10), name = "test", showlegend=true)))
	add_plotly_listener!(p,"plotly_click", "
	(e) => {

	console.log(e)
    let dt = e.points[0]
	PLOT.value = [dt.x, dt.y]
	PLOT.dispatchEvent(new CustomEvent('input'))
}
	")
	p
end

# ╔═╡ ce29fa1f-0c52-4d38-acbd-0a96cb3b9ce6
asdasd

# ╔═╡ c3e29c94-941d-4a52-a358-c4ffbfc8cab8
md"""
## @bind filtering
"""

# ╔═╡ b0473b9a-2db5-4d03-8344-b8eaf8428d6c
points = [(rand(),rand()) for _ in 1:10000]

# ╔═╡ 73945da3-af45-41fb-9c5d-6fbba6362256
@bind limits let
	p = Plot(
		scatter(x = first.(points), y = last.(points), mode = "markers")
	)|> PlutoPlot
	add_plotly_listener!(p, "plotly_relayout", "
	e => {
	//console.log(e)
	let layout = PLOT.layout
	let asd = {xaxis: layout.xaxis.range, yaxis: layout.yaxis.range}
	PLOT.value = asd
	PLOT.dispatchEvent(new CustomEvent('input'))
	}
	")
end

# ╔═╡ ea9faecf-ecd7-483b-99ad-ede08ba05383
visible_points = let
	if ismissing(limits)
		points
	else
		xrange = limits["xaxis"]
		yrange = limits["yaxis"]
		func(x,y) = x >= xrange[1] && x <= xrange[2] && y >= yrange[1] && y <= yrange[2]
		filter(x -> func(x...), points)
	end
end

# ╔═╡ 684ef6d7-c1ae-4af3-b1bd-f54bc29d7b53
length(visible_points)

# ╔═╡ f8f7b530-1ded-4ce0-a7d9-a8c92afb95c7
md"""
## Multiple Listeners
"""

# ╔═╡ c3b1a198-ef19-4a54-9c32-d9ea32a63812
let
	p = PlutoPlot(Plot(rand(10), Layout(uirevision = 1)))
	add_plotly_listener!(p, "plotly_relayout", htl_js("""
function(e) {

	console.log('listener 1')

}
	"""))
	add_plotly_listener!(p, "plotly_relayout", htl_js("""
function(e) {

	console.log('listener 2')

}
	"""))
	@htl "$p"
end

# ╔═╡ e9fc2030-c2f0-48e9-a807-424039e796b2
let
	p = PlutoPlot(Plot(rand(10), Layout(uirevision = 1)))
	add_plotly_listener!(p, "plotly_relayout", htl_js("""
function(e) {

	console.log('listener 1')

}
	"""))
	add_plotly_listener!(p, "plotly_relayout", htl_js("""
function(e) {

	console.log('listener 2')

}
	"""))
	p.plotly_listeners
end

# ╔═╡ de101f40-27db-43ea-91ed-238502ceaaf7
md"""
## JS Listener
"""

# ╔═╡ 6c709fa0-7a53-4554-ab2a-d8181267ec93
lololol = 1

# ╔═╡ 671296b9-6743-48d6-9c4d-1beac2b505b5
let
	lololol
	p = PlutoPlot(Plot(rand(10), Layout(uirevision = 1)))
	add_js_listener!(p, "mousedown", htl_js("""
function(e) {

	console.log('MOUSEDOWN!')

}
	"""))
end

# ╔═╡ 6128ff76-3f1f-4144-bb3d-f44678210013
md"""
## flexbox
"""

# ╔═╡ a5823eb2-3aaa-4791-bdc8-196eac2ccf2e
@htl """
<div style='height: 400px; display: flex'>
$(Plot(rand(10)) |> PlutoPlot)

$(Plot(rand(10)) |> PlutoPlot)
</div>
"""

# ╔═╡ aaf0fe61-d5e6-4d93-8a22-7f97f1249b35
md"""
## flexbox + uirevision
"""

# ╔═╡ 6e12592d-01fe-455a-a19c-7544258b9791
voila = 1

# ╔═╡ 36c4a5b1-03f2-4f5f-b9af-822a8f7c8cdf
let
	voila
	@htl """
<div style='height: 550px; display: flex; flex-direction: column;'>
<div style='display: flex; flex: 1 1 0'>
$(Plot(rand(10), Layout(uirevision = 1)) |> PlutoPlot)

$(Plot(rand(10)) |> PlutoPlot)
</div>
<div style='display: flex; flex: 1 1 0'>
$(Plot(rand(10)) |> PlutoPlot)

$(Plot(rand(10)) |> PlutoPlot)
</div>
</div>
"""
end

# ╔═╡ 8b1ab8a6-d2a7-4a15-9690-d83ebaed5c19
html"""
<style>
	.js-plotly-plot {
		flex-grow: 1;
		flex-shrink: 1;
	}
</style>
"""

# ╔═╡ 38a81414-0bcd-4d71-af1d-fe154d2ae09a
md"""
## custom class
"""

# ╔═╡ 2dd5534f-ce46-4770-b0f3-6e16005b3a90
cl = ["test_css_class", "lol"]

# ╔═╡ f69c6955-800c-461e-b464-cab4989913f6
let
	p = PlutoPlot(Plot(rand(10)))
	for cn ∈ cl
		add_class!(p, cn)
	end
	p
end

# ╔═╡ bfe5f717-4702-4316-808a-726fefef9e7e
html"""
<style>
	.test_css_class {
		border: 2px;
		border-style: solid;
	}
</style>
"""

# ╔═╡ cb3f5ee4-5504-4337-8a8d-d45784f54c85
md"""
## Sphere
"""

# ╔═╡ e0271a15-08b5-470f-a2d2-6f064cd3a2b2
function sphere(radius, origin = [0,0,0]; N = 50)
u = range(0, 2π; length = N)
v = range(0, π, length = N)
x = [radius * cos(u) * sin(v) + origin[1] for u ∈ u, v ∈ v]
y = [radius * sin(u) * sin(v) + origin[2] for u ∈ u, v ∈ v]
z = [radius * cos(v) + origin[3] for u ∈ u, v ∈ v]
	surface(;x,y,z)
end

# ╔═╡ cb1f840f-8d99-4076-9554-7d8ba56e9865
@bind M Slider(0:90)

# ╔═╡ 22245242-80a6-4a5b-815e-39b469002f84
let
	s = sphere(1)
	r = 2
	i = 0
	# M = 90
	x = [r*cosd(M)]
	y = [r*sind(M)]
	z = [0]
	sat = scatter3d(;x,y,z, mode = "markers")
	data = [s,sat]
	plot(data, Layout(
		uirevision = 1,
		scene = attr(
		xaxis_range = [-3,3],
		yaxis_range = [-3,3],
		zaxis_range = [-3,3],
		aspectmode = "cube",
	),
	))
end

# ╔═╡ ebdc0ebc-da58-49c8-a992-5924045c2cac
md"""
## Plot Dates
"""

# ╔═╡ 6d74e3fa-1806-40b0-9995-c1555519603d
let
	x = Date.(2000:2020)
	y = rand(length(x))
	plot(scatter(;x,y,name = "values", showlegend = true), Layout(
		template = "none"
	))
end

# ╔═╡ 7054b9ce-cd00-42cf-b81c-bc27b29f4714
md"""
## Heatmap/Contour tests
"""

# ╔═╡ 5b5293bf-81b3-4e80-995a-f15f91971bd4
md"""
This is testing the expected behavior of heatmap and contour plots as per [https://github.com/JuliaPlots/PlotlyJS.jl/issues/355](https://github.com/JuliaPlots/PlotlyJS.jl/issues/355)
"""

# ╔═╡ 3612826c-0af0-4fef-aafe-112a2948f669
let
	plasma = [[0.0, "#0c0786"],
	          [0.1, "#40039c"],
	          [0.2, "#6a00a7"],
	          [0.3, "#8f0da3"],
	          [0.4, "#b02a8f"],
	          [0.5, "#cb4777"],
	          [0.6, "#e06461"],
	          [0.7, "#f2844b"],
	          [0.8, "#fca635"],
	          [0.9, "#fcce25"],
	          [1.0, "#eff821"]]
	plot(heatmap(z=[1 24 30; 28 18 40; 36 54  7], colorscale=plasma), Layout(
		template = "none",
		yaxis = attr(
			scaleanchor = "x",
			scaleratio = 1,
		)
	))
end

# ╔═╡ ce90eb60-98ed-4901-9218-ed8466bb03c7
let
	a, b = -1., 1.
	c, d= 0., 1.
	nx, ny= 200, 100
	hx = (b-a)/(nx-1)
	x = a:hx:b
	hy = (d-c)/(ny-1)
	y=c:hy:d
	X = x' .* ones(length(y))
	Y = ones(length(x))'.*y
	Z= sin.(pi*pi*(X.*X+Y.*Y));
	# print(size(Z))
	plot(contour(x=x, y=y, z=Z), Layout(width=650, height=350, template = "none"))
end

# ╔═╡ 2fa13939-eba2-4d25-b461-56be79fc1db6
md"""
# Re-execute errored cells
"""

# ╔═╡ 5a324fba-1033-4dcf-b10c-1fa4f231355c
md"""
Since we have functions defined at the bottom of the notebook, when first opening the notebook some of the tests above will error and would need to be re-executed after all the function definitions are loaded.

To do this, we put at the bottom of the notebook a javascript function that re-executes all the errored cells
"""

# ╔═╡ 9f2c0123-7e1a-43b7-861a-d059bb28f776
# # Not really needed anymore after putting the tests at the bottom

# @htl """
# <script>
# const jlerrors = document.querySelectorAll('jlerror')
# for (const err of jlerrors) {
# 	const cell = err.closest('pluto-cell')
# 	const runbut = cell.querySelector('button.runcell')
# 	runbut.click()
# }
# </script>
# """

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
AbstractPlutoDingetjes = "6e696c72-6542-2067-7265-42206c756150"
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
HypertextLiteral = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
PlotlyBase = "a03496cd-edff-5a9b-9e67-9cda94a718b5"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
AbstractPlutoDingetjes = "~1.1.4"
HypertextLiteral = "~0.9.3"
LaTeXStrings = "~1.3.0"
PlotlyBase = "~0.8.18"
PlutoUI = "~0.7.37"
"""

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

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

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

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

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

[[ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "12fc73e5e0af68ad3137b886e3f7c1eacfca2640"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.17.1"

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

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

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

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

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

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.3"

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

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

[[PlotlyBase]]
deps = ["ColorSchemes", "Dates", "DelimitedFiles", "DocStringExtensions", "JSON", "LaTeXStrings", "Logging", "Parameters", "Pkg", "REPL", "Requires", "Statistics", "UUIDs"]
git-tree-sha1 = "180d744848ba316a3d0fdf4dbd34b77c7242963a"
uuid = "a03496cd-edff-5a9b-9e67-9cda94a718b5"
version = "0.8.18"

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

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

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

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

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

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# ╔═╡ Cell order:
# ╠═7b8a4360-aaa4-11ec-2efc-b1778f6d3a8c
# ╠═810fb486-10b5-460f-a25a-1a7c9d84e256
# ╠═676d318f-b4a4-4949-a5db-1c3a5fd9fa68
# ╠═fc52e423-1370-4ca9-95dc-090815278a4a
# ╠═a9e7c04a-2d60-4727-9e24-a6dd9a0ee73c
# ╠═7bd46437-8af0-4a15-87e9-1508869e1600
# ╠═4a18fa5d-7c73-468b-bed2-3acff51e3981
# ╠═0ae3f943-4f9b-4cfb-aa76-3bcdc7dc9963
# ╠═90fd960f-65b3-4d8c-b8a8-42d3be8c770f
# ╠═e9d43bc6-390e-43c3-becb-d1584202da41
# ╟─16f4b455-086b-4a8b-8767-26fb00a77aad
# ╠═03bf1bc5-37a9-4b02-bff7-f8b42500c4fc
# ╟─c3dcd4d9-5e57-4189-a7f8-524afd6db1e6
# ╠═214cae09-fb98-4ca8-8475-62563e31f665
# ╠═6a4a5cc2-dca5-4f5d-a7e2-9b1f2fbaa406
# ╠═441b20b3-ef9a-4d8a-a6b0-6b6be151a3dd
# ╠═271fd3a7-8347-407d-92d0-2d49758cb3f1
# ╠═6da3c910-a350-4a7e-b481-88942c97686b
# ╠═1686debe-6d74-4ec5-bf25-12346c8045c2
# ╠═fdc22972-b8aa-4202-bb4e-bfff92574814
# ╟─907d51fd-9aaf-43d0-a83b-879cae330a0b
# ╠═10da78b9-9a67-4cd8-9453-c01ea4baabeb
# ╠═ea88edae-c1a1-4cd3-95da-fd6d5cf337ff
# ╠═5e5cf476-609c-4669-bf5d-c3fc3b75b2fe
# ╠═0c84069c-8362-43e7-8f0f-10f445ddc7fd
# ╟─83f2fe95-a6c2-42ec-ac86-72a4b4ec95c3
# ╠═4b1688b2-677e-41be-9446-e395925a7311
# ╠═49cf85b1-bf09-49de-8468-4e240dc621fa
# ╠═997f1421-b2f7-40c2-bc5b-f8a21cb4b04a
# ╟─e6b52b32-def4-4d71-80ca-e43530b1e704
# ╟─77fe2c5d-f3dd-4779-92a4-e0ceadb639a9
# ╠═2380a265-700d-4fed-a52e-f6fa1ce41391
# ╠═65638528-ed6b-4e67-9cda-70390db81f64
# ╠═bc727ded-8675-420d-806e-0b49357118e5
# ╠═f9c0a331-1f1c-4648-9c24-5e9e16d6be18
# ╠═b0d77b4f-da8f-4a0b-a244-043b2e3bdfae
# ╠═64ce91b4-aaa3-45ec-b4d6-f24457167667
# ╠═b8e1b177-6686-4b58-8c4c-991d9c148520
# ╟─4e296bdd-cbd4-4d43-a769-0b4a80d7dec9
# ╠═628c6e1f-03eb-43a2-8092-a2f61cf6bcbd
# ╠═ebcc9c42-9928-4a20-a307-02ee6ef726d0
# ╟─fa975cb6-4ec1-419a-bcd6-527c0762a533
# ╟─8b57581f-65b3-4edf-abe3-9dfa4ed82ed5
# ╠═f5491a94-5ea8-4459-b6ee-5d37f2ba6188
# ╠═92f5a728-6c57-47b7-9929-e2e19f91da2f
# ╟─656e982b-d805-4a75-b3e6-53a4444e5374
# ╠═0f088a21-7d5f-43f7-b99f-688338b61dc6
# ╟─4ebd0ae4-9f4f-42b2-980e-a25550d01b6b
# ╠═de2a547b-3ccd-4f56-96c0-81a7d9b2d272
# ╠═35e643ab-e3ea-427b-85f2-685b6b6103b8
# ╠═3cc48426-1c39-4afe-bf3e-7f4aa2197fca
# ╠═67137f84-a284-4615-b7a8-729f0a412939
# ╠═83906aab-d4ac-4c2b-b7a4-718edb0c2a18
# ╟─0215aea2-eb79-449e-8dee-a32ca3c5d5f9
# ╠═4c6a1004-52ca-40c1-915a-081c0a3c5fbf
# ╟─87af8c6f-3d6d-44c6-87bb-588e01829339
# ╠═50077612-a858-48a0-a187-a9de1489f34f
# ╠═49fe75d5-844d-46d5-a251-7023706a7f92
# ╠═73290c37-481a-4f7d-a92d-038766702890
# ╠═c3fcc72a-389c-456a-aba5-cbfb4a798c9e
# ╠═f6a63433-553c-4857-b767-33465eb22934
# ╠═18e74f8f-39b6-4c8f-a06f-214d4e9dc6fb
# ╟─8a047414-cd5d-4491-a143-eb30578928ce
# ╠═f9d1e69f-7a07-486d-b43a-334c1c77790a
# ╠═d42d4694-e05d-4e0e-a198-79a3a5cb688a
# ╟─acba5003-a456-4c1a-a53f-71a3bec30251
# ╟─0c30855c-6542-4b1a-9427-3a8427e75210
# ╠═8bf75ceb-e4ae-4c6c-8ab0-a81350f19bc7
# ╠═de0cb780-ff4e-4236-89c4-4c3163337cfc
# ╠═dd23fe10-a8d5-461a-85a8-e03468cdcd97
# ╠═ccf62e33-8fcf-45d9-83ed-c7de80800b76
# ╠═1460ece1-7828-4e93-ac37-e979b874b492
# ╠═18c80ea2-0df4-40ea-bd87-f8fee463161e
# ╠═ce29fa1f-0c52-4d38-acbd-0a96cb3b9ce6
# ╟─c3e29c94-941d-4a52-a358-c4ffbfc8cab8
# ╠═b0473b9a-2db5-4d03-8344-b8eaf8428d6c
# ╠═73945da3-af45-41fb-9c5d-6fbba6362256
# ╠═684ef6d7-c1ae-4af3-b1bd-f54bc29d7b53
# ╠═ea9faecf-ecd7-483b-99ad-ede08ba05383
# ╠═f8f7b530-1ded-4ce0-a7d9-a8c92afb95c7
# ╠═c3b1a198-ef19-4a54-9c32-d9ea32a63812
# ╠═e9fc2030-c2f0-48e9-a807-424039e796b2
# ╟─de101f40-27db-43ea-91ed-238502ceaaf7
# ╠═6c709fa0-7a53-4554-ab2a-d8181267ec93
# ╠═671296b9-6743-48d6-9c4d-1beac2b505b5
# ╟─6128ff76-3f1f-4144-bb3d-f44678210013
# ╠═a5823eb2-3aaa-4791-bdc8-196eac2ccf2e
# ╟─aaf0fe61-d5e6-4d93-8a22-7f97f1249b35
# ╠═6e12592d-01fe-455a-a19c-7544258b9791
# ╠═36c4a5b1-03f2-4f5f-b9af-822a8f7c8cdf
# ╠═8b1ab8a6-d2a7-4a15-9690-d83ebaed5c19
# ╟─38a81414-0bcd-4d71-af1d-fe154d2ae09a
# ╠═2dd5534f-ce46-4770-b0f3-6e16005b3a90
# ╠═f69c6955-800c-461e-b464-cab4989913f6
# ╠═bfe5f717-4702-4316-808a-726fefef9e7e
# ╟─cb3f5ee4-5504-4337-8a8d-d45784f54c85
# ╠═e0271a15-08b5-470f-a2d2-6f064cd3a2b2
# ╠═cb1f840f-8d99-4076-9554-7d8ba56e9865
# ╠═22245242-80a6-4a5b-815e-39b469002f84
# ╟─ebdc0ebc-da58-49c8-a992-5924045c2cac
# ╠═6d74e3fa-1806-40b0-9995-c1555519603d
# ╟─7054b9ce-cd00-42cf-b81c-bc27b29f4714
# ╟─5b5293bf-81b3-4e80-995a-f15f91971bd4
# ╠═3612826c-0af0-4fef-aafe-112a2948f669
# ╠═ce90eb60-98ed-4901-9218-ed8466bb03c7
# ╟─2fa13939-eba2-4d25-b461-56be79fc1db6
# ╟─5a324fba-1033-4dcf-b10c-1fa4f231355c
# ╠═9f2c0123-7e1a-43b7-861a-d059bb28f776
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
module Metamaterials

using Devices, Clipper
import Devices: Meta, GDSMeta, AbstractPolygon, Coordinate

export double_spiral!, rectangular_meander!, interdigit_capacitor!,
    lumped_elem_resonator!, lumped_elem_resonator_with_Ck!, negative,
    cross_polygon

include("double_spiral.jl")
include("interdigit_capacitor.jl")
include("rectangular_meander.jl")

"""
    lumped_elem_resonator!{T,S}(c::Cell{T,S}, Lr_width, Lr_height,
        Cr_overlap, conductor_width, conductor_gap, meta::Meta=GDSMeta(0,0);
        in_length=zero(T), out_length=zero(T))
Create a positive trace of symmetric lumped-element resonator with
specified double spiral coil dimensions (`Lr_width`, `Lr_height`), capacitor
overlap (`Cr_overlap`). The conductor properties are defined with
`conductor_width` and `conductor_gap`. The length of the conductor that goes
in and out of the resonator can be specified by keyword arguments
`in_length` and `out_length`.
"""
function lumped_elem_resonator!(c::Cell{T,S}, Lr_width, Lr_height, Cr_overlap,
    conductor_width, conductor_gap, meta::Meta=GDSMeta(0,0);
    in_length=zero(T), out_length=zero(T)) where {T,S}

    # width and height of the resonator as a whole
    width, height = Lr_width, 2 * Lr_height + Cr_overlap + conductor_gap

    # constants for rectangular coils
    connect_length = Cr_overlap / 2 + conductor_gap

    # constants for interdigitial capacitor
    terminal_width = zero(T)
    terminal_length = Lr_width - conductor_gap - conductor_width

    # cells
    upper_coil = Cell{T,S}(uniquename("upper_coil_Lr"))
    double_spiral!(upper_coil, Lr_width, Lr_height, conductor_width,
        conductor_gap, meta, in_length=connect_length, out_length=out_length)

    lower_coil = Cell{T,S}(uniquename("lower_coil_Lr"))
    double_spiral!(lower_coil, Lr_width, Lr_height, conductor_width,
        conductor_gap, meta, clockwise=false, in_length=connect_length,
        out_length=in_length)

    cap = Cell{T,S}(uniquename("capacitor_Cr"))
    interdigit_capacitor!(cap, terminal_width, terminal_length, Cr_overlap,
        conductor_width, conductor_gap, conductor_gap, meta)

    # rendering cells
    pointoffset = Point(-width, Lr_height + connect_length + in_length)

    trans_upper_coil = Translation(pointoffset...)
    upper_coil_polys = [trans_upper_coil(p) for
        p in polygon.(elements(upper_coil))]

    trans_lower_coil = Translation(pointoffset...) ∘ Rotation(π)
    lower_coil_polys = [trans_lower_coil(p) for
        p in polygon.(elements(lower_coil))]

    capoffset = pointoffset + Point(conductor_width + conductor_gap,
        connect_length)
    trans_cap = Translation(capoffset...) ∘ Rotation(-π/2)
    cap_polys = [trans_cap(p) for p in polygon.(elements(cap))]

    for p in upper_coil_polys
        render!(c, p, meta)
    end
    for p in cap_polys
        render!(c, p, meta)
    end
    for p in lower_coil_polys
        render!(c, p, meta)
    end
    return c
end

    """ lumped_elem_resonator_with_Ck!{T,S}(c::Cell{T,S}, Lr_width, Lr_height,
            Cr_overlap, Ck_overlap, Ck_terminal_width, conductor_width,
            conductor_gap, meta::Meta=GDSMeta(0,0); ground_length=zero(T),
            Ck_terminal_width2=Ck_terminal_width)
    Create a positive trace of symmetric lumped-element resonator with coupling
    capacitor. The coil dimensions (`Lr_width`, `Lr_height`) and capacitor
    overlap (`Cr_overlap`) defines the resonator. The terminal width
    (`Ck_terminal_width`) and the overlap (`Ck_overlap`) characterizes the
    coupling capacitor. The conductor properties are defined with
    `conductor_width` and `conductor_gap`. The width of the terminal from the
    origin and the length of the conductor into the ground can be specified with
    keyword arguments `Ck_terminal_width2` and `ground_length`, respectively.
    """
    function lumped_elem_resonator_with_Ck!(c::Cell{T,S}, Lr_width, Lr_height,
        Cr_overlap, Ck_overlap, Ck_terminal_width, conductor_width,
        conductor_gap, meta::Meta=GDSMeta(0,0); ground_length=zero(T),
        Ck_terminal_width2=Ck_terminal_width) where {T,S}

        # resonator
        res = Cell{T,S}(uniquename("lumped_elem_res_LrCr"))
        lumped_elem_resonator!(res, Lr_width, Lr_height, Cr_overlap,
            conductor_width, conductor_gap, meta,
            in_length=Ck_terminal_width+conductor_gap,
            out_length=ground_length)

        # capacitor
        cap = Cell{T,S}(uniquename("coupling_cap_Ck"))
        Ck_terminal_length = Lr_width - 2 * conductor_width
        interdigit_capacitor!(cap, Ck_terminal_width, Ck_terminal_length,
            Ck_overlap, conductor_width, conductor_gap, conductor_gap, meta,
            terminal_width2=Ck_terminal_width2)

        # offset point
        p_x = -Ck_terminal_length/2
        p_y = Ck_overlap+2*conductor_gap+Ck_terminal_width+Ck_terminal_width2
        pointoffset = Point(p_x, p_y)

        Ck_ref = CellReference(cap, pointoffset, rot=-π/2)
        res_ref = CellReference(res,
        Point(Lr_width/2, Ck_overlap + 2 * conductor_gap + Ck_terminal_width2))

        push!(c.refs, Ck_ref)
        push!(c.refs, res_ref)
        return c
    end

    """
        negative{T,S}(c::Cell{T,S}, mask::AbstractPolygon,
            meta::Meta=GDSMeta(0,0))
    Take whatever is drawn in the cell off from the mask polygon. i.e. create a
    negative cell drawing within the extent of mask. The references of the
    returned cell is discarded and only the flattened cell remains.
    """
    function negative(c::Cell{T,S}, mask::AbstractPolygon{U},
        meta::Meta=GDSMeta(0,0)) where {T,S,U<:Coordinate}
        c = flatten!(c)
        poly = polygon.(elements(c))

        neg = clip(Clipper.ClipTypeDifference, [mask], poly)[1]
        cell = Cell{T,S}(uniquename("negative_$(name(c))"))
        render!(cell, neg, meta)
        return cell
    end
    function negative(c::Cell{T,S}, mask::AbstractVector{U},
        meta::Meta=GDSMeta(0,0)) where {T,S,U<:AbstractPolygon}
        c = flatten!(c)
        poly = polygon.(elements(c))

        neg = clip(Clipper.ClipTypeDifference, mask, poly)[1]
        cell = Cell{T,S}(uniquename("negative_$(name(c))"))
        render!(cell, neg, meta)
        return cell
    end

end # module

ext: jl


text:
export headtohead

using DataFrames, Chain, DataFrameMacros
using Cascadia: nodeText
using TableScraper: scrape_tables
# using Infiltrator
using Tables

function headtohead(name1, name2; verbose=true)
    # get the front page
    front_page_url = "https://www.goratings.org/en/"
    front_page = scrape_tables(front_page_url, identity)[2] |> DataFrame

    player1_url = @chain front_page begin
        @transform :url = front_page_url*:Name[1].attributes["href"][7:end]
        @subset nodeText(:Name) == name1
        select(_, :url).url[1]
    end

    opp1 = @chain player1_url begin
        scrape_tables(_)[2]
    end

    # println("got here")
    # @infiltrate
    opp1 = DataFrame(opp1)
    # println("didn't got here")

    opp1a = @chain opp1 begin
        @subset :Opponent1 == name2
    end

    opp2 = @chain opp1a begin
        @groupby(:Result)
        @combine(:nrow = length(:Opponent1))
    end

    wins = 0
    losses = 0

    lv = opp2.Result .== "Loss"
    if any(lv)
        losses = opp2[lv, :nrow][1]
    end

    wv = opp2.Result .== "Win"
    if any(wv)
        wins = opp2[wv, :nrow][1]
    end

    if verbose
        println("## Head to head: $name1 vs $name2"); println("")
    end
    res = DataFrame(
        # Player1=[string(Player(name1))],
        Player1 = [name1],
        Record = ["$wins : $losses"],
        Player2 = [name2],
        # Player2=[string(Player(name2))],
    )

    res, select(@transform(opp1a, :Player1=name1), :Date, :Player1, :Opponent1=>:Player2, :Color, :Result)
end


text:
typealias ASet{T} Base.AbstractSet{T}

typealias ODict ObjectIdDict

immutable ObjectIdSet{T} <: ASet{T}
  dict::ObjectIdDict
  ObjectIdSet() = new(ObjectIdDict())
end

Base.eltype{T}(::ObjectIdSet{T}) = T

ObjectIdSet() = ObjectIdSet{Any}()

Base.push!{T}(s::ObjectIdSet{T}, x::T) = (s.dict[x] = nothing; s)
Base.delete!{T}(s::ObjectIdSet{T}, x::T) = (delete!(s.dict, x); s)
Base.in(x, s::ObjectIdSet) = haskey(s.dict, x)

(::Type{ObjectIdSet{T}}){T}(xs) = push!(ObjectIdSet{T}(), xs...)

ObjectIdSet(xs) = ObjectIdSet{eltype(xs)}(xs)

Base.collect(s::ObjectIdSet) = collect(keys(s.dict))
Base.similar(s::ObjectIdSet, T::Type) = ObjectIdSet{T}()

@forward ObjectIdSet.dict Base.length

@iter xs::ObjectIdSet -> keys(xs.dict)

typealias OSet ObjectIdSet

immutable ObjectArraySet{T} <: ASet{T}
  xs::Vector{T}
  ObjectArraySet() = new(T[])
end

Base.in{T}(x::T, s::ObjectArraySet{T}) = any(y -> x ≡ y, s.xs)
Base.push!(s::ObjectArraySet, x) = (x ∉ s && push!(s.xs, x); s)

(::Type{ObjectArraySet{T}}){T}(xs) = push!(ObjectArraySet{T}(), xs...)

ObjectArraySet(xs) = ObjectArraySet{eltype(xs)}(xs)

Base.collect(xs::ObjectArraySet) = xs.xs
Base.similar(s::ObjectArraySet, T::Type) = ObjectArraySet{T}()

@forward ObjectArraySet.xs Base.length

@iter xs::ObjectArraySet -> xs.xs

typealias OASet{T} ObjectArraySet{T}


text:
using Test

# 2-category
############

A, B, C, D = Ob(FreeCategory2, :A, :B, :C, :D)
f, g, F, G = [ Hom(sym, A, B) for sym in [:f,:g,:F,:G] ]
h, k, H, K = [ Hom(sym, B, C) for sym in [:h,:k,:H,:K] ]

# Domains and codomains
α, β = Hom2(:α, f, g), Hom2(:β, g, h)
@test dom(α) == f
@test codom(α) == g
@test dom(dom(α)) == A
@test codom(dom(α)) == B
@test dom(compose(α,β)) == f
@test codom(compose(α,β)) == h
@test_throws SyntaxDomainError composeH(α,β)

α, β = Hom2(:α, f, g), Hom2(:β, h, k)
@test dom(composeH(α,β)) == compose(f,h)
@test codom(composeH(α,β)) == compose(g,k)

# Infix notation (Unicode)
α, β = Hom2(:α, f, g), Hom2(:β, g, h)
@test unicode(compose(f,h)) == "f⋅h"
@test unicode(compose(α,β)) == "α⋅β"

α, β = Hom2(:α, f, g), Hom2(:β, h, k)
@test unicode(composeH(α,β)) == "α*β"

# Double category
#################

A, B, C, D, X, Y = Ob(FreeDoubleCategory, :A, :B, :C, :D, :X, :Y)
f, g, h, k = HomH(:f, A, X), HomH(:g, B, Y), HomH(:h, X, C), HomH(:k, Y, D)
l, r, rr = HomV(:ϕ, A, B), HomV(:r, X, Y), HomV(:rr, C, D)
α, β = Hom2(:α, f, g, l, r), Hom2(:β, h, k, r, rr)
@test composeH(α, β) == α⋆β
αβ = composeH(α, β)
@test top(αβ) == composeH(f, h)
@test bottom(αβ) == composeH(g, k)

# Mmonoidal double category
###########################

A, B, C, D = Ob(FreeSymmetricMonoidalDoubleCategory, :A, :B, :C, :D)
E, F, G, H = Ob(FreeSymmetricMonoidalDoubleCategory, :E, :F, :G, :H)
t, b, l, r = HomH(:t, A, B), HomH(:b, C, D), HomV(:l, A, C), HomV(:r, B, D)
t′, b′, l′, r′ = HomH(:t′, E, F), HomH(:b′, G, H), HomV(:l′, E, G), HomV(:r′, F, H)
s = Hom2(:s, t, b, l, r)
s′ = Hom2(:s′, t′, b′, l′, r′)

# Domains and codomains
@test dom(otimes(t,b)) == otimes(dom(t),dom(b))
@test codom(otimes(t,b)) == otimes(codom(t),codom(b))
@test dom(otimes(l,r)) == otimes(dom(l),dom(r))
@test codom(otimes(l,r)) == otimes(codom(l),codom(r))
@test top(otimes(s,s′)) == otimes(top(s),top(s′))
@test bottom(otimes(s,s′)) == otimes(bottom(s),bottom(s′))

@test dom(braidV(A,B)) == otimes(A,B)
@test codom(braidV(A,B)) == otimes(B,A)
@test top(braidH(t,b)) == otimes(t,b)
@test bottom(braidH(t,b)) == otimes(b,t)
@test left(braidH(t,b)) == braidV(A,C)
@test right(braidH(t,b)) == braidV(B,D)
@test σV(A, B) == braidV(A,B)
@test σH(t, b) == braidH(t,b)

# Associativity and unit
I = munit(FreeSymmetricMonoidalDoubleCategory.Ob)
@test otimes(A,I) == A
@test otimes(I,A) == A
@test otimes(otimes(A,B),A) == otimes(A,otimes(B,A))
@test otimes(otimes(t,b),t) == otimes(t,otimes(b,t))


text:
@testset "Kronecker operators" begin
    TAs = ((Float64,  Array),
           (BigFloat, Array))
    if CUDA.has_cuda_gpu()
        TAs = (TAs..., (Float32, CuArray))
    end

    for (T, A) in TAs
        rng = MersenneTwister(37)
        a = adapt(A, rand(rng, T, 3, 2))
        b = adapt(A, rand(rng, T, 4, 5))
        c = adapt(A, rand(rng, T, 1, 7))

        for args in ((a, Eye{T}(5)),
                     (Eye{T}(2), b),
                     (a, b),
                     (Eye{T}(3), Eye{T}(2), c),
                     (Eye{T}(2), b, Eye{T}(7)),
                     (a, Eye{T}(4), Eye{T}(7)),
                     (a, b, c))

            K = adapt(A, collect(Kron(adapt(Array, args)...)))
            d = adapt(A, rand(SVector{2, T}, size(K, 2), 6))
            e = adapt(A, rand(SVector{2, T}, size(K, 2)))
            @test Array(Kron(args...) * e) ≈ Array(K * e)
            @test Array(Kron(args...) * d) ≈ Array(K * d)

            if isbits(T)
                f = rand(rng, T, size(K, 2), 3, 2)
                f = adapt(A, reinterpret(reshape, SVector{2, T},
                                         PermutedDimsArray(f, (3, 1, 2))))
                @test Array(Kron(args...) * f) ≈ Array(K * f)
            end

            @test adapt(Array, Kron(args...)) == Kron(adapt.(Array, args)...)
        end
    end
end


text:
include("load_stuff.jl")
unc,ms,mp,ps,pp=load_data()
function nxt(t,TB,TE)
    if t< TE
        return t+1
    end
    if t>= TE
        return TB
    end
end
function MP!(ms::ms_type,mp::mp_type,unc::u_type,ps::ps_type,pp::pp_type)
    m=Model(Gurobi.Optimizer)

    # */ ------------------ investment variables --------------------------------/* #
    @variable(m, f>=0) # investment-only cost
    @variable(m, x0ph[ms.PH,ms.ZH,ms.I0]>=0)   # newly installed units capacity on energy hub
    @variable(m, x0pp[ms.PP,ms.ZP,ms.I0]>=0)   # newly installed units capacity on production platform
    @variable(m, x0po[ms.PO,ms.ZO,ms.I0]>=0)   # newly installed units capacity onshore
    @variable(m, x0l[ms.L,ms.I0]>=0)    # newly installed line capacity
    @variable(m, x[unc.tx,ms.I]>=0)    # rhs paramters for operational subproblem (e.g., accumulated capacity)
    @variable(m, β[ms.I]>=0)    # operational cost of node i

    # */ ----------------- operational variables --------------------------- /* #
    @variable(m, ϕ[unc.tc,ms.I]>=0)    # operational cost dependent on uncertain cost cᵢ (now is certain)
    @variable(m, c0[ms.I]>=0)  # operational cost independent of certain costs c

    # */ ---------------- electricity variables ---------------------------- /* #
    @variable(m, pD[ps.ZP,ps.sH,ms.I]>=0)            # total power demand on platform z in period t
    @variable(m, pG[ps.G,ps.ZP,ps.sH,ms.I]>=0)       # power output of gas turbine g in period t
    @variable(m, pResG[ps.G,ps.ZP,ps.sH,ms.I]>=0)    # power reserved for spinning reserve of gas turbine g on platform z in period t
    @variable(m, pSEP[ps.Dsep,ps.ZP,ps.sH,ms.I]>=0)  # power required by separator on platform z in period t
    @variable(m, pRW[ps.RW,ps.ZH,ps.sH,ms.I]>=0)     # wind power output on platform z in period t
    @variable(m, pRS[ps.RS,ps.ZH,ps.sH,ms.I]>=0)     # solar power output on platform z in period t
    @variable(m, pSEI[ps.SE,ps.ZP,ps.sH,ms.I]>=0)    # power injection to electricity store e on platform z in period t
    @variable(m, pSEO[ps.SE,ps.ZP,ps.sH,ms.I]>=0)    # power withdrawn from selectricity store e on platform z in period t
    @variable(m, pResSE[ps.SE,ps.ZP,ps.sH,ms.I]>=0)  # power reserved in electricity store e on platform z in period t
    @variable(m, qSE[ps.SE,ps.ZP,ps.sH,ms.I]>=0)     # energy storage level in electricity store e on platform z in period t
    @variable(m, pGShed[ps.Z,ps.sH,ms.I]>=0)         # generation shed on platform z in period t
    @variable(m, pLShed[ps.ZP,ps.sH,ms.I]>=0)        # load shed on platform z in period t
    @variable(m, pZO[ps.B,ps.ZO,ps.sH,ms.I]>=0)      # power output of onshore bus
    @variable(m, pL[ps.L,ps.sH,ms.I])                # power flow in line l in period t

    # */ ------------------- heat variables -------------------------------- /* #
    @variable(m, pEB[ps.EB,ps.ZP,ps.sH,ms.I]>=0) # heat power output from electric boiler on platform z in period t
    @variable(m, pHGShed[ps.ZP,ps.sH,ms.I]>=0)   # heat generation shed on platform z in period t
    @variable(m, pHLShed[ps.ZP,ps.sH,ms.I]>=0)   # heat load shed on platfom z in period t
    @variable(m, pSEPH[ps.Dsep,ps.ZP,ps.sH,ms.I]>=0) # heat power required by separator

    # */ ------------------ natural gas variables -------------------------- /* #
    @variable(m, pC[ps.C,ps.ZP,ps.sH,ms.I]>=0) # power required by recompressor
    @variable(m, pGI[ps.GI,ps.ZP,ps.sH,ms.I]>=0) # power required by gas injection compressor

    # */ ------------------ water variables -------------------------- /* #
    @variable(m, pO[ps.PO,ps.ZP,ps.sH,ms.I]>=0) # power required by oil pump

    # */ ------------------ oil variables -------------------------- /* #
    @variable(m, pWI[ps.PWI,ps.ZP,ps.sH,ms.I]>=0) # power required by water injection pump
    @variable(m, pWL[ps.PWL,ps.ZP,ps.sH,ms.I]>=0) # power required by water lift pump

    # */ ------------------ hydrogen variables ---------------------------- /* #
    @variable(m, pF[ps.F,ps.ZH,ps.sH,ms.I]>=0)   # power output of fuel cell f on offshore energy hub zh in period t
    @variable(m, pE[ps.E,ps.ZH,ps.sH,ms.I]>=0)   # power input to electrolyser on offshore energy hub zh in period t
    @variable(m, vEF[ps.E,ps.ZH,ps.sH,ms.I]>=0)  # hydrogen produced from electrolyser and supply fuel cell on offshore energy hub zh in period t
    @variable(m, vSHy[ps.SHy,ps.ZH,ps.sH,ms.I]>=0)      # storage level of hydrogen storage on offshore energy hub zh in period t
    @variable(m, vSHyI[ps.SHy,ps.ZH,ps.sH,ms.I]>=0)     # hydrogen injection into storage facility on offshore energy hub zh in period t
    @variable(m, vSHyO[ps.SHy,ps.ZH,ps.sH,ms.I]>=0)     # hydrogen withdrawn from the storage facility on offshore energy hub zh in period t
    @variable(m, x0[unc.tx,ms.I])                # values of rhs parameters in the subproblem (e.g., capacity)

    # */ ------------------ objective function ----------------------------- /* #
    @objective(m, Min, f + sum(mp.π[i]*β[i] for i in ms.I) ) # investment plus operational cost (10⁶£)

    # */ ------------------ master constraints ----------------------------- /* #
    @constraint(m, investment_cost, f>=exp10(-6)*sum(mp.π0[i0]*(sum(mp.ci[mp.map_php[p],i0]*x0ph[p,z,i0] for z in ms.ZH for p in ms.PH) # investment cost
                                                    +sum(mp.ci[mp.map_ppp[p],i0]*x0pp[p,z,i0] for z in ms.ZP for p in ms.PP)
                                                    +sum(mp.ci[mp.map_pop[p],i0]*x0po[p,z,i0] for z in ms.ZO for p in ms.PO)
                                                    +sum(mp.ci[mp.map_ll[l],i0]*x0l[l,i0] for l in ms.L)) for i0 in ms.I0)
                                      +exp10(-6)*sum(mp.π[i]*(sum(mp.cf[mp.map_ppp[p]]*x[(p,z),i] for z in ms.ZP for p in ms.PP)    # fixed O&M cost
                                                    +sum(mp.cf[mp.map_php[p]]*x[(p,z),i] for z in ms.ZH for p in ms.PH)
                                                    +sum(mp.cf[mp.map_pop[p]]*x[(p,z),i] for z in ms.ZO for p in ms.PO)
                                                    +sum(mp.cf[mp.map_ppp[p]]*x[(p,z),i] for z in ms.ZP for p in ms.PP)
                                                    +sum(mp.cf[mp.map_ll[l]]*x[(l,l),i] for l in ms.L)) for i in ms.I))
    @constraint(m, accu_cap_zh[p=ms.PH,z=ms.ZH,i=ms.I], x[(p,z),i]==mp.xhph[mp.map_ph[p],mp.map_zh[z],i]+sum(x0ph[p,z,i0] for i0 in mp.map[i]))
    @constraint(m, accu_cap_zp[p=ms.PP,z=ms.ZP,i=ms.I], x[(p,z),i]==mp.xhpp[mp.map_pp[p],mp.map_zp[z],i]+sum(x0pp[p,z,i0] for i0 in mp.map[i]))
    @constraint(m, accu_cap_zo[p=ms.PO,z=ms.ZO,i=ms.I], x[(p,z),i]==mp.xhpo[mp.map_po[p],mp.map_zo[z],i]+sum(x0po[p,z,i0] for i0 in mp.map[i]))
    @constraint(m, accu_cap_line[l=ms.L,i=ms.I], x[(l,l),i] == mp.xhl[mp.map_l[l],i] + sum(x0l[l,i0] for i0 in mp.map[i]) )
    @constraint(m, max_cap_zh[p=ms.PH,z=ms.ZH,i=ms.I], sum(x0ph[p,z,i0] for i0 in mp.map[i])<=mp.xm0[mp.map_php[p]])
    @constraint(m, max_cap_zp[p=ms.PP,z=ms.ZP,i=ms.I], sum(x0pp[p,z,i0] for i0 in mp.map[i])<=mp.xm0[mp.map_ppp[p]])
    @constraint(m, max_cap_zo[p=ms.PO,z=ms.ZO,i=ms.I], sum(x0po[p,z,i0] for i0 in mp.map[i])<=mp.xm0[mp.map_pop[p]])
    @constraint(m, max_cap_line[l=ms.L,i=ms.I],sum(x0l[l,i0] for i0 in mp.map[i])<=mp.xm0[mp.map_ll[l]])
    for i in ms.I for j in unc.th
        fix(x[(j,"-"),i],unc.h[i,unc.map_h[j]];force=true)   # set rhs paramters to dependent on subproblem i (e.g., CO₂ budget parameter)
    end end

    # */ ---------------- subproblems constraints --------------------------- /* #
    @constraint(m, cost_emission[i=ms.I], β[i] >= c0[i] + unc.c[i,unc.map_c["co2_price"]]*ϕ["co2_price",i]) # compute generation cost (10⁶£) of subproblem i
    @constraint(m, cost_operation[i=ms.I], c0[i] == exp10(-6)*(sum(sum((sum((pp.CG+pp.CFuel)/pp.ηG*pG[g,z,h,i] for g in ps.G)+
                                           pp.CLShed*pLShed[z,h,i]+pp.CHLShed*pHLShed[z,h,i])*pp.Ht for h in ps.sH)*pp.WS for z in ps.ZP)+
                                           sum(sum(sum(pp.CZOE[pp.map_zo[z]]*pZO[b,z,h,i]*pp.EG/pp.ηG for b in ps.B)*pp.Ht for h in ps.sH)*pp.WS for z in ps.ZO)))

    # */ ----------------------- CO_2 emission--------------------------- /* #
    @constraint(m, emisison[i=ms.I], ϕ["co2_price",i] == exp10(-6)*sum(sum(sum(pG[g,z,h,i]*pp.EG/pp.ηG for g in ps.G)*pp.Ht for h in ps.sH)*pp.WS for z in ps.ZP))

    # */ ---------------- electricity constraints ---------------------------------- /* #
    @constraint(m, onshore_bus_cap[b=ps.B, z=ps.ZO, s=ps.S, h=ps.Hs[s],i=ms.I], pZO[b,z,h,i] <= x0[(b,z),i])
    @constraint(m, turbine_cap[g=ps.G, z=ps.ZP, s=ps.S, h=ps.Hs[s],i=ms.I], pG[g,z,h,i]+pResG[g,z,h,i] <= x0[(g,z),i])
    @constraint(m, turbine_ramp_up[g=ps.G, z=ps.ZP, s=ps.S, h=ps.Hs[s],i=ms.I], pG[g,z,h,i]+pResG[g,z,h,i]-pG[g,z,nxt(h,ps.Hs[s][1],ps.Hs[s][end]),i]-pResG[g,z,nxt(h,ps.Hs[s][1],ps.Hs[s][end]),i]
                <=pp.GGR*x0[(g,z),i])
    @constraint(m, turbine_ramp_down[g=ps.G, z=ps.ZP, s=ps.S, h=ps.Hs[s],i=ms.I], pG[g,z,h,i]+pResG[g,z,h,i]-pG[g,z,nxt(h,ps.Hs[s][1],ps.Hs[s][end]),i]-pResG[g,z,nxt(h,ps.Hs[s][1],ps.Hs[s][end]),i]
                >=-pp.GGR*x0[(g,z),i])
    @constraint(m, solar[r=ps.RS, z=ps.ZH, s=ps.S, h=ps.Hs[s],i=ms.I], pRS[r,z,h,i]==pp.RPS[h,pp.map_zh[z]]*x0[(r,z),i])
    @constraint(m, wind[r=ps.RW, z=ps.ZH, s=ps.S, h=ps.Hs[s],i=ms.I], pRW[r,z,h,i]==pp.RPW[h,pp.map_zh[z]]*x0[(r,z),i])
    @constraint(m, e_store_balance[se=ps.SE, z=ps.ZP, s=ps.S, h=ps.Hs[s],i=ms.I], qSE[se,z,nxt(h,ps.Hs[s][1],ps.Hs[s][end]),i] == qSE[se,z,h,i]+pp.Ht*(pp.ηSE*pSEI[se,z,h,i]-pSEO[se,z,h,i]))
    @constraint(m, e_store_charge_cap[se=ps.SE, z=ps.ZP, s=ps.S, h=ps.Hs[s],i=ms.I], pSEI[se,z,h,i] <= x0[(se,z),i]*pp.ΗSE)
    @constraint(m, e_store_dicharge_cap[se=ps.SE, z=ps.ZP, s=ps.S, h=ps.Hs[s],i=ms.I], pSEO[se,z,h,i]+pResSE[se,z,h,i] <= x0[(se,z),i]*pp.ΗSE)
    @constraint(m, e_store_energy_cap[se=ps.SE, z=ps.ZP, s=ps.S, h=ps.Hs[s],i=ms.I], qSE[se,z,h,i] <= x0[(se,z),i])
    @constraint(m, e_store_discharge_energy_relation[se=ps.SE, z=ps.ZP, s=ps.S, h=ps.Hs[s],i=ms.I], pp.Ht*(pResSE[se,z,h,i]+pSEO[se,z,h,i]) <= qSE[se,z,h,i])
    @constraint(m, KCL_zp[z=ps.ZP, s=ps.S, h=ps.Hs[s],i=ms.I], sum(pG[g,z,h,i] for g in ps.G) +
                    + sum(pp.ηL[pp.map_l[l]]*pp.AE[pp.map_zpz[z],pp.map_l[l]]*pL[l,h,i] for l in ps.L)+ sum(pSEO[se,z,h,i] for se in ps.SE) + pLShed[z,h,i]
                    == pD[z,h,i] + pGShed[z,h,i] + sum(pSEI[se,z,h,i] for se in ps.SE))
    @constraint(m, KCL_zh[z=ps.ZH, s=ps.S, h=ps.Hs[s],i=ms.I], sum(pRS[r,z,h,i] for r in ps.RS) + sum(pRW[r,z,h,i] for r in ps.RW)
                    + sum(pp.ηL[pp.map_l[l]]*pp.AE[pp.map_zhz[z],pp.map_l[l]]*pL[l,h,i] for l in ps.L) + sum(pF[f,z,h,i] for f in ps.F)
                    == pGShed[z,h,i] + sum(pE[e,z,h,i] for e in ps.E))
    @constraint(m, KCL_zo[z=ps.ZO, s=ps.S, h=ps.Hs[s],i=ms.I], sum(pp.ηL[pp.map_l[l]]*pp.AE[pp.map_zoz[z],pp.map_l[l]]*pL[l,h,i] for l in ps.L) + sum(pZO[b,z,h,i] for b in ps.B)
                    == pGShed[z,h,i])
    @constraint(m, line_cap1[l=ps.L, s=ps.S, h=ps.Hs[s],i=ms.I], pL[l,h,i] <= x0[(l,l),i])
    @constraint(m, line_cap2[l=ps.L, s=ps.S, h=ps.Hs[s],i=ms.I], pL[l,h,i] >= -x0[(l,l),i])
    @constraint(m, reserve[z=ps.ZP, s=ps.S, h=ps.Hs[s],i=ms.I], pp.σRes*pD[z,h,i] <= sum(pResG[g,z,h,i] for g in ps.G) + sum(pResSE[se,z,h,i] for se in ps.SE))
    @constraint(m, electricity_demand[z=ps.ZP, s=ps.S, h=ps.Hs[s],i=ms.I], pD[z,h,i] == sum(pEB[eb,z,h,i] for eb in ps.EB) +
                    sum(pC[c,z,h,i] for c in ps.C) + sum(pO[po,z,h,i] for po in ps.PO) + sum(pGI[gi,z,h,i] for gi in ps.GI) +
                    sum(pWI[pwi,z,h,i] for pwi in ps.PWI) + sum(pWL[pwl,z,h,i] for pwl in ps.PWL))

    # */ ------------------------- heat contraints----------------------------------- /* #
    @constraint(m, E_boiler_cap[eb=ps.EB,z=ps.ZP, s=ps.S, h=ps.Hs[s],i=ms.I], pEB[eb,z,h,i] <= x0[(eb,z),i])
    @constraint(m, heat_balance[z=ps.ZP, s=ps.S, h=ps.Hs[s],i=ms.I], sum(pp.ηGh*pG[g,z,h,i] for g in ps.G) + sum(pp.ηEB*pEB[eb,z,h,i] for eb in ps.EB) + pHLShed[z,h,i]
                        == sum(pSEPH[dsep,z,h,i] for dsep in ps.Dsep) + pHGShed[z,h,i])
    @constraint(m, separator[z=ps.ZP, s=ps.S, h=ps.Hs[s],i=ms.I], sum(pSEPH[dsep,z,h,i] for dsep in ps.Dsep) == pp.ρSEP*(pp.VOD[h,pp.map_zp[z]])*x0[("demand_scaling","-"),i])

    # */ ------------------------- natural gas contraints---------------------------- /* #
    @constraint(m, comrepssor[z=ps.ZP, s=ps.S, h=ps.Hs[s],i=ms.I], sum(pC[c,z,h,i] for c in ps.C) == (max((pp.VD[h,pp.map_zp[z]]-pp.VGI[h,pp.map_zp[z]]),0)*0.0051*pp.θNG))
    @constraint(m, inj_compressor[z=ps.ZP, s=ps.S, h=ps.Hs[s],i=ms.I], sum(pGI[gi,z,h,i] for gi in ps.GI) == pp.VGI[h,pp.map_zp[z]]*0.011*pp.θNG)

    # */ ------------------------- oil contraints----------------------------------- /* #
    @constraint(m, power_oil_pump[z=ps.ZP, s=ps.S, h=ps.Hs[s],i=ms.I], sum(pO[po,z,h,i] for po in ps.PO) == pp.κPO*pp.VOD[h,pp.map_zp[z]]*x0[("demand_scaling","-"),i])

    # */ ------------------------- water contraints----------------------------------- /* #
    @constraint(m, power_water_inj[z=ps.ZP, s=ps.S, h=ps.Hs[s],i=ms.I], sum(pWI[pwi,z,h,i] for pwi in ps.PWI) == pp.κPWI*pp.VWI[h,pp.map_zp[z]])
    @constraint(m, power_water_lift[z=ps.ZP, s=ps.S, h=ps.Hs[s],i=ms.I], sum(pWL[pwl,z,h,i] for pwl in ps.PWL) == pp.κPWL*pp.VWL[h,pp.map_zp[z]])

    # */ ------------------------- hydrogen contraints--------------------------------/* #
    @constraint(m, electrolyser_cap[e=ps.E, z=ps.ZH, s=ps.S, h=ps.Hs[s],i=ms.I], pE[e,z,h,i] <= x0[(e,z),i])
    @constraint(m, fuel_cell_cap[f=ps.F, z=ps.ZH, s=ps.S, h=ps.Hs[s],i=ms.I], pF[f,z,h,i] <= x0[(f,z),i])
    @constraint(m, hydrogen_strore_cap[shy=ps.SHy,z=ps.ZH, s=ps.S, h=ps.Hs[s],i=ms.I], vSHy[shy,z,h,i] <= x0[(shy,z),i])
    @constraint(m, fuel_cell_ramp_up[f=ps.F,z=ps.ZH, s=ps.S, h=ps.Hs[s],i=ms.I], pF[f,z,h,i]-pF[f,z,nxt(h,ps.Hs[s][1],ps.Hs[s][end]),i]<=pp.FR*x0[(f,z),i])
    @constraint(m, fuel_cell_ramp_down[f=ps.F,z=ps.ZH, s=ps.S, h=ps.Hs[s],i=ms.I], pF[f,z,h,i]-pF[f,z,nxt(h,ps.Hs[s][1],ps.Hs[s][end]),i]>=-pp.FR*x0[(f,z),i])
    @constraint(m, hydrogen_store_balance[shy=ps.SHy,z=ps.ZH, s=ps.S, h=ps.Hs[s],i=ms.I], vSHy[shy,z,nxt(h,ps.Hs[s][1],ps.Hs[s][end]),i]- vSHy[shy,z,h,i]==
                    vSHyI[shy,z,h,i]-vSHyO[shy,z,h,i])
    @constraint(m, hydrogen_balance1[z=ps.ZH, s=ps.S, h=ps.Hs[s],i=ms.I], pp.Ht*sum(pE[e,z,h,i] for e in ps.E) == pp.ηES*sum(vSHyI[shy,z,h,i] for shy in ps.SHy) + sum(vEF[e,z,h,i] for e in ps.E)*pp.ηEF)
    @constraint(m, hydrogen_balance2[z=ps.ZH, s=ps.S, h=ps.Hs[s],i=ms.I], sum(pp.Ht*pF[f,z,h,i]/(pp.ηF*pp.θHY) for f in ps.F) == sum(vEF[e,z,h,i] for e in ps.E) + sum(vSHyO[shy,z,h,i] for shy in ps.SHy))

    # */ ----------------------------- CO2 budget ------------------------------------ /* #
    @constraint(m, co2_budget[i=ms.I], (sum(pG[g,z,h,i]*pp.EG/pp.ηG for g in ps.G for h in ps.sH for z in ps.ZP))*pp.Ht*pp.WS <=
                        pp.ΜE)

    @constraint(m, pass_cap[n=1:length(unc.tx)-2,i=ms.I], x0[unc.tx[n],i]*exp10(-3) == x[unc.tx[n],i])
    @constraint(m, pass_limits[n=unc.th,i=ms.I], x0[(n,"-"),i]   == x[(n,"-"),i])

    return m
end

p=MP!(ms,mp,unc,ps,pp)
optimize!(p)

## output results to directory "plots"
# i=0;investment_cost=[];operational_cost=[];emission=[];co2tax=[];energy_loss=[];location=[];tech=[];line=[];electricity_shed=[];heat_shed=[];co2tax2=[];battery_cap=[];
# wind_cap=[];solar_cap=[];electrolyser_cap=[];fuel_cell_cap=[];hydrogen_store_cap=[];converter_cap=[];co2tax3=[];line_cap=[];line=[];location_EH=[];E_boiler_cap=[];
# co2tax4=[];power_demand=[];heat_demand=[];wind_power=[];gas_power=[];gas_heat=[];onshore_power=[];fuel_cell_power=[];co2tax5=[]; location_ZO=[]; onshore_cap=[];co2budget=[];
# co2budget2=[];co2budget3=[];co2budget4=[];co2budget5=[];el_battery=[];el_boiler=[];el_separator=[];el_inj_comp=[];el_pump_oil=[];el_comp=[];el_wi=[];el_wl=[];el_el=[];el_fc=[];
# el_gs1=[];el_gs2=[];el_gs3=[];el_gs4=[];el_gs5=[];el_gs6=[];el_gs7=[];el_gs8=[];el_gs9=[];el_gs10=[];el_gs11=[];el_gs12=[];el_gs13=[];el_gs14=[];el_gs15=[];
# el_hgs1=[];el_hgs2=[];el_hgs3=[];el_hgs4=[];el_hgs5=[];el_OCGT=[];el_l=[];
# while i<=250
#     i=i+1
#     setfield!(unc,:c,[20.0*(i-1) 20.0*(i-1)])
#     p=MP!(ms,mp,unc,ps,pp)
#     optimize!(p)
#
#         push!(co2tax,20.0*(i-1))
#         push!(investment_cost,value.(p[:f]))
#         push!(operational_cost, value.(p[:β][1]))
#         push!(emission, value.( p[:ϕ]["co2_price",1]))
#         push!(energy_loss, ((sum((1/pp.ηG-1-pp.ηGh)value.(p[:pG]["OCGT",z,h,1])
#         +(1-pp.ηSE)value.(p[:pSEI]["battery",z,h,1])
#         +(1-pp.ηEB)value.(p[:pEB]["E_boiler",z,h,1])
#         +(1-0.9)value.(p[:pSEPH]["separator",z,h,1])
#         +(1-0.8)value.(p[:pGI]["gas_inj_compressor",z,h,1]+p[:pC]["compressor",z,h,1])
#         +(1-0.65)value.(p[:pO]["oil_exp_pump",z,h,1])
#         +(1-0.75)value.(p[:pWI]["water_inj_pump",z,h,1]+p[:pWL]["water_lift_pump",z,h,1])
#         +value.(p[:pHGShed][z,h,1]) for z in ps.ZP for h in ps.sH)
#         +sum(value.(p[:pE]["electrolyser",z,h,1]-pp.θHY*(p[:vEF]["electrolyser",z,h,1]+p[:vSHyI]["hydrogen_store",z,h,1])
#         + (1/pp.ηF-1)*p[:pF]["fuel_cell",z,h,1]) for z in ps.ZH for h in ps.sH))*pp.WS
#         + sum(value.(p[:pGShed][z,h,1] for z in ps.Z for h in ps.sH))*pp.WS
#         + sum(abs((1-pp.ηL[pp.map_l[l]])*value.(p[:pL][l,h,1])) for l in ps.L for h in ps.sH)*pp.WS)/1000) # in GWh
#
#         push!(el_OCGT, (sum((1/pp.ηG-1-pp.ηGh)value.(p[:pG]["OCGT",z,h,1]) for z in ps.ZP for h in ps.sH)*pp.WS)/1000)
#         push!(el_battery, (sum((1-pp.ηSE)value.(p[:pSEI]["battery",z,h,1]) for z in ps.ZP for h in ps.sH)*pp.WS)/1000)
#         push!(el_boiler, (sum((1-pp.ηEB)value.(p[:pEB]["E_boiler",z,h,1]) for z in ps.ZP  for h in ps.sH)*pp.WS)/1000)
#         push!(el_separator, (sum((1-0.9)value.(p[:pSEPH]["separator",z,h,1]) for z in ps.ZP for h in ps.sH)*pp.WS)/1000)
#         push!(el_inj_comp, (sum((1-0.8)value.(p[:pGI]["gas_inj_compressor",z,h,1]) for z in ps.ZP for h in ps.sH)*pp.WS)/1000)
#         push!(el_comp, (sum((1-0.8)value.(p[:pC]["compressor",z,h,1]) for z in ps.ZP for h in ps.sH)*pp.WS)/1000)
#         push!(el_pump_oil, (sum((1-0.65)value.(p[:pO]["oil_exp_pump",z,h,1]) for z in ps.ZP for h in ps.sH)*pp.WS)/1000)
#         push!(el_wi, (sum((1-0.75)value.(p[:pWI]["water_inj_pump",z,h,1]) for z in ps.ZP for h in ps.sH)*pp.WS)/1000)
#         push!(el_wl, (sum((1-0.75)value.(p[:pWL]["water_lift_pump",z,h,1]) for z in ps.ZP for h in ps.sH)*pp.WS)/1000)
#         push!(el_hgs1, sum(value.(p[:pHGShed]["cluster1",h,1]) for h in ps.sH)*pp.WS/1000)
#         push!(el_hgs2, sum(value.(p[:pHGShed]["cluster2",h,1]) for h in ps.sH)*pp.WS/1000)
#         push!(el_hgs3, sum(value.(p[:pHGShed]["cluster3",h,1]) for h in ps.sH)*pp.WS/1000)
#         push!(el_hgs4, sum(value.(p[:pHGShed]["cluster4",h,1]) for h in ps.sH)*pp.WS/1000)
#         push!(el_hgs5, sum(value.(p[:pHGShed]["cluster5",h,1]) for h in ps.sH)*pp.WS/1000)
#         push!(el_gs1, sum(value.(p[:pGShed]["cluster1",h,1]) for h in ps.sH)*pp.WS/1000)
#         push!(el_gs2, sum(value.(p[:pGShed]["cluster2",h,1]) for h in ps.sH)*pp.WS/1000)
#         push!(el_gs3, sum(value.(p[:pGShed]["cluster3",h,1]) for h in ps.sH)*pp.WS/1000)
#         push!(el_gs4, sum(value.(p[:pGShed]["cluster4",h,1]) for h in ps.sH)*pp.WS/1000)
#         push!(el_gs5, sum(value.(p[:pGShed]["cluster5",h,1]) for h in ps.sH)*pp.WS/1000)
#         push!(el_gs6, sum(value.(p[:pGShed]["EH1",h,1]) for h in ps.sH)*pp.WS/1000)
#         push!(el_gs7, sum(value.(p[:pGShed]["EH2",h,1]) for h in ps.sH)*pp.WS/1000)
#         push!(el_gs8, sum(value.(p[:pGShed]["EH3",h,1]) for h in ps.sH)*pp.WS/1000)
#         push!(el_gs9, sum(value.(p[:pGShed]["EH4",h,1]) for h in ps.sH)*pp.WS/1000)
#         push!(el_gs10, sum(value.(p[:pGShed]["EH5",h,1]) for h in ps.sH)*pp.WS/1000)
#         push!(el_gs11, sum(value.(p[:pGShed]["onshore1",h,1]) for h in ps.sH)*pp.WS/1000)
#         push!(el_gs12, sum(value.(p[:pGShed]["onshore2",h,1]) for h in ps.sH)*pp.WS/1000)
#         push!(el_gs13, sum(value.(p[:pGShed]["onshore3",h,1]) for h in ps.sH)*pp.WS/1000)
#         push!(el_gs14, sum(value.(p[:pGShed]["onshore4",h,1]) for h in ps.sH)*pp.WS/1000)
#         push!(el_gs15, sum(value.(p[:pGShed]["onshore5",h,1]) for h in ps.sH)*pp.WS/1000)
#         push!(el_el,sum(value.(p[:pE]["electrolyser",z,h,1]-pp.θHY*(p[:vEF]["electrolyser",z,h,1]+p[:vSHyI]["hydrogen_store",z,h,1])) for h in ps.sH for z in ps.ZH)*pp.WS/1000)
#         push!(el_fc, sum((1/pp.ηF-1)*value.(p[:pF]["fuel_cell",z,h,1]) for z in ps.ZH for h in ps.sH)*pp.WS/1000)
#         push!(el_l, sum(abs((1-pp.ηL[pp.map_l[l]])*value.(p[:pL][l,h,1])) for l in ps.L for h in ps.sH)*pp.WS/1000)
#
#
#         for z in ms.ZP
#             push!(co2tax2,20.0*(i-1))
#             push!(location,z)
#             push!(electricity_shed,sum(value.(p[:pLShed][z,h,1] for h in ps.sH)))
#             push!(heat_shed,sum(value.(p[:pHLShed][z,h,1] for h in ps.sH))/2880)
#             push!(power_demand, sum(value.(p[:pD][z,h,1] for h in ps.sH))/2880)
#             push!(heat_demand, sum(value.(p[:pSEPH][dsep,z,h,1] for dsep in ps.Dsep for h in ps.sH))/2880)
#             push!(gas_power, sum(value.(p[:pG][g,z,h,1] for g in ps.G for h in ps.sH))/2880)
#             push!(gas_heat, sum(value.(pp.ηGh*p[:pG][g,z,h,1] for g in ps.G for h in ps.sH))/2880)
#             push!(battery_cap,value.(p[:x0pp]["battery",z,1])*exp10(3))
#             push!(E_boiler_cap,value.(p[:x0pp]["E_boiler",z,1])*exp10(3))
#
#         end
#
#         for z in ms.ZH
#             push!(co2tax4,20.0*(i-1))
#             push!(location_EH,z)
#             push!(wind_cap,value.(p[:x0ph]["wind",z,1])*exp10(3))
#             push!(solar_cap,value.(p[:x0ph]["solar",z,1])*exp10(3))
#             push!(electrolyser_cap,value.(p[:x0ph]["electrolyser",z,1])*exp10(3))
#             push!(fuel_cell_cap,value.(p[:x0ph]["fuel_cell",z,1])*exp10(3))
#             push!(hydrogen_store_cap,value.(p[:x0ph]["hydrogen_store",z,1])) # in tonne
#             push!(converter_cap,value.(p[:x0ph]["converter",z,1])*exp10(3))
#             push!(wind_power, sum(value.(p[:pRW][w,z,h,1] for w in ps.RW for h in ps.sH))/2880)
#             push!(fuel_cell_power, sum(value.(p[:pF][f,z,h,1] for f in ps.F for h in ps.sH))/2880)
#
#         end
#
#         for z in ms.ZO
#             push!(co2tax5,20.0*(i-1))
#             push!(location_ZO,z)
#             push!(onshore_power,sum(value.(p[:pZO][b,z,h,1] for b in ps.B for h in ps.sH))/2880)
#             push!(onshore_cap, maximum(value.(p[:pZO][ps.B,z,ps.sH,1])))
#         end
#
#         for l in ps.L
#             push!(co2tax3,20.0*(i-1))
#             push!(line_cap,value.(p[:x0l][l,1])*exp10(3))
#             push!(line,l)
#         end
# end
#
# cost_D=DataFrame(co2tax=co2tax,investment_cost=investment_cost, operational_cost=operational_cost)
# emission_D=DataFrame(co2tax=co2tax, emission=emission, energy_loss=energy_loss, el_battery=el_battery, el_boiler=el_boiler, el_separator=el_separator, el_inj_comp=el_inj_comp,
# el_pump_oil=el_pump_oil,el_comp=el_comp, el_wi=el_wi,el_wl=el_wl,el_el=el_el,el_fc=el_fc,
# el_gs1=el_gs1,el_gs2=el_gs2,el_gs3=el_gs3,el_gs4=el_gs4,el_gs5=el_gs5,el_gs6=el_gs6,el_gs7=el_gs7,el_gs8=el_gs8,el_gs9=el_gs9,el_gs10=el_gs10,el_gs11=el_gs11,el_gs12=el_gs12,
# el_gs13=el_gs13,el_gs14=el_gs14,el_gs15=el_gs15,
# el_hgs1=el_hgs1,el_hgs2=el_hgs2,el_hgs3=el_hgs3,el_hgs4=el_hgs4,el_hgs5=el_hgs5,el_OCGT=el_OCGT,el_l=el_l)
# shed_D=DataFrame(co2tax=co2tax2,location=location,electricity_shed=electricity_shed,heat_shed=heat_shed)
# tech_D=DataFrame(co2tax=co2tax2,location=location_EH,wind=wind_cap,solar=solar_cap,electrolyser=electrolyser_cap,fuel_cell=fuel_cell_cap,hydrogen=hydrogen_store_cap,converter=converter_cap,battery=battery_cap,E_boiler=E_boiler_cap)
# line_D=hcat(DataFrame(line=line),DataFrame(co2tax=co2tax3),DataFrame(line_cap=line_cap))
# power_demand_D=DataFrame(co2tax=co2tax2,power_demand=power_demand,location=location)
# heat_demand_D=DataFrame(co2tax=co2tax2,heat_demand=heat_demand,location=location)
# gas_power_D=DataFrame(co2tax=co2tax2,OCGT=gas_power,location=location)
# power_supply_D=DataFrame(co2tax=co2tax4,wind_power=wind_power,fuel_cell_power=fuel_cell_power,location=location_EH)
# power_onshore_D=DataFrame(co2tax=co2tax5,onshore_power=onshore_power,location=location_ZO,onshore_cap=onshore_cap)
#
# XLSX.writetable(pwd()*"/plots/co2tax_sensitivity_5000.xlsx",overwrite=true,
# cost=(collect(DataFrames.eachcol(cost_D)), DataFrames.names(cost_D)),
# emission=(collect(DataFrames.eachcol(emission_D)), DataFrames.names(emission_D)),
# shed=(collect(DataFrames.eachcol(shed_D)), DataFrames.names(shed_D)),
# tech=(collect(DataFrames.eachcol(tech_D)), DataFrames.names(tech_D)),
# line=(collect(DataFrames.eachcol(line_D)), DataFrames.names(line_D)),
# power_demand=(collect(DataFrames.eachcol(power_demand_D)), DataFrames.names(power_demand_D)),
# heat_demand=(collect(DataFrames.eachcol(heat_demand_D)), DataFrames.names(heat_demand_D)),
# gas_power=(collect(DataFrames.eachcol(gas_power_D)), DataFrames.names(gas_power_D)),
# power_supply=(collect(DataFrames.eachcol(power_supply_D)), DataFrames.names(power_supply_D)),
# power_onshore=(collect(DataFrames.eachcol(power_onshore_D)), DataFrames.names(power_onshore_D)));


text:
"""
    CSVtoDataFrame(directory::String)

Read all the Boletin Estadístico files in a directory, tranform them into a DataFrame.
Return a Dictionary with pairs (filename, dataframe)

"""
function CSVtoDataFrame(directory::String,from::Date=Date(1900))

    #Definición de formato de fechas utilizado por Banco de España en su boletín estadístico
    bde_months = ["ENERO", "FEBRERO", "MARZO", "ABIRL", "MAYO", "JUNIO", "JULIO", "AGOSTO", "SEPTIEMBRE", "OCTUBRE", "NOVIEMBRE", "DICIEMBRE"]
    bde_monts_abbrev = ["ENE","FEB","MAR","ABR","MAY","JUN", "JUL","AGO","SEP","OCT","NOV","DIC"]

    Dates.LOCALES["BdE"] = Dates.DateLocale(bde_months,bde_monts_abbrev, [""],[""],)

    BdE_MonthYear = DateFormat("uuu yyyy", "BdE") #Definición del formato de fechas utilizado por BdE
    BdE_DayMonthYear = DateFormat("dd uuu yyyy", "BdE")

    File_Dict = Dict{String,DataFrame}() #Inicializo el Diccionario con los datos

    #Lectura datos ficheros BdE. La carpeta debe estar en una carpeta llamada "be04"
    for file in readdir(directory)

        #Lectura fichero transformando Latin1 a UTF8
        tempdf = CSV.File(open(read,joinpath(directory,file), enc"Latin1"), header=4, datarow=7, footerskip=2, delim=',',
                        decimal='.',quotechar='"',missingstrings=["_","..."]) |> DataFrame

        #Conversión de datos (distinguiendo cuando hay año, mes y año, día, mes y año)
        if typeof(tempdf[1,1])== Int
           tempdf[!,1] = Date.(tempdf[:,1],12,31)

        elseif isdigit(tempdf[1,1][1])
            tempdf[!,1] = Date.(tempdf[:,1],BdE_DayMonthYear)

        elseif typeof(tempdf[1,1])==String
            tempdf[!,1] = lastdayofmonth.(Date.(tempdf[:,1],BdE_MonthYear))

        end

        #Filtrado por fecha según from
        tempdf = tempdf[(tempdf[:,1]).>=from,:]

        #Creación de Dict para guardar cada dataframe asociado al nombre del fichero según Banco de España
        filename = file[1:length(file)-4]
        push!(File_Dict,filename => tempdf)

    end #for

    File_Dict

end
"""
Transpose a Dataframe
"""

function transposedf(df::DataFrame; columna1::Symbol = :Concepto)

    DataFrame([[names(df[:,2:end])]; collect.(eachrow(df[:,2:end]))], [columna1;Symbol.(df[:,1])])

end #function

"""

Graba un dataframe con codificación latin1 para capturar correctamente los caracteres de un fichero de Banco de España

"""

function DataFrametoCSV(df::DataFrame, directory::String, filename::String)

    UTF8_IO = IOBuffer()

    CSV.write(UTF8_IO,df, delim=';',decimal=',',quotestrings=true)

    file = open(joinpath(directory,filename),"w")
    encoder = StringEncoder(file, enc"Latin1")

    write(encoder,read(seekstart(UTF8_IO),String))
    close(encoder) #codifica el stream de datos en Latin1
    close(file) #graba el stream de datos

end


text:
using Null_models_I_and_II
using Test

@testset "Null_models_I_and_II.jl" begin
    # Write your tests here.
end


text:
# ---
# title: 1619. Mean of Array After Removing Some Elements
# id: problem1619
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/mean-of-array-after-removing-some-elements/description/>
# hidden: true
# ---
#
# Given an integer array `arr`, return _the mean of the remaining integers after
# removing the smallest`5%` and the largest `5%` of the elements._
#
# Answers within `10-5` of the **actual answer** will be considered accepted.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3]
#     Output: 2.00000
#     Explanation: After erasing the minimum and the maximum values of this array, all elements are equal to 2, so the mean is 2.
#
#
# **Example 2:**
#
#
#
#     Input: arr = [6,2,7,5,1,2,0,3,10,2,5,0,5,5,0,8,7,6,8,0]
#     Output: 4.00000
#
#
# **Example 3:**
#
#
#
#     Input: arr = [6,0,7,0,7,5,7,8,3,4,0,7,8,1,6,8,1,1,2,4,8,1,9,5,4,3,8,5,10,8,6,6,1,0,6,10,8,2,3,4]
#     Output: 4.77778
#
#
# **Example 4:**
#
#
#
#     Input: arr = [9,7,8,7,7,8,4,4,6,8,8,7,6,8,8,9,2,6,0,0,1,10,8,6,3,3,5,1,10,9,0,7,10,0,10,4,1,10,6,9,3,6,0,0,2,7,0,6,7,2,9,7,7,3,0,1,6,1,10,3]
#     Output: 5.27778
#
#
# **Example 5:**
#
#
#
#     Input: arr = [4,8,4,10,0,7,1,3,7,8,8,3,4,1,6,2,1,1,8,0,9,8,0,3,9,10,3,10,1,10,7,3,2,1,4,9,10,7,6,4,0,8,5,1,2,1,6,2,5,0,7,10,9,10,3,7,10,5,8,5,7,6,7,6,10,9,5,10,5,5,7,2,10,7,7,8,2,0,1,1]
#     Output: 5.29167
#
#
#
#
# **Constraints:**
#
#   * `20 <= arr.length <= 1000`
#   * `arr.length` ******is a multiple** of `20`.
#   * `0 <= arr[i] <= 105`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@doc """
    CenteredDifference{N}(n, order, step, len, [coeff_func])
    CenteredDifference{N}(n, order, steps, len, [coeff_func])

See also: [`UpwindDifference`](@ref)
"""
CenteredDifference

@doc """
    calculate_weights(n::Int, x₀::Real, x::Vector)

Return a vector `c` such that `c⋅f.(x)` approximates ``f^{(n)}(x₀)`` for smooth `f`.

The points `x` need not be evenly spaced.

The stencil `c` has the highest approximation order possible given values of `f` at `length(x)` points.  More precisely, if `x` has length `m`, there is a function `g` such that ``g(y) = f(y) + O(y-x₀)^{m-n+?}`` and ``c⋅f.(x) = g'(x₀)``.

The algorithm is due to [Fornberg](https://doi.org/10.1090/S0025-5718-1988-0935077-0), with a [modification](http://epubs.siam.org/doi/pdf/10.1137/S0036144596322507) to improve stability.
"""
calculate_weights

@doc """
    DerivativeOperator{T<:Real,N,Wind,T2,S1,S2<:SVector,T3,F}(<fields>)

Represent a finite difference derivative operator.

These operators implement the `DiffEqBase.AbstractDiffEqLinearOperator` interface.  Therefore `eltype` returns `T`.

These operators can be contracted over an arbitrary dimension, given by the type parameter `N`.

The `CentredDifference` and `UpwindDifference` types serve as hooks to hang constructors, but no objects are actually constructed with those types.  Their constructors return structures of type `DerivativeOperator`, with stencils constructed by `calculate_weights`, and other fields filled out appropriately.

The finite-difference methods are defined for `DerivativeOperator`.  In particular, `*`, which actually takes derivatives.

The key data are three stencils of coefficients, stored in `SVector`s.  The interior stencil, stored in `stencil_coefs::S1`, is the normal one used in the interior of the grid.  The others, `low_boundary_coefs::S2` and `high_boundary_coefs::S2`, are used where the normal stencil would jut out of the grid boundary.  These can have a different length than the interior stencil, hence the two types.

When the operator is applied by `mul!`, these stencils are multiplied by the vector in the `coefficients` field.  It is set to `coeff_func.(1:len)` in `UpwindDifference`, which seems odd, because `len` is not the stencil length.  Scalar multiplication is absorbed into this vector, as required by the `DiffEqBase.AbstractDiffEqLinearOperator` interface.

The `coefficients` field appears to be a more general `DifferentialEquations` thing.  There is an `update_coefficients!` method.

The term “left boundary” is used interchangeably with “low boundary”, and “right” with “high”.

# Type Parameters

Presumably, most of these are to force method specialization on:a) the stencil lengths, which are given by the parameters of the `SVector` types `S1` and `S2`, and b) even or arbitrary grid spacing.

- `T<:Real`: Function range type.
- `N`: Contraction dimension
- `Wind`: Flag `true` for upwind operators, `false` for centered.
- `T2`: Concrete type of `dx` field, which might be `T` or `AbstractVector{T}`.
- `S1`: Concrete type of the interior stencil.  (When is this not an `SVector`?)
- `S2<:SVector`: Concrete type of boundary coefficients.
- `T3`: `Vector{T}` if a `coeff_func` was supplied, `Nothing` otherwise.
- `F`: `typeof(coeff_func)` if a `coeff_func` was supplied, `Nothing` otherwise.

# Fields

- `derivative_order::Int`: The `n` in `f⁽ⁿ⁾(x)`.
- `approximation_order::Int`: The degree of polynomial for which the operator is numerically exact.
- `dx::T2`: Grid step.
- `len::Int`: The number of interior points on the grid.
- `stencil_length::Int`: Length of the interior stencil
- `stencil_coefs::S1`: An `SVector` of the interior stencil coefficients.
- `boundary_stencil_length::Int`:
- `boundary_point_count::Int`:
- `low_boundary_coefs::S2`: An array of `SVector` stencils for the interior points where `stencil_coefs` juts over the left boundary.
- `high_boundary_coefs::S2`: An array of `SVector` stencils for the interior points where `stencil_coefs` juts over the right boundary.
- `coefficients::T3`: Multiplier for the stencil.
- `coeff_func::F`:  God only knows what this function is for.

See also: [`AbstractBC`](@ref), [`CenteredDifference`](@ref), [`UpwindDifference`](@ref)
"""
DerivativeOperator

@doc """
    abstract type AbstractBC{T}

`T` is the range type of the discretised function.

Boundary condition operators extrapolate the discretised function
as illustrated in the README, adding a ghost node at each end such
that an interpolated polynomial satisfies the boundary condition.

Some of these are affine operators with a constant term, not linear
operators, so using `*` and `\\` is a minor abuse of notation.  For example,
when there is a non-zero Dirichlet boundary value, the value at the
ghost node is constant, and does not scale proportional to the
interior function values.  Because boundary condition operators are
not linear, they can not be concretised as matrices.

`GhostDerivativeOperator` somehow concretises as a pair of matrices.

The implementation of `*(::AffineBC, ::Vector)` returns a
`BoundaryPaddedVector`.  This is a structure that stores the ghost
values separately from the original vector.  Presumably this is to
avoid allocations, although the `mul!` method for `GhostDerivativeOperator`
and `Matrix` method allocates a new `BoundaryPaddedArray` for each
column.

In higher dimensional cases, it might return a
`BoundaryPaddedArray`.

The `*` operator for a derivative multiplying a boundary condition
packages them up as a `GhostDerivativeOperator`.  This can not
return a sparse matrix, because the operator is affine not linear.
When a vector is multiplied by one of those, it uses `*(::AbstractBC,
::AbstractVector)` to pad with ghost nodes, then
`mul!(::DerivativeOperator, ::AbstractVector)` to evaluate the
derivative.

See also: [`DerivativeOperator`](@ref), [`PeriodicBC`](@ref) [`NeumannBC`](@ref), [`DirichletBC`](@ref), [`RobinBC`](@ref)
"""
AbstractBC

@doc """
    PeriodicBC{T}()

Only one ghost node is added at each end.  This means that periodic
functions will not be consistent at the interior and near the
boundary with high order stencils.  It is likely to cause other
limitations too.

Currently, periodic boundary conditions are only implemented in one
dimension.

# Type parameters

- `T` is the domain type of the discretised function.

See also: [`AbstractBC`](@ref)
"""
PeriodicBC

@doc """
NeumannBC(α::NTuple{2,T}, dx::Union{AbstractVector{T}, T}, order = 1) where T = RobinBC((zero(T), one(T), α[1]), (zero(T), one(T), α[2]), dx, order)
DirichletBC(αl::T, αr::T) where T = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), 1.0, 2.0 )
#specialized constructors for Neumann0 and Dirichlet0
Dirichlet0BC(T::Type) = DirichletBC(zero(T), zero(T))
Neumann0BC(dx::Union{AbstractVector{T}, T}, order = 1) where T = NeumannBC((zero(T), zero(T)), dx, order)

See also: [`AbstractBC`](@ref)
"""
RobinBC
@doc (@doc RobinBC) NeumannBC
@doc (@doc RobinBC) Neumann0BC
@doc (@doc RobinBC) DirichletBC
@doc (@doc RobinBC) Dirichlet0BC


text:
"""
$(TYPEDEF)

A scalar equation for optimization.

# Fields
$(FIELDS)

# Examples

```
@variables x y z
@parameters σ ρ β

op = σ*(y-x) + x*(ρ-z)-y + x*y - β*z
os = OptimizationSystem(eqs, [x,y,z],[σ,ρ,β])
```
"""
struct OptimizationSystem <: AbstractSystem
    """Vector of equations defining the system."""
    op::Operation
    """Unknown variables."""
    states::Vector{Variable}
    """Parameters."""
    ps::Vector{Variable}
    """
    Name: the name of the system
    """
    name::Symbol
    """
    systems: The internal systems
    """
    systems::Vector{OptimizationSystem}
end

function OptimizationSystem(op, states, ps;
                            name = gensym(:OptimizationSystem),
                            systems = OptimizationSystem[])
    OptimizationSystem(op, convert.(Variable,states), convert.(Variable,ps), name, systems)
end

function calculate_gradient(sys::OptimizationSystem)
    expand_derivatives.(gradient(equations(sys), [dv() for dv in states(sys)]))
end

function generate_gradient(sys::OptimizationSystem, vs = states(sys), ps = parameters(sys); kwargs...)
    grad = calculate_gradient(sys)
    return build_function(grad, convert.(Variable,vs), convert.(Variable,ps);
                          conv = AbstractSysToExpr(sys),kwargs...)
end

function calculate_hessian(sys::OptimizationSystem)
    expand_derivatives.(hessian(equations(sys), [dv() for dv in states(sys)]))
end

function generate_hessian(sys::OptimizationSystem, vs = states(sys), ps = parameters(sys);
                          sparse = false, kwargs...)
    hes = calculate_hessian(sys)
    if sparse
        hes = sparse(hes)
    end
    return build_function(hes, convert.(Variable,vs), convert.(Variable,ps);
                          conv = AbstractSysToExpr(sys),kwargs...)
end

function generate_function(sys::OptimizationSystem, vs = states(sys), ps = parameters(sys); kwargs...)
    vs′ = convert.(Variable,vs)
    ps′ = convert.(Variable,ps)
    return build_function(equations(sys), vs′, ps′;
                          conv = AbstractSysToExpr(sys),kwargs...)
end

equations(sys::OptimizationSystem) = isempty(sys.systems) ? sys.op : sys.op + reduce(+,namespace_operation.(sys.systems))
namespace_operation(sys::OptimizationSystem) = namespace_operation(sys.op,sys.name,nothing)

"""
```julia
function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem,
                                          parammap=DiffEqBase.NullParameters();
                                          u0=nothing, lb=nothing, ub=nothing,
                                          hes = false, sparse = false,
                                          checkbounds = false,
                                          linenumbers = true, parallel=SerialForm(),
                                          kwargs...) where iip
```

Generates an OptimizationProblem from an OptimizationSystem and allows for automatically
symbolically calculating numerical enhancements.
"""
function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem,
                                          parammap=DiffEqBase.NullParameters();
                                          u0=nothing, lb=nothing, ub=nothing,
                                          hes = false, sparse = false,
                                          checkbounds = false,
                                          linenumbers = true, parallel=SerialForm(),
                                          kwargs...) where iip
    dvs = states(sys)
    ps = parameters(sys)

    f = generate_function(sys,checkbounds=checkbounds,linenumbers=linenumbers,
                              parallel=parallel,expression=Val{false})
    u0 = varmap_to_vars(u0,dvs)
    p = varmap_to_vars(parammap,ps)
    lb = varmap_to_vars(lb,dvs)
    ub = varmap_to_vars(ub,dvs)
    OptimizationProblem(f,p;u0=u0,lb=lb,ub=ub,kwargs...)
end

"""
```julia
function DiffEqBase.OptimizationProblemExpr{iip}(sys::OptimizationSystem,
                                          parammap=DiffEqBase.NullParameters();
                                          u0=nothing, lb=nothing, ub=nothing,
                                          hes = false, sparse = false,
                                          checkbounds = false,
                                          linenumbers = true, parallel=SerialForm(),
                                          kwargs...) where iip
```

Generates a Julia expression for an OptimizationProblem from an
OptimizationSystem and allows for automatically symbolically
calculating numerical enhancements.
"""
struct OptimizationProblemExpr{iip} end

function OptimizationProblemExpr{iip}(sys::OptimizationSystem,
                                          parammap=DiffEqBase.NullParameters();
                                          u0=nothing, lb=nothing, ub=nothing,
                                          hes = false, sparse = false,
                                          checkbounds = false,
                                          linenumbers = false, parallel=SerialForm(),
                                          kwargs...) where iip
    dvs = states(sys)
    ps = parameters(sys)

    f = generate_function(sys,checkbounds=checkbounds,linenumbers=linenumbers,
                              parallel=parallel,expression=Val{true})
    u0 = varmap_to_vars(u0,dvs)
    p = varmap_to_vars(parammap,ps)
    lb = varmap_to_vars(lb,dvs)
    ub = varmap_to_vars(ub,dvs)
    quote
        f = $f
        p = $p
        u0 = $u0
        lb = $lb
        ub = $ub
        OptimizationProblem(f,p;u0=u0,lb=lb,ub=ub,kwargs...)
    end
end