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

using Markdown
using InteractiveUtils

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

# ╔═╡ 88b46d2e-220e-11eb-0f7f-b3f523f0214e
begin
	using Plots
	using PlutoUI
	using LaTeXStrings
	using Roots
end

# ╔═╡ 6f554dc0-220c-11eb-341f-1b66de841c86
md"""
## Nonlinear dynamics: stability and bifurcations
"""

# ╔═╡ 816de40c-220c-11eb-2d3a-23e6267cd529
md"""
In this module on climate modelling we will be studying various **differential equations** that model the evolution of aspects of the climate over time.

The simplest such models are **ordinary differential equations** (ODEs), where one or a few continuous variables evolve continuously in time, with a model that specifies their instantaneous rate of change as a function of their current values.

The general form of an ODE (that doesn't depend explicitly on time, i.e. that is **autonomous**), is

$$\frac{dx(t)}{dt} = f(x(t)),$$

or

$$\dot{x}(t) = f(x(t)),$$

with an initial condition $x(t=0) = x_0$.

Here, $\dot{x}(t)$ denotes the derivative of the function $t \mapsto x(t)$ at time $t$.
"""

# ╔═╡ d21f6358-220c-11eb-00b8-03fe0746fdc9
md"""
The simplest numerical method to solve such an equation is the **(forward) Euler method**, in which we convert this equation into an explicit time-stepping routine: we take a small time step of length $h$ and approximate the derivative as

$$\frac{dx(t)}{dt} \simeq \frac{x(t + h) - x(t)}{h},$$

giving the approximation

$$x(t + h) \simeq x(t) + h \, f(x(t))$$

for the approximation at the next time step.
"""

# ╔═╡ 735bb47e-220d-11eb-1ba6-8d591a935857
md"""
## Solving the ODE: Euler method
"""

# ╔═╡ 79b148b6-220d-11eb-2785-25d05068aeed
md"""
Let's use this to simulate a simple nonlinear ODE that describes the dynamics of a population of bacteria. The bacteria will grow by reproduction at a rate $\lambda$ provided there is sufficient food, in which case we would have $\dot{x} = \lambda x$. But the available food will actually always limit the sustainable population to a value $K$. A simple model for this is as follows:

$$\dot{x} = \lambda \, x \, (K - x).$$

When $x$ is close to $0$, the growth rate is $\lambda$, but that rate decreases as $x$ increases.

This is sometimes called the [**logistic** differential equation](https://en.wikipedia.org/wiki/Logistic_function#Logistic_differential_equation) (although the name does not seem particularly helpful).
"""

# ╔═╡ d43a869e-220d-11eb-009c-a119de57e7da
md"""
Our goal is to use computational thinking, but we will actually not be interested so much in the exact dynamics in time, but rather in the **qualitative** features of the behaviour of the system. For example, at long times (formally $t \to \infty$) does the population get arbitrarily large? Or does it, for example, oscillate around a particular value? Or does it converge to a particular size? This forms the subject of **nonlinear dynamics** or **dynamical systems** theory.
"""

# ╔═╡ 088eabb2-220e-11eb-1ac0-df419b87f39a
md"""
Let's simulate the system using the Euler method to try to guess the answer to this question. Note that there are many much more sophisticated methods for solving ODEs collected in the DifferentialEquations.jl package. We should never use the Euler method in practice, but should use a tried and tested library instead, and algorithms that provide much better accuracy in the solutions, if we are interested in faithful numerical results.
"""

# ╔═╡ a44cc180-22d4-11eb-2129-211d024c4921
md"""
We'll rescale the variables to the simplest form:

$$\dot{x} = x \, (1 - x).$$
"""

# ╔═╡ 6aa1d2f4-220e-11eb-06e0-c346f74aa018
logistic(x) = x * (1 - x)

# ╔═╡ b6d0dd18-22d4-11eb-2083-0b02de030579
md"""
Let's simulate this with Euler and plot the trajectory $x(t)$ as a function of time $t$:
"""

# ╔═╡ f055be76-22d4-11eb-268a-bf70c7d8f1a1
md"""
We see that for this particular initial condition, the solution seems to settle down to a fixed value after some time, and then remains at that value thereafter.

Such a value is called a **fixed point** or a **stationary point** of the ODE.
"""

# ╔═╡ 90ccb392-2216-11eb-1fd8-83b7d7c16b54
md"""
## Qualitative behaviour: Fixed points and their stability
"""

# ╔═╡ bef9c2e4-220e-11eb-24d8-bd618d2985ea
md"""
 Let's see if what happens for other initial conditions:
"""

# ╔═╡ 17317eb8-22d5-11eb-0f46-2bf1bddd36bb
@bind x0 Slider(-1.0:0.001:20, default=0.0, show_value=true)

# ╔═╡ ae862e10-22d5-11eb-3d75-1d748cf86944
md"""
To get an overview of the behaviour we can also draw all the results on a single graph:
"""

# ╔═╡ 446c564e-220f-11eb-191a-6b419e790f3f
md"""
We see that all the curves starting near to $x_0=1.0$ seem to converge to 1 at long times. If the system starts *exactly* at 0 then it stays there forever. However, if it starts close to 0, on either side, then it moves *away* from 0 (on that same side of 0) -- starting from a negative value $x$ becomes ever more negative. (Even though negative populations have no meaning in the original interpretation as the dynamics of a population, we can still ask study the dynamics of the equation with negative initial conditions, since it may model other systems too.)

The special values $x^*_1=1$ and $x^*_2=0$ are called **stationary points** or **fixed points** of the differential equation. If we start at $x^*_i$, then the derivative there is $f'(x^*_i) = 0$, and hence we cannot move away from $x^*_i$! The fixed points can be found as zeros or **roots** of the function $f$, i.e. values $x^*$ such that $f(x^*) = 0$.
"""

# ╔═╡ 38894c32-2210-11eb-26ff-d9796edd7871
md"""
We see, though, that the two types of fixed points are **qualitatively different**: trajectories that start close to $x^*_1 = 1$ move *towards* $x^*_1$, whereas trajectories that start close to $x^*_2 = 0$ move *away* from it. We say that $x^*_1$ is a **stable fixed point** and $x^*_2$ is an **unstable fixed point**.

In general it is not possible to find analytical formulas for the position and stability of fixed points; instead, we can use numerical **root-finding algorithms**, for example the Newton method: see e.g. course 18.330 next semester. There are various Julia packages for root finding; we will use `Roots.jl`.
"""

# ╔═╡ a196a20a-2216-11eb-11f9-6789bb2a3ece
md"""
## State space: Vector field and phase portrait
"""

# ╔═╡ c15a53e8-2216-11eb-1460-8510aea0f570
md"""
If we want to find the whole trajectory for a given initial condition then we need to solve the equations, either numerically or analytically.

However, we may want less information about the system, for example the **long-time** or **asymptotic** dynamics. It turns out that we can obtain some information about that *without* explicitly solving the ODE! This is the **qualitative approach** to studying nonlinear systems.

Instead of drawing trajectories $x(t)$ as a function of time $t$, as we did above, let's use a different graphical representation, where we draw **state space** or **phase space**: This is the set ("space") of all possible values of the dependent variables ("states"). For the above ODE there is only a single dependent variable, $x$, so the state space is the real line, $\mathbb{R}$.

At each possible value of $x$, the ODE gives us information about the rate of change of $x(t)$ at that point. Let's draw an **arrow** at that point, pointing in the direction that a particle placed at that point would move: to the right if $\dot{x} > 0$ and to the left if $\dot{x} < 0$.
"""

# ╔═╡ 42c28b16-2218-11eb-304a-e534353fa12b
md"""
This vector field indeed gives us a *qualitative* picture of the dynamics. It does not tell us how fast the dynamics will occur in each region, but it indicates what the *tendency* is. We have coded the fixed points according to their stability; this may be calculated using the derivative evaluated at the fixed point, $f'(x^*)$, since this derivative controls the behaviour of nearby initial conditions $x^* + \delta x$.
The unstable fixed point is shown as a green square, and the stable fixed point as a grey circle; we will use this convention throughout the notebook.
"""

# ╔═╡ e6ed3366-2246-11eb-3a86-f375162c746d
md"""
## Bifurcations
"""

# ╔═╡ ec3c5d9c-2246-11eb-1ec5-e3d9c8a7fa23
md"""
Now suppose that there is a **parameter** $\mu$ in the system that can be varied. For each value of $\mu$ we have a *different* ODE

$$\dot{x} = f_\mu(x).$$

For example,
$$\dot{x} = \mu + x^2.$$

Let's draw the state space for each different value of $\mu$:
"""

# ╔═╡ 3fa21738-2247-11eb-09d5-4dcf0e7377fc
g(μ, x) = μ + x^2

# ╔═╡ 19eb987e-2248-11eb-2982-1930b5840aed
@bind λ Slider(-1.0:0.05:1, show_value=true)

# ╔═╡ 49ec803a-22d7-11eb-3f04-db8f8c119a7d
md"""
Now let's collect all the vector fields into a single plot. We *rotate* the vector field to now be vertical, thinking of the dynamics of $x$ as occurring along the vertical direction. The horizontal axis now represents the different possible values of the parameter $\mu$:
"""

# ╔═╡ 27190ac0-230b-11eb-3367-af2bbbf57e5e
md"""
We see that at the **critical value** $\mu_c = 0$ there is a **qualitative change in behaviour** in the system: for $\mu_c < 0$ there are two fixed points, whereas for $\mu_c > 0$ there are no fixed points at all.  In this particular ODE the two fixed points collide in a **saddle--node** or **fold** bifurcation.
"""

# ╔═╡ 48ef665e-22d0-11eb-2c72-b722fafd296e
md"""
## Bistability and hysteresis
"""

# ╔═╡ 51eeb46c-22d0-11eb-1643-c92e725d9da0
md"""
Now let's look at the dynamics of the following system:

$$\dot{x} = \mu + x - x^3.$$
"""

# ╔═╡ 44c1738a-2247-11eb-3232-f30b46fd7ce6
h(μ, x) = μ + x - x^3

# ╔═╡ a56ad0bc-22d7-11eb-0418-d36f3ef14109
md"""
We see that there is a range of values of $\mu$ for which there are *three coexisting fixed points*, two stable and one unstable. Since there are two stable fixed points in which the system can remain, we say that the system is **bistable**.
"""

# ╔═╡ 5b9c66ce-224b-11eb-02a4-b9a5a173ee6a
md"""
Now that we understand what the plots mean and the dynamics, let's plot just the fixed points $x^*(\mu)$ as a function of $\mu$. Such a plot is called a **bifurcation diagram**:
"""

# ╔═╡ aa81663c-230b-11eb-3c0e-39759f36deb6
md"""
The pieces of curve are called **branches**.
"""

# ╔═╡ 55d4d188-2249-11eb-0c47-dfddce9b2027
md"""
## Hysteresis
"""

# ╔═╡ 5ac33518-2249-11eb-2d57-7de185d43e5f
md"""
Suppose we now think about slowly varying the parameter $\mu$. If we change the parameter $\mu$ by a little, the system is no longer at a fixed point, since the position of the fixed point moves when $\mu$ changes. However, the system will theb **relax** by following the dynamics at the new value of $\mu$, and will rapidly converge to the new fixed point nearby.

For example, starting at $\mu=-2$, the system will stay on the lower black (stable) **branch** until $\mu=0.4$ or so. At that point, two fixed points collide and annihilate each other! After that there is no longer a fixed point nearby. However, there is another fixed point much further up that will now attract all trajectories, so the system rapidly transitions to that fixed point.

Now suppose we decrease the parameter again. The system will now track the *upper* branch until $\mu=-0.4$ or so, when again it will jump back down.

For each parameter value $\mu$ in the interval $[-0.4, 0.4]$ there is **bistability**, i.e. **coexistence** of *two* fixed points with the same value of $\mu$ (together with a third, unstable fixed point that is not observable).

The fact that the system tracks different stable branches depending on where we started, i.e. on the history, is known as **hysteresis**.
"""

# ╔═╡ 4c73705e-230c-11eb-3c90-b14536d78808
md"""
Hysteretic behaviour like this is found in many scientific and engineering contexts,  including switches in biology, for example genetic switches, and in the historical dynamics of the earth's climate.
"""

# ╔═╡ 2ced9e26-22d0-11eb-34a5-bf80822f5c43
md"""
## Slow--fast systems
"""

# ╔═╡ 3213d406-22d0-11eb-0e89-71b42ea5cada
md"""
What are we actually doing when we let the parameter $\mu$ vary? Effectively we now have a system with *two* equations, for example

$$\dot{x} = \mu + x - x^3;$$
$$\dot{\mu} = \epsilon,$$

where $\mu$ varies at some slow speed $\epsilon$. On a time scale much shorter than $1 / \epsilon$, the dynamics of $x$ "does not know" that $\mu$ is changing, so it will converge to a fixed point $x^*(\mu)$ for the current value of $\mu$. [An associated term is **adiabatic approximation**.] However, $\mu$ does gradually change, so the value of $x$ will effectively "slide along" the curve $x(t) \simeq x^*(\mu(t))$, tracking the curve of fixed points as $\mu$ changes.

Once $\mu$ reaches a critical value $\mu_c$, however, there is no longer a nearby fixed point, and the dynamics will rapidly transition to the far away alternative fixed point.

If we now reverse the dynamics of $\mu$, we slide back along the upper branch.
"""

# ╔═╡ d0918f0c-22d4-11eb-383a-553e33114c52
md"""
## Function library
"""

# ╔═╡ 319c2ba8-220e-11eb-1ee0-4dd51ff4cd25
euler_step(f, x, h) = x + h * f(x)

# ╔═╡ 42de6dd6-220e-11eb-2a48-a3c69b10a03b
function euler(f, x0, h, t_final)
	ts = [0.0]
	xs = [x0]

	x = x0
	t = 0.0

	while t < t_final
		x = euler_step(f, x, h)
		t += h

		push!(xs, x)
		push!(ts, t)
	end

	return (ts=ts, xs=xs)  # a named tuple
end

# ╔═╡ 6e90ea26-220e-11eb-0c65-bf52b3d2e195
results = euler(logistic, 0.5, 0.01, 20.0)

# ╔═╡ 96e22792-220e-11eb-2729-63964507b5f2
begin
	plot(results.ts, results.xs, size=(400, 300), leg=false, xlabel=L"t", ylabel=L"x(t)", lw=3)
	scatter!([(results.ts[1], results.xs[1])])
	ylims!(0.4, 1.1)
end

# ╔═╡ 2300b29a-22d5-11eb-3c99-bdec0e5a2685
let
	p = plot(xlabel=L"t", ylabel=L"x(t)", size = (400, 300), leg=false, ylim=(-1, 2))

	results = euler(logistic, x0, 0.01, 10.0)

	plot!(results.ts, results.xs, alpha=1, lw=3)
	scatter!([0], [x0])

	hline!([0], ls=:dash)
	p
end

# ╔═╡ f1d7fd8e-220e-11eb-0d2a-cf9e829c2bb5
let
	p = plot(xlabel=L"t", ylabel=L"x(t)", size = (400, 300), leg=false, ylim=(-1, 2))

	for x0 in -0.5:0.1:2.0
		results = euler(logistic, x0, 0.01, 10.0)

		plot!(results.ts, results.xs, alpha=0.5, lw=2)
	end

	p
end

# ╔═╡ 546168ce-2218-11eb-198b-1da9f8a9a242
derivative(f, x, h=0.001) = ( f(x + h) - f(x - h) ) / (2h)

# ╔═╡ 46f1e3a2-2217-11eb-184f-f5649b892004
"Draw 1D vector field using centred arrows"
function horiz_vector_field(f)

	xlo, xhi = -2, 2.5
	arrow_size = 0.07
	tol=1e-5  # tolerance to check for fixed point

	p = plot(size=(400, 100), xlim=(xlo, xhi), ylim=(-0.5, 0.5), leg=false, yticks=[])

	for x in xlo:0.2:xhi
		d = f(x)  # derivative

		if d > tol
			plot!([x-arrow_size, x+arrow_size], [0, 0], arrow=true, c=:blue, alpha=0.5, lw=2)

		elseif d < -tol
			plot!([x+arrow_size, x-arrow_size], [0, 0], arrow=true, c=:red, alpha=0.5, lw=2)

		end

	end

	roots = find_zeros(f, -10, 10)

	for root in roots
		stability = sign(derivative(f, root))

		if stability > 0
			scatter!([root], [0], c=:green, alpha=0.5, m=:square, ms=5)  # unstable
		else
			scatter!([root], [0], c=:black, alpha=0.5, ms=6)  # stable
		end

	end

	p
end

# ╔═╡ 2b1b49b0-2247-11eb-1d25-57f3b4f46d04
"Draw vertical vector field of 1D ODE"
function vector_field!(p, μ, f)

	xlo, xhi = -2, 2
	arrow_size = 0.07
	tol=1e-5  # tolerance to check for fixed point

	for x in xlo:0.2:xhi
		d = f(x)  # derivative

		if d > tol

			plot!([μ, μ], [x-arrow_size, x+arrow_size], arrow=(3.0, 2.0), arrowstyle=:triangle, c=:blue, alpha=0.4, lw=2)

		elseif d < -tol
			plot!([μ, μ], [x+arrow_size, x-arrow_size], arrow=true, c=:red, alpha=0.4, lw=2)
		end

	end

	roots = find_zeros(f, -10, 10)

	for root in roots
		stability = sign(derivative(f, root))

		if stability > 0
			scatter!([μ], [root], c=:green, alpha=0.5, m=:square, ms=4)  # unstable
		else
			scatter!([μ], [root], c=:black, alpha=0.5, ms=5)  # stable
		end

	end

	p
end

# ╔═╡ d4fa7bec-2371-11eb-3e9f-69d552546331
function vector_field(f)
	p = plot(leg=false, xticks=[], size=(100, 300), xrange=(-1, 1))
	vector_field!(p, 0, f)
end

# ╔═╡ 9a833edc-2217-11eb-0701-99862b410bfa
vector_field(logistic)

# ╔═╡ 25578fa6-2248-11eb-2838-57fbfb928fc4
vector_field(x -> g(λ, x))

# ╔═╡ e55d2780-2247-11eb-01e0-fbdc94bba264
function bifurcation_diagram(h)
	p = plot(leg=false, ratio=1)

	for μ in -2:0.15:2
		vector_field!(p, μ, x -> h(μ, x))
	end

	xlabel!("μ")
	ylabel!("fixed points and dynamics with given μ")

	return p
end

# ╔═╡ b9ee4db4-22d7-11eb-38ef-511e30df16cc
bifurcation_diagram(g)

# ╔═╡ f040d1c4-2247-11eb-0bf0-090e90a86a85
bifurcation_diagram(h)

# ╔═╡ e27dace4-2248-11eb-2ae1-953639d8c944
function fixed_points(f)
	p = plot(leg=false, ratio=1)

	for μ in -2:0.01:2

		roots = find_zeros(x -> f(μ, x), -10, 10)

		for root in roots
			stability = sign(derivative(x -> f(μ, x), root))

			if stability > 0
				scatter!([μ], [root], c=:green, alpha=0.5, m=:square)  # unstable
			else
				scatter!([μ], [root], c=:black, alpha=0.5)  # stable
			end
		end

	end

	xlabel!(L"\mu")
	ylabel!(L"x^*(\mu)")

	return p
end

# ╔═╡ 0c3c88fe-2249-11eb-03af-bf9c393fadd4
fixed_points(h)

# ╔═╡ Cell order:
# ╠═88b46d2e-220e-11eb-0f7f-b3f523f0214e
# ╟─6f554dc0-220c-11eb-341f-1b66de841c86
# ╟─816de40c-220c-11eb-2d3a-23e6267cd529
# ╟─d21f6358-220c-11eb-00b8-03fe0746fdc9
# ╟─735bb47e-220d-11eb-1ba6-8d591a935857
# ╟─79b148b6-220d-11eb-2785-25d05068aeed
# ╟─d43a869e-220d-11eb-009c-a119de57e7da
# ╟─088eabb2-220e-11eb-1ac0-df419b87f39a
# ╟─a44cc180-22d4-11eb-2129-211d024c4921
# ╠═6aa1d2f4-220e-11eb-06e0-c346f74aa018
# ╟─b6d0dd18-22d4-11eb-2083-0b02de030579
# ╠═6e90ea26-220e-11eb-0c65-bf52b3d2e195
# ╟─f055be76-22d4-11eb-268a-bf70c7d8f1a1
# ╟─96e22792-220e-11eb-2729-63964507b5f2
# ╟─90ccb392-2216-11eb-1fd8-83b7d7c16b54
# ╟─bef9c2e4-220e-11eb-24d8-bd618d2985ea
# ╠═17317eb8-22d5-11eb-0f46-2bf1bddd36bb
# ╟─2300b29a-22d5-11eb-3c99-bdec0e5a2685
# ╟─ae862e10-22d5-11eb-3d75-1d748cf86944
# ╟─f1d7fd8e-220e-11eb-0d2a-cf9e829c2bb5
# ╟─446c564e-220f-11eb-191a-6b419e790f3f
# ╟─38894c32-2210-11eb-26ff-d9796edd7871
# ╟─a196a20a-2216-11eb-11f9-6789bb2a3ece
# ╟─c15a53e8-2216-11eb-1460-8510aea0f570
# ╠═9a833edc-2217-11eb-0701-99862b410bfa
# ╟─42c28b16-2218-11eb-304a-e534353fa12b
# ╟─e6ed3366-2246-11eb-3a86-f375162c746d
# ╟─ec3c5d9c-2246-11eb-1ec5-e3d9c8a7fa23
# ╠═3fa21738-2247-11eb-09d5-4dcf0e7377fc
# ╠═19eb987e-2248-11eb-2982-1930b5840aed
# ╠═25578fa6-2248-11eb-2838-57fbfb928fc4
# ╟─49ec803a-22d7-11eb-3f04-db8f8c119a7d
# ╠═b9ee4db4-22d7-11eb-38ef-511e30df16cc
# ╟─27190ac0-230b-11eb-3367-af2bbbf57e5e
# ╟─48ef665e-22d0-11eb-2c72-b722fafd296e
# ╟─51eeb46c-22d0-11eb-1643-c92e725d9da0
# ╠═44c1738a-2247-11eb-3232-f30b46fd7ce6
# ╠═f040d1c4-2247-11eb-0bf0-090e90a86a85
# ╟─a56ad0bc-22d7-11eb-0418-d36f3ef14109
# ╟─5b9c66ce-224b-11eb-02a4-b9a5a173ee6a
# ╠═0c3c88fe-2249-11eb-03af-bf9c393fadd4
# ╟─aa81663c-230b-11eb-3c0e-39759f36deb6
# ╟─55d4d188-2249-11eb-0c47-dfddce9b2027
# ╟─5ac33518-2249-11eb-2d57-7de185d43e5f
# ╟─4c73705e-230c-11eb-3c90-b14536d78808
# ╟─2ced9e26-22d0-11eb-34a5-bf80822f5c43
# ╟─3213d406-22d0-11eb-0e89-71b42ea5cada
# ╟─d0918f0c-22d4-11eb-383a-553e33114c52
# ╟─319c2ba8-220e-11eb-1ee0-4dd51ff4cd25
# ╟─42de6dd6-220e-11eb-2a48-a3c69b10a03b
# ╟─546168ce-2218-11eb-198b-1da9f8a9a242
# ╟─46f1e3a2-2217-11eb-184f-f5649b892004
# ╟─2b1b49b0-2247-11eb-1d25-57f3b4f46d04
# ╟─d4fa7bec-2371-11eb-3e9f-69d552546331
# ╟─e55d2780-2247-11eb-01e0-fbdc94bba264
# ╟─e27dace4-2248-11eb-2ae1-953639d8c944

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
# solving the wave equation
# \frac{\partial ^2 u}{\partial t ^2} - c^2 [ \frac{\partial ^2 u}{\partial z^2} + \frac{\partial ^2 u}{\partial x^2}] + f = 0

import DPFEHM
#using Printf
import Optim
import PyPlot
import Zygote
## FUNCTIONS ####################

# generate source wavelet (Ricker wavelet)
function ricker(f, t; Ã=0)  # generate ricker wavelet
    return (1 .- 2 .* pi.^2 .* f.^2 .* (t.-Ã).^2) .* exp.(-pi.^2 .* f.^2 .* (t.-Ã).^2)
end

# get norm to generate objective function
function getNorm(c, f, nz, nx, nt, dt, ridx, ur)
    usave = DPFEHM.getuIters(c, f, nz, nx, nt, dz, dx, dt)
    l2norm = 0
    for i=1:length(ridx)
        l2norm += sum((ur[i, :] .- usave[ridx[i], :]).^2)
    end
    c_matrix = reshape(c, nz, nx)
    grad_reg = sum((c_matrix[2:end, :] - c_matrix[1:end - 1, :]) .^ 2) + sum((c_matrix[:, 2:end] - c_matrix[:, 1:end - 1]) .^ 2)
    cbin = (c .- 2000) ./ 500
    clow = c .< 2250
    chigh = c .>= 2250
    bin_reg = (sum((cbin .* clow) .^ 2) + sum((cbin .* (chigh .- 1)) .^ 2)) / (4 * length(c))
    return l2norm / 1e-6 + grad_reg / (1e4 * 500 ^ 2) + 10 * bin_reg
end

# plot shot record
function plotShot(r, nt, dt, nr)
    w, h = figaspect(1.6)
    figure(figsize=[1.5*w, 1.5*h])
    boundary = 0.95*maximum(abs.(r))
    imshow(r', extent=[1, nr, dt*nt, 0], vmin=-boundary, vmax=boundary, aspect="auto", cmap="seismic")
    cb = colorbar()
    cb.set_label("Amplitude")
    ylabel("Time (s)")
    xlabel("Receiver number")
    title("Shot record")
end

## PROBLEM SETUP ####################
# NOTE: z is depth, x is horizontal distance. Indexing is (z, x)
# time steps
nt = 500    # number of time steps (unitless)
dt = 0.001  # how long each time step is (seconds)

# spatial steps
nz = 100    # number of spatial steps in z (unitless)
nx = 100    # number of spatial steps in x (unitless)
dz = 5      # how long each spatial step is in z (meters)
dx = 5      # how long each spatial step is in x (meters)

#z = range(0, length=nz, step=dz)   # can be helpful for plotting
#x = range(0, length=nx, step=dx)   # can be helpful for plotting

# wavelet
t = range(0, stop=nt*dt, length=nt)
wave = 1000 .*ricker(45, t; Ã=0.06)

# create forcing term (ricker wavelet)
f = zeros(nz, nx, nt)
halfx = floor(Int, nx/2)
halfz = floor(Int, nz/2)
f[1, halfx, :] = wave
f = reshape(f, nz * nx, nt)

# get starting (constant background) velocity model
c = 2000 .*ones(nz, nx) # m/s
c = reshape(c, nz * nx)

# truth velocity model (scatterer in center of domain)
ctrue = 2000 .*ones(nz, nx) # m/s
#ctrue[halfz, halfx] = 2800 # m/s
sz = 10 # scatterer size
ctrue[halfz-sz:halfz+sz, halfx-sz:halfx+sz] .= 2500 # m/s
#ctrue[halfz-sz-25:halfz+sz-25, halfx-sz:halfx+sz] .= 3200 # m/s
ctrue = reshape(ctrue, nz * nx)

# receiver locations along surface
xloc = Array(1:nx)
zloc = 5 .*ones(size(xloc))
ridx = floor.(Int, DPFEHM.linearIndex.(nz.*ones(size(xloc)), zloc, xloc))

# get wavefield in background medium (can comment out)
@time u = DPFEHM.getuIters(c, f, nz, nx, nt, dz, dx, dt)
fig, axs = PyPlot.subplots(2, 13, figsize=(30, 5))
axs = permutedims(axs, [2, 1])
j = 1
#for i = 2:10:size(u, 2)
for i = 2:20:size(u, 2)
    global j
    axs[j].imshow(reshape(u[:, i], nz, nx), cmap="seismic")
    j += 1
end
fig.tight_layout()
display(fig)
println()
PyPlot.close(fig)

# get wavefield in truth velocity model
@time utrue = DPFEHM.getuIters(ctrue, f, nz, nx, nt, dz, dx, dt)
fig, axs = PyPlot.subplots(2, 13, figsize=(30, 5))
axs = permutedims(axs, [2, 1])
j = 1
#for i = 2:10:size(u, 2)
for i = 2:20:size(utrue, 2)
    global j
    axs[j].imshow(reshape(utrue[:, i], nz, nx), cmap="seismic")
    j += 1
end
fig.tight_layout()
display(fig)
println()
PyPlot.close(fig)

# shot record (data recorded at receivers)
shotRecord = utrue[ridx, :] # get the shot record

# gradient functions
g3(c, f) = getNorm(c, f, nz, nx, nt, dt, ridx, shotRecord)
g3(c) = getNorm(c, f, nz, nx, nt, dt, ridx, shotRecord)

# Full waveform inversion (with gradient descent)
niter = 3                           # number of iterations
csave = zeros(size(c, 1), niter)    # save c at every iteration
csave[:, 1] = c                     # initial c at first iteration
stepLength=1e10 .*ones(niter)       # step length
l2save = zeros(niter)               # save loss function at every iteration

# update model
#for i = 1:2 # for one iteration (aka RTM)
#=
for i = 1:niter-1
    l2save[i] = g3(csave[:, i])
    @show i
    grad = Zygote.gradient(g3, csave[:, i])
    @show extrema(grad[1])
    csave[:, i+1] = csave[:, i] .- stepLength[i].*reshape(grad[1], nz*nx)
    #@printf("Done with iteration %i\n", i)
end
=#

function grad!(storage, c)
    storage .= Zygote.gradient(g3, c)[1]
    return storage
end
opt = Optim.optimize(g3, grad!, c, Optim.LBFGS(), Optim.Options(show_trace=true, iterations=10))
#opt = Optim.optimize(g3, grad!, ctrue, Optim.LBFGS(), Optim.Options(show_trace=true, iterations=5))
csave[:, 2] = opt.minimizer
fig, axs = PyPlot.subplots(1, 2)
axs[1].imshow(reshape(ctrue, nz, nx))
axs[2].imshow(reshape(opt.minimizer, nz, nx))
display(fig)
println()
PyPlot.close(fig)
fig, axs = PyPlot.subplots(1, 2)
axs[1].imshow(reshape(ctrue, nz, nx), vmin=0.9 * minimum(ctrue), vmax=1.1 * maximum(ctrue))
axs[2].imshow(reshape(opt.minimizer, nz, nx), vmin=0.9 * minimum(ctrue), vmax=1.1 * maximum(ctrue))
display(fig)
println()
PyPlot.close(fig)
uopt = DPFEHM.getuIters(opt.minimizer, f, nz, nx, nt, dz, dx, dt)
fig, axs = PyPlot.subplots(10, 10, figsize=(10, 10), sharex=true, sharey=false)
for i in 1:length(ridx)
    axs[i].plot(shotRecord[i, :], label="true", alpha=0.5)
    axs[i].plot(uopt[ridx[i], :], label="inverse result", alpha=0.5)
end
axs[1].legend()
fig.tight_layout()
display(fig)
println()
PyPlot.close(fig)
@show sum((uopt[ridx, :] .- shotRecord) .^ 2)

ext: jl


text:
"""
Computes shortest distance of every node from a given source.

Input:
-----
* G::SparseMatrixCSC
* root::Int

Output:
------
* sssp_dist::Vector{Int32}

"""
function sssp(G::SparseMatrixCSC, root::Int)
	rows = G.colptr - 1
    cols = G.rowval - 1
    nodes = length(rows) - 1
    edges = length(cols)
    rows = map(Int32, rows)
    cols = map(Int32, cols)
	vals = ones(Cint, edges)
    sssp_dist = Array(Cuint, nodes)
	ccall((:sssp, gunrock_lib), Void, (Ptr{Cuint}, Int, Int, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Int), sssp_dist, nodes, edges, rows, cols, vals, root - 1)
	sssp_dist
end


text:
#include("Init.jl")
include("Agents/AbstractTypes/AbstractType.jl")
import PyPlot
using Statistics
using Random

mutable struct Bid <: Container
    bidPrice::Array{Float64,1}
    bidQuantity::Array{Float64,1}
    bidCondition::Symbol #Accepted, unaccepted
    expectedProfit::Float64
    expectedProfit_Type::Array{Float64,1}
    attempt::Int64
    lastProfit::Float64

    Update::Function

    function Bid(bidPrice::Array{Float64,1},unitCapacity::Float64)
        self = new()
        self.bidPrice = bidPrice
        self.bidQuantity = unitCapacity .* ones(length(self.bidPrice))
        self.expectedProfit = 0.0
        self.expectedProfit_Type = zeros(length(self.bidPrice))
        self.attempt = 0
        self.lastProfit = 0.0

        self.Update = function(typeProfit::Array{Float64,1},learningRate::Float64)
            newProfit = sum(typeProfit)
            self.lastProfit = sum(typeProfit)
            self.expectedProfit_Type = self.expectedProfit_Type .* (1-learningRate) .+ typeProfit .* learningRate
            self.expectedProfit = self.expectedProfit * (1-learningRate) + newProfit * learningRate
        end

        return self
    end
end

mutable struct TestTechnology <: Technology
    type::Symbol
    marginalCost::Float64
    commitment::Bool
    offerPrice::Float64
    offerQuantity::Float64
    acceptedQuantity::Float64
    owner::ActiveAgent

    function TestTechnology(type::Symbol,marginalCost::Float64)
        self = new()
        self.type = type
        self.marginalCost = marginalCost
        self.commitment = true
        self.offerPrice = 0.0
        self.offerQuantity = 0.0
        self.acceptedQuantity = 0.0

        return self
    end
end

mutable struct TestGenCo <: ActiveAgent
    name::String
    bidSpace::Array{Bid,1}
    commitment::Array{Bool,1}
    technologies::Array{Technology,1}
    baseBid::Bid
    lastBid::Bid

    Learn::Function
    SelectBid::Function
    SubmitBid::Function

    function TestGenCo(technologies::Array{TestTechnology,1},testGenCoName::String)
        self = new()
        self.name = testGenCoName
        self.technologies = technologies
        self.bidSpace = Array{Bid,1}(undef,0)
        self.commitment = [technology.commitment for technology in technologies]
        baseBidPrice = zeros(length(self.technologies))
        for techCount in 1:length(self.technologies)
            baseBidPrice[techCount] = self.technologies[techCount].marginalCost
        end
        self.baseBid = Bid(baseBidPrice,800.0)

        allMarkupLevel = collect(Iterators.product(ntuple(i-> 0:2:10, length(self.technologies))...))[:]
        for markupLevel in allMarkupLevel
            newBidPrice = deepcopy(baseBidPrice)
            [newBidPrice[i] = newBidPrice[i] *( markupLevel[i] *0.1 + 1) for i = 1:length(self.technologies)]
            #newBidPrice = newBidPrice .* (1 .+ newBidPrice*0.1)
            newBid = Bid(newBidPrice,800.0)
            push!(self.bidSpace,newBid)
        end

        self.Learn = function(newPrice::Float64,learningRate::Float64)
            unitProfit = (newPrice .- self.baseBid.bidPrice) #.* Int.(self.commitment)
            acceptedQuantity = [technology.acceptedQuantity for technology in self.technologies]
            #unitProfit = max.(newPrice .- self.baseBid.bidPrice,0)
            #println(unitProfit)
            #newProfit = sum(unitProfit .* acceptedQuantity)
            typeProfit = unitProfit .* acceptedQuantity
            self.lastBid.Update(typeProfit,learningRate)
            sort!(self.bidSpace, by = x -> x.expectedProfit,rev = true)
            [technology.acceptedQuantity = 0.0 for technology in self.technologies]
            #return newProfit
        end

        self.SelectBid = function(explorationRate::Float64)
            randomNumber = rand(1)[1]
            bidSelect = self.bidSpace[1]
            if randomNumber < explorationRate
                bidSelect = self.bidSpace[rand(2:length(self.bidSpace))]
            end
            bidSelect.attempt += 1
            return bidSelect
        end

        self.SubmitBid = function(bid::Bid)
            for techCount in 1:length(self.technologies)
                self.technologies[techCount].offerPrice = bid.bidPrice[techCount]
                self.technologies[techCount].offerQuantity = bid.bidQuantity[techCount]
            end
            self.lastBid = bid
        end

        return self
    end
end

#MO,allGenCos,systemLoad,allPowerFlows,technologyPool,systemState,systemHistory = InitSystem(systemSettings,genCoHorizon)

#testGenCo3 = TestGenCo(deepcopy([baseTech,midTech,peakTech,peakTech]))

function ClearMarket(allTechnologies::Array{TestTechnology,1}, load::Float64)
    onlineTechnologies = allTechnologies
    #NOTE:start clear the market
    marketPrice = 0
    if length(onlineTechnologies) == 0
        println("No Technology!")
        marketPrice = 3000
    end

    if length(onlineTechnologies) > 0
        sort!(onlineTechnologies, by = x -> x.offerPrice)
        generatorCount = 1 #initialize generator number
         #this is a variable indicating the energy gap of one consumer when switching generators
        maxLoad = deepcopy(load)
        residualTechnology = 0

        while generatorCount <= length(onlineTechnologies)
            thisQuantity = onlineTechnologies[generatorCount].offerQuantity
            maxLoad += - thisQuantity
            #onlineTechnologies[generatorCount].commitment = true
            onlineTechnologies[generatorCount].acceptedQuantity = thisQuantity
            if maxLoad <= 0
                residualTechnology = -maxLoad
                marketPrice = onlineTechnologies[generatorCount].offerPrice
                onlineTechnologies[generatorCount].acceptedQuantity += maxLoad
                break
            end
            generatorCount += 1
        end
        #println(maxLoad)
        #println("Generator count ", generatorCount)
        if maxLoad > 0
            marketPrice = 3000
            residualTechnology = 0
        end
    end
    return marketPrice,onlineTechnologies[generatorCount].type,marketPrice/onlineTechnologies[generatorCount].marginalCost
end

function PlotTrade(allTestGenCos::Array{TestGenCo,1},load::Float64)
    allTechnologies = Array{TestTechnology,1}(undef,0)
    for testGenCo in allTestGenCos
        for technology in testGenCo.technologies
            push!(allTechnologies,deepcopy(technology))
        end
    end
    sort!(allTechnologies, by = x -> x.offerPrice)
    [println(technology.offerPrice) for technology in allTechnologies]
    allGenCoNames = [string(technology.owner.name) for technology in allTechnologies]
    println(allGenCoNames)
    allOfferPrice = [technology.offerPrice for technology in allTechnologies]
    allMarginalCost = [technology.marginalCost for technology in allTechnologies]
    allMarkupLevel = [technology.offerPrice/technology.marginalCost for technology in allTechnologies]
    allOfferQuantity = [technology.offerQuantity for technology in allTechnologies]
    allOfferQuantity_Cum = cumsum(allOfferQuantity)
    pushfirst!(allOfferQuantity_Cum,0.0)
    pushfirst!(allOfferPrice,allOfferPrice[1])
    pushfirst!(allMarginalCost,allMarginalCost[1])

    println(allOfferPrice)
    println(allMarginalCost)

    allTechNames = [string(technology.type,"(",round((technology.offerPrice/technology.marginalCost-1),digits = 1)*100,"%)") for technology in allTechnologies]
    println(allTechNames)
    PyPlot.step(allOfferQuantity_Cum, allOfferPrice, linestyle="-",label = "Bid")
    PyPlot.step(allOfferQuantity_Cum, allMarginalCost, linestyle=":",label = "Marginal cost")
    PyPlot.plot(load .* ones(100), collect(1:100),label = "Load")
    for txt_tech in zip(allOfferQuantity_Cum[1:end-1],allOfferQuantity_Cum[2:end], allOfferPrice[2:end], allTechNames)
        PyPlot.text(x=(txt_tech[1]+txt_tech[2])/2, y=txt_tech[3]+1, s=txt_tech[4])
    end

    for txt_tech in zip(allOfferQuantity_Cum[1:end-1],allOfferQuantity_Cum[2:end], allOfferPrice[2:end], allGenCoNames)
        PyPlot.text(x=(txt_tech[1]+txt_tech[2])/2, y=txt_tech[3]+5, s=txt_tech[4])
    end

    PyPlot.xlim([0,sum(allOfferQuantity)])
    PyPlot.xlabel("Load [MW]",fontsize=16)
    PyPlot.ylabel("Price [€/MWh]",fontsize=16)
    PyPlot.xticks(fontsize=16)
    PyPlot.yticks(fontsize=16)
    PyPlot.legend(fontsize = 16)
    PyPlot.grid(true)
end

function RunBidding(allTestGenCos::Array{TestGenCo,1},totalEpisode::Int64,load::Float64)
    profitHistory = zeros(length(allTestGenCos),totalEpisode)
    marginalMarkupHistory = zeros(length(allTestGenCos),totalEpisode)
    marketPriceHistory = zeros(totalEpisode)
    for episode in 1:totalEpisode
        progressIndex = 1-episode/totalEpisode
        explorationRate = progressIndex^(1/20)
        exploitationRate = 1-explorationRate
        allTechnologies = Array{TestTechnology,1}(undef,0)
        for testGenCo in allTestGenCos
            testGenCo.SubmitBid(testGenCo.SelectBid(explorationRate))
            append!(allTechnologies,testGenCo.technologies)
        end

        shuffle!(allTechnologies)
        #if changeSequence[episode] == true
        #    reverse!(allTechnologies)
        #end
        #=
        randomNumber = rand(1)[1]
        if randomNumber > 0.5
            reverse!(allTechnologies)
        end
        =#

        #testGenCo3.SubmitBid(testGenCo3.SelectBid(explorationRate))
        #allTechnologies = vcat(testGenCo1.technologies,testGenCo2.technologies)
        #allTechnologies = testGenCo1.technologies
        marketPrice,marginalGenerator,markup = ClearMarket(allTechnologies, load)

        marketPriceHistory[episode] = marketPrice
        #println(marketPrice)

        for genCoCount in 1:length(allTestGenCos)
            marginalMarkupHistory[genCoCount,episode] = allTestGenCos[genCoCount].bidSpace[1].bidPrice[end]/allTestGenCos[genCoCount].baseBid.bidPrice[end]
            allTestGenCos[genCoCount].Learn(marketPrice,0.1)
            profitHistory[genCoCount,episode] = allTestGenCos[genCoCount].bidSpace[1].expectedProfit
        end
        #testGenCo3.Learn(marketPrice,0.05)

        #profitHistory_3[episode] = testGenCo3.bidSpace[1].expectedProfit
    end
    optimalGenCos = deepcopy(allTestGenCos)
    allTechnologies_optimal = Array{TestTechnology,1}(undef,0)
    for testGenCo in optimalGenCos
        testGenCo.SubmitBid(testGenCo.bidSpace[1])
        append!(allTechnologies_optimal,testGenCo.technologies)
    end
    optimalMarketPrice,marginalGenerator,markup = ClearMarket(allTechnologies_optimal, load)
    return optimalMarketPrice,marginalGenerator,markup,profitHistory,marginalMarkupHistory,marketPriceHistory
end

#baseTech = TestTechnology(:Base, 12.5)
#midTech = TestTechnology(:Mid,21.5)
#peakTech = TestTechnology(:Peak,45.67)

baseTech = TestTechnology(:Base, 10.0)
midTech = TestTechnology(:Mid,30.0)
peakTech = TestTechnology(:Peak,70.0)
testTechnologyPool = [baseTech,midTech,peakTech]

testGenCo1 = TestGenCo(deepcopy([baseTech,midTech,peakTech]),"GenCo1")
testGenCo2 = TestGenCo(deepcopy([baseTech,midTech,peakTech]),"GenCo2")
testGenCo3 = TestGenCo(deepcopy([baseTech,midTech,peakTech]),"GenCo3")
testGenCo4 = TestGenCo(deepcopy([baseTech,midTech,peakTech]),"GenCo4")
testGenCo5 = TestGenCo(deepcopy([baseTech,midTech,peakTech]),"GenCo5")
#allTestGenCos = [testGenCo1,testGenCo2,testGenCo3,testGenCo4,testGenCo5]
allTestGenCos = [testGenCo1,testGenCo2,testGenCo3]
for testGenCo in allTestGenCos
    for technology in testGenCo.technologies
        technology.owner = testGenCo
    end
end
totalInstalledCapacity = length(allTestGenCos) * 3 * 800.0
loadList = collect(100.0:200.0:(totalInstalledCapacity-100))

newAllTestGenCos = deepcopy(allTestGenCos)

totalEpisode = 100000
simulationTimes = 20

optimalMarketPriceList = Array{Real,1}(undef,0)
marginalGeneratorList = Array{Symbol,1}(undef,0)

profitHistory_all = zeros(length(allTestGenCos),totalEpisode*simulationTimes)
marginalMarkupHistory_all = zeros(length(allTestGenCos),totalEpisode*simulationTimes)
marketPriceHistory_all = zeros(totalEpisode*simulationTimes)
averageMarkupHistory_all = Array{Float64,1}(undef,0)
allTestGenCos_history = Array{Array{TestGenCo,1},1}(undef,0)

for load in loadList
    println(load)
    allTestGenCos = deepcopy(newAllTestGenCos)
    markupList = Array{Float64,1}(undef,0)
    randomSeedList = rand(1:1000,simulationTimes)
    for randomSeedCount in 1:length(randomSeedList)
        Random.seed!(randomSeedList[randomSeedCount])
        println(randomSeedCount)
        optimalMarketPrice,marginalGenerator,markup,profitHistory,marginalMarkupHistory,marketPriceHistory = RunBidding(allTestGenCos,totalEpisode,load)
        thisStart = (randomSeedCount-1)*totalEpisode+1
        thisStop = randomSeedCount * totalEpisode
        profitHistory_all[:,thisStart:thisStop] = profitHistory
        marginalMarkupHistory_all[:,thisStart:thisStop] = marginalMarkupHistory
        marketPriceHistory_all[thisStart:thisStop] = marketPriceHistory
        push!(optimalMarketPriceList,optimalMarketPrice)
        push!(marginalGeneratorList,marginalGenerator)
        push!(markupList,markup)
    end
    push!(allTestGenCos_history,allTestGenCos)
    averageMarkup = Statistics.mean(markupList.-1)
    push!(averageMarkupHistory_all,averageMarkup)
end

PyPlot.figure(6)
PyPlot.plot(averageMarkupHistory_all.*100,marker = "o",label = "markup percentage")

PyPlot.xticks(collect(0:(length(averageMarkupHistory_all)-1)),loadList,rotation = 80)
PyPlot.title("Markup percentage vs load level")
PyPlot.xlabel("Load (MW)",fontsize=16)
PyPlot.ylabel("Markup percentage [%]",fontsize=16)

PyPlot.yticks(fontsize=16)
PyPlot.grid(true)
loadDivisionList = (loadList.+100) .% 800.0 .== 0
for loadDivision in collect(1:length(averageMarkupHistory_all))[loadDivisionList]
    PyPlot.plot((loadDivision-0.5) * ones(100),collect(0:99))
end

allTechnologies_plot = Array{TestTechnology,1}(undef,0)
for testGenCo in allTestGenCos
    append!(allTechnologies_plot,testGenCo.technologies)
end
sort!(allTechnologies_plot, by = x -> x.marginalCost)
#allMarginalCost = [technology.marginalCost for technology in allTechnologies_plot]
allTechNames = [string(technology.type) for technology in allTechnologies_plot]
for txt_tech in zip(collect(1:length(averageMarkupHistory_all))[loadDivisionList], allTechNames)
    PyPlot.text(x=txt_tech[1]-4, y=5, s=txt_tech[2])
    PyPlot.text(x=txt_tech[1]-4, y=3, s="(800 MW)")
end
PyPlot.legend(loc = 2)

PyPlot.twinx()
averageElectricityPrice = zeros(length(averageMarkupHistory_all))
for markupCount in 1:length(averageMarkupHistory_all)
    if loadList[markupCount] <= length(allTestGenCos) * 800
        averageElectricityPrice[markupCount] = (averageMarkupHistory_all[markupCount] +1) * baseTech.marginalCost
    elseif loadList[markupCount] <= length(allTestGenCos) * 2 * 800
        averageElectricityPrice[markupCount] = (averageMarkupHistory_all[markupCount] +1) * midTech.marginalCost
    else
        averageElectricityPrice[markupCount] = (averageMarkupHistory_all[markupCount] +1) * peakTech.marginalCost
    end
end
PyPlot.plot(averageElectricityPrice,marker = ".",label = "electricity price [€/MWh]",color = "green")
PyPlot.ylabel("Electricity price",fontsize=16,color = "green")
PyPlot.legend(loc = 1)
#PyPlot.xticks(collect(0:(length(averageMarkupHistory_all)-1)),loadList,rotation = 80)
#=
expectedProfit_Type = zeros(3,length(loadList))
expectedProfit_std = zeros(3,length(loadList))
for allTestGenCosCount in 1:length(allTestGenCos_history)
    thisExpect = zeros(3)
    tempStd = Array{Float64,2}(undef,3,0)
    for testGenCo in allTestGenCos_history[allTestGenCosCount]
        thisExpect = thisExpect + testGenCo.bidSpace[1].expectedProfit_Type
        tempStd = hcat(tempStd,thisExpect)

    end
    println(tempStd)
    for i in 1:3
        expectedProfit_std[i,allTestGenCosCount] = Statistics.std(tempStd[i,:])
    end

    thisExpect = thisExpect ./ length(allTestGenCos)# ./ sum(thisExpect) * 50
    expectedProfit_Type[:,allTestGenCosCount] = thisExpect
end

PyPlot.bar(collect(0:length(loadList)-1),expectedProfit_Type[1,:],0.5,alpha = 0.5,label = "Base")
PyPlot.bar(collect(0:length(loadList)-1),expectedProfit_Type[2,:],0.5,bottom = expectedProfit_Type[1,:],alpha = 0.5,label = "Mid")
PyPlot.bar(collect(0:length(loadList)-1),expectedProfit_Type[3,:],0.5,bottom = expectedProfit_Type[1,:] + expectedProfit_Type[2,:],alpha = 0.5,label = "Peak")
PyPlot.legend()
=#

#=

PyPlot.figure(5)
explorationRateHistory = (1 .- collect(1:totalEpisode)./totalEpisode).^(1/20)
PyPlot.plot(explorationRateHistory)
PyPlot.xlabel("Episode",fontsize=16)
PyPlot.ylabel("Exploration rate",fontsize=16)
PyPlot.xticks(fontsize=16)
PyPlot.yticks(fontsize=16)
#PyPlot.legend(fontsize = 16)
PyPlot.grid(true)

PyPlot.figure(1)
for genCoCount in 1:length(allTestGenCos)
    PyPlot.plot(1:(totalEpisode*simulationTimes),profitHistory_all[genCoCount,:],label = allTestGenCos[genCoCount].name)
    PyPlot.xlabel("Episode",fontsize=16)
    PyPlot.ylabel("Maximum expected profit [€]",fontsize=16)
    PyPlot.xticks(fontsize=16)
    PyPlot.yticks(fontsize=16)
    PyPlot.legend(fontsize = 16)
    PyPlot.title("Max. expected profit vs. simulation")
    PyPlot.grid(true)
end

PyPlot.twinx()
explorationRateHistory = (1 .- collect(1:totalEpisode)./totalEpisode).^(1/20)
PyPlot.plot(repeat(explorationRateHistory,simulationTimes),color = "red")
PyPlot.ylim([0,2])
PyPlot.ylabel("Exploration rate",color = "red")
#PyPlot.plot(1:length(profitHistory_2),profitHistory_3,label = "GenCo3")

PyPlot.figure(2)
for genCoCount in 1:length(allTestGenCos)
    PyPlot.plot(1:(totalEpisode*simulationTimes),marginalMarkupHistory_all[genCoCount,:],label = allTestGenCos[genCoCount].name)
    PyPlot.xlabel("Episode",fontsize=16)
    PyPlot.ylabel("Markup",fontsize=16)
    PyPlot.xticks(fontsize=16)
    PyPlot.yticks(fontsize=16)
    PyPlot.title("Markup percentage of the peak tech. in the optimal bid")
    PyPlot.legend(fontsize = 16)
    PyPlot.grid(true)
end

PyPlot.twinx()
explorationRateHistory = (1 .- collect(1:totalEpisode)./totalEpisode).^(1/20)
PyPlot.plot(repeat(explorationRateHistory,simulationTimes),color = "red")
PyPlot.ylim([0,2])
PyPlot.ylabel("Exploration rate",color = "red")

PyPlot.figure(3)
PyPlot.plot((markupList.-1) .* 100)
PyPlot.xticks(collect(0:(simulationTimes-1)),string.(marginalGeneratorList))
PyPlot.title("Markup percentage of the marginal generator")
PyPlot.xlabel("Marginal technology type",fontsize=16)
PyPlot.ylabel("Markup percentage",fontsize=16)
PyPlot.xticks(fontsize=16)
PyPlot.yticks(fontsize=16)

PyPlot.grid(true)

PyPlot.figure(4)
plotGenCos = deepcopy(allTestGenCos)
for testGenCo in plotGenCos
    testGenCo.SubmitBid(testGenCo.bidSpace[1])
end
PlotTrade(plotGenCos,loadList[1])
=#


text:
module FlutterBird

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Flutterbird" => Ahorn.EntityPlacement(
        Maple.Flutterbird
    ),
)

sprite = "scenery/flutterbird/idle00.png"

colors = [
    (137, 251, 255, 255) ./ 255,
    (240, 252, 108, 255) ./ 255,
    (244, 147, 255, 255) ./ 255,
    (147, 186, 255, 255) ./ 255
]

function Ahorn.selection(entity::Maple.Flutterbird)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y, jx=0.5, jy=1.0)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Flutterbird, room::Maple.Room)
    rng = Ahorn.getSimpleEntityRng(entity)
    color = rand(rng, colors)

    Ahorn.drawSprite(ctx, sprite, 0, 0, jx=0.5, jy=1.0, tint=color)
end

end


text:
function display_map(map_matrix::Matrix{Char}; path=nothing)
    height, width = size(map_matrix)
    map_colors = Matrix{RGB}(undef, height, width)
    for i in 1:height, j in 1:width
        c = map_matrix[i, j]
        if c == '.'
            x = RGB(1.0, 1.0, 1.0)
        elseif c == 'G'
            x = RGB(1.0, 1.0, 1.0)
        elseif c == 'S'
            x = RGB(0.4, 0.4, 0.0)
        elseif c == 'W'
            x = RGB(0.0, 0.0, 1.0)
        elseif c == 'T'
            x = RGB(0.0, 0.7, 0.0)
        elseif c == '@'
            x = RGB(0.0, 0.0, 0.0)
        elseif c == 'O'
            x = RGB(0.0, 0.0, 0.0)
        else
            x = RGB(0.0, 0.0, 0.0)
        end
        map_colors[i, j] = x
    end
    if !isnothing(path)
        g = GridGraph(map_matrix)
        for ed in path
            s, d = src(ed), dst(ed)
            (is, js) = node_coord(g, s)
            (id, jd) = node_coord(g, d)
            map_colors[is, js] = RGB(1., 0., 0.)
            map_colors[id, jd] = RGB(1., 0., 0.)
        end
    end
    return map_colors
end


text:
module GoogleProtoBuf
    include("descriptor_pb.jl")
end
module GoogleProtoBufCompiler
    include("plugin_pb.jl")
end


text:
function motor_ramp_test(uri=_first_available())
    play(uri) do cf
        cf.commander.send_setpoint(0, 0, 0, 0)
        thrust = 20000
        for i = 1:20
            cf.commander.send_setpoint(0, 0, 0, thrust)
            thrust += i > 10 ? -500 : 500
            sleep(0.1)
        end
        cf.commander.send_setpoint(0, 0, 0, 0)
        sleep(0.1)
    end
end


text:
function preview_coverage(; dump_coverage_io::IO = stdout,
                            min_padding::Integer = 1,
                            max_padding::Integer = 24)
    setup()
    all_tracked_files = tracked_files()
    for x in all_tracked_files
        backup_filename = string(x, ".backup")
        coverage_filename = string(x, ".microcov")
        coverage_lines = _compute_coverage_file_contents(x;
                                                         min_padding = min_padding,
                                                         max_padding = max_padding,
                                                         source_code_filename = Symbol(backup_filename))
        # println(dump_coverage_io, "\n")
        println(dump_coverage_io, "# Preview of $(coverage_filename):")
        for x in coverage_lines
            println(dump_coverage_io, x)
        end
        # println(dump_coverage_io, "\n")
    end
    return nothing
end


text:
using Plots
import Plots.plot

"""
Plot a volume image

#### Arguments

* `v`: A VolumeImage type
* `threshold`: Minimum value to plot, values smaller than this are plotted as `minsize`
* `min_val`: Force a minimum value for color and size scaling
* `max_val`: Force a maximum value for color and size scaling
* `minsize`: Minimum size a marker can be
* `maxsize`: Maximum size a marker can be
* `exclude`: Values not to plot
* `title`: Figure title
* `elp`: Path to elp file to overlay channel names
* `colorbar`: Should a colorbar be plotted

#### Returns

* A Plots.jl figure

"""
function plot(v::VolumeImage; kwargs...)

    x = AbstractFloat[xi / (1 * u"mm") for xi in v.x]
    y = AbstractFloat[yi / (1 * u"mm") for yi in v.y]
    z = AbstractFloat[zi / (1 * u"mm") for zi in v.z]

    plot_src(dropdims(v.data, dims = 4), x, y, z; kwargs...)
end

function plot_src(
    d::Array{A,3},
    x::Vector{A},
    y::Vector{A},
    z::Vector{A};
    threshold::Real = -Inf,
    min_val::Real = Inf,
    max_val::Real = -Inf,
    minsize::Real = 2,
    maxsize::Real = 6,
    exclude::A = 0.0,
    title::S = "",
    elp::AbstractString = "",
    colorbar::Bool = true,
    kwargs...,
) where {A<:AbstractFloat,S<:AbstractString}

    # cols = [colorant"darkblue", colorant"orange", colorant"darkred"]

    scaleval = maxsize / maximum(d)

    plot_labels = false
    if elp != ""
        e = read_elp(elp)
        plot_labels = true
    end

    #
    # First facet

    x_tmp = AbstractFloat[]
    y_tmp = AbstractFloat[]
    c_tmp = AbstractFloat[]
    s_tmp = AbstractFloat[]
    t = copy(d)
    t = maximum(t, dims = 3)
    t = dropdims(t, dims = 3)
    for x_i = 1:size(t, 1)
        for y_i = 1:size(t, 2)
            val = t[x_i, y_i]
            if val != exclude
                push!(x_tmp, x[x_i])
                push!(y_tmp, y[y_i])
                if val > threshold
                    push!(s_tmp, max(scaleval * val, minsize))
                    push!(c_tmp, val)
                else
                    push!(s_tmp, minsize)
                    push!(c_tmp, minsize)
                end
            end
        end
    end
    if max_val > maximum(c_tmp)
        @debug("Manually specifying maximum plotting value")
        push!(s_tmp, max_val)
        push!(c_tmp, max_val)
        push!(x_tmp, -200)
        push!(y_tmp, -200)
    end
    if min_val < minimum(c_tmp)
        @debug("Manually specifying minimum plotting value")
        push!(s_tmp, min_val)
        push!(c_tmp, min_val)
        push!(x_tmp, -200)
        push!(y_tmp, -200)
    end
    p1 = plot(
        x_tmp,
        y_tmp,
        zcolor = c_tmp,
        ms = s_tmp,
        legend = false,
        l = :scatter,
        lab = "Source",
        colorbar = false,
        markerstrokewidth = 0.1,
        xlabel = "Left - Right (mm)",
        ylabel = "Posterior - Anterior (mm)",
        xlims = (-100, 100),
        ylims = (-120, 90);
        kwargs...,
    )
    if plot_labels
        plotlist = [
            "Fpz",
            "Fp2",
            "AF8",
            "F8",
            "FT8",
            "T8",
            "TP8",
            "P10",
            "PO8",
            "O2",
            "Oz",
            "O1",
            "PO7",
            "P9",
            "TP7",
            "T7",
            "FT7",
            "F7",
            "AF7",
            "Fp1",
        ]
        for elec in e
            if something(findfirst(isequal(elec.label), plotlist), 0) > 0
                annotate!(
                    p1,
                    elec.coordinate.x - 5,
                    1.1 * elec.coordinate.y - 2,
                    elec.label,
                    colorbar = false,
                )
            end
        end
    end

    #
    # Second facet

    x_tmp = AbstractFloat[]
    y_tmp = AbstractFloat[]
    c_tmp = AbstractFloat[]
    s_tmp = AbstractFloat[]
    t = copy(d)
    t = maximum(t, dims = 1)
    t = dropdims(t, dims = 1)
    for x_i = 1:size(t, 1)
        for y_i = 1:size(t, 2)
            val = t[x_i, y_i]
            if val != exclude
                push!(x_tmp, y[x_i])
                push!(y_tmp, z[y_i])
                if val > threshold
                    push!(s_tmp, max(scaleval * val, minsize))
                    push!(c_tmp, val)
                else
                    push!(s_tmp, minsize)
                    push!(c_tmp, minsize)
                end
            end
        end
    end
    if max_val > maximum(c_tmp)
        @debug("Manually specifying maximum plotting value")
        push!(s_tmp, max_val)
        push!(c_tmp, max_val)
        push!(x_tmp, -200)
        push!(y_tmp, -200)
    end
    if min_val < minimum(c_tmp)
        @debug("Manually specifying minimum plotting value")
        push!(s_tmp, min_val)
        push!(c_tmp, min_val)
        push!(x_tmp, -200)
        push!(y_tmp, -200)
    end
    p2 = plot(
        x_tmp,
        y_tmp,
        zcolor = c_tmp,
        ms = s_tmp,
        legend = false,
        l = :scatter,
        title = title,
        lab = "Source",
        colorbar = false,
        markerstrokewidth = 0.1,
        xlabel = "Posterior - Anterior (mm)",
        ylabel = "Inferior - Superior (mm)",
        xlims = (-120, 90),
        ylims = (-70, 100);
        kwargs...,
    )
    if plot_labels
        plotlist = ["Iz", "Oz", "POz", "Pz", "CPz", "Cz", "FCz", "Fz", "AFz", "Fpz"]
        for elec in e
            if something(findfirst(isequal(elec.label), plotlist), 0) > 0
                annotate!(p2, elec.coordinate.y - 5, elec.coordinate.z, elec.label)
            end
        end
    end

    #
    # Third facet

    x_tmp = AbstractFloat[]
    y_tmp = AbstractFloat[]
    c_tmp = AbstractFloat[]
    s_tmp = AbstractFloat[]
    t = copy(d)
    t = maximum(t, dims = 2)
    t = dropdims(t, dims = 2)
    for x_i = 1:size(t, 1)
        for y_i = 1:size(t, 2)
            val = t[x_i, y_i]
            if val != exclude
                push!(x_tmp, x[x_i])
                push!(y_tmp, z[y_i])
                if val > threshold
                    push!(s_tmp, max(scaleval * val, minsize))
                    push!(c_tmp, val)
                else
                    push!(s_tmp, minsize)
                    push!(c_tmp, minsize)
                end
            end
        end
    end
    if max_val > maximum(c_tmp)
        @debug("Manually specifying maximum plotting value")
        push!(s_tmp, max_val)
        push!(c_tmp, max_val)
        push!(x_tmp, -200)
        push!(y_tmp, -200)
    end
    if min_val < minimum(c_tmp)
        @debug("Manually specifying minimum plotting value")
        push!(s_tmp, min_val)
        push!(c_tmp, min_val)
        push!(x_tmp, -200)
        push!(y_tmp, -200)
    end
    p3 = plot(
        x_tmp,
        y_tmp,
        zcolor = c_tmp,
        ms = s_tmp,
        legend = false,
        l = :scatter,
        lab = "",
        markerstrokewidth = 0.1,
        colorbar = colorbar,
        xlabel = "Left - Right (mm)",
        ylabel = "Inferior - Superior (mm)",
        xlims = (-100, 100),
        ylims = (-70, 100);
        kwargs...,
    )
    if plot_labels
        plotlist = ["T7", "C5", "C3", "C1", "Cz", "C2", "C4", "C6", "T8"]
        for elec in e
            if something(findfirst(isequal(elec.label), plotlist), 0) > 0
                annotate!(p3, elec.coordinate.x - 5, elec.coordinate.z, elec.label)
            end
        end
    end

    l = @layout([a b c])
    return plot(p1, p2, p3, layout = l)
end


text:
####
#### COLOR:
####  set_color!, set_colors!, set_fill!, set_fills!, set_alpha!
####

"""
    set_color!(obj, col)

Internal functions to set the color value `col` (after parsing) to the appropriate
field of object `obj`.
"""
set_color!(o::Union{TextStyle,LineStyle,MarkerStyle,BarStyle}, c::Color) = (o.color = c)

set_color!(o::Union{Figure,Legend}, c::Option{Color}) = (o.bgcolor = c)

set_color!(o::Hist2D, c::Color) = set_color!(o.barstyle, c)
set_color!(o::Union{Ticks,StraightLine2D,Box2D,Scatter3D}, c::Color) = set_color!(o.linestyle, c)
set_color!(o::Colorbar, c::Color) = set_color!(o.ticks, c)

set_textcolor!(o::Ticks, c::Color) = set_color!(o.labels.textstyle, c)
set_textcolor!(o::Colorbar, c::Color) = set_color!(o.ticks, c)
set_textcolor!(o::Union{Figure,Axis,Title}, c::Color) = set_color!(o.textstyle, c)

"""
    set_fill!(obj, col)

Internal functions to set the fill color value `v` (after parsing) to the appropriate
field of object `obj`.
"""
set_fill!(o::Union{FillStyle,BarStyle}, c::Color) = (o.fill = c)
set_fill!(o::Fill2D, c::Color) = set_fill!(o.fillstyle, c)
set_fill!(o::Hist2D, c::Color) = set_fill!(o.barstyle, c)

function set_fill!(o::Box2D, c::Color)
    isnothing(o.fillstyle) && (o.fillstyle = FillStyle())
    set_fill!(o.fillstyle, c)
end

"""
    set_colors!(obj, cols)

Internal function to set the color values `cols`.
If a single value is passed, all fields will be assigned to that value.
"""
function set_colors!(vs::Vector, c::Union{Color, Vector{<:Color}}, field::Option{Symbol}=nothing)
    c isa Vector || (c = fill(c, length(vs)))
    # check dimensions match
    length(c) == length(vs) || throw(DimensionMismatch("colors // dimensions don't match"))
    # assign
    if isnothing(field)
        for (i, col) ∈ enumerate(c)
            set_color!(vs[i], col)
        end
    else
        for (i, col) ∈ enumerate(c)
            set_color!(getfield(vs[i], field), col)
        end
    end
    return nothing
end

"""
    set_fills!(obj, cols)

Internal functions to set the fill color values `cols`.
If a single value is passed, all fields will be assigned to that value.
"""
function set_fills!(vs::Vector, c::Union{Color, Vector{<:Color}}, field::Option{Symbol}=nothing)
    c isa Vector || (c = fill(c, length(vs)))
    # check dimensions match
    length(c) == length(vs) || throw(DimensionMismatch("fills // dimensions don't match"))
    # assign
    if isnothing(field)
        for (i, col) ∈ enumerate(c)
            set_fill!(vs[i], col)
        end
    else
        for (i, col) ∈ enumerate(c)
            set_fill!(getfield(vs[i], field), col)
        end
    end
    return nothing
end

"""
    set_alpha!(obj, α)

Internal function to set the alpha value of `obj.field` to `α`. There must be a color
value available, it will be reinterpreted with the given alpha value.
"""
function set_alpha!(o::Union{Fill2D,Hist2D,Box2D}, α::Float64, parent::Symbol)
    eval(:($o.$parent.fill = coloralpha($o.$parent.fill, $α)))
    return nothing
end
set_alpha!(o::Union{Fill2D,Box2D}, α::Float64) = set_alpha!(o, α, :fillstyle)
set_alpha!(o::Hist2D, α::Float64) = set_alpha!(o, α, :barstyle)
set_alpha!(o::Union{Figure,Legend}, α::Float64) = (o.bgcolor = coloralpha(o.bgcolor, α); ∅)

####
#### TEXT
####  set_font!, set_hei!
####

"""
    set_font!(obj, font)

Internal function to set the font associated with an object `obj` to the value `font`.
"""
function set_font!(ts::TextStyle, font::String)
    @assert get_backend() == GLE "font // only GLE backend supported"
    ts.font = get(GLE_FONTS, font) do
        throw(OptionValueError("font", font))
    end
    return nothing
end
set_font!(o::Ticks, f::String)    = set_font!(o.labels.textstyle, f)
set_font!(o::Colorbar, f::String) = set_font!(o.ticks, f)
set_font!(o, f::String)           = set_font!(o.textstyle, f)

"""
    set_hei!(obj, fontsize)

Internal function to set the font associated with an object `obj` to the value `font`.
"""
set_hei!(ts::TextStyle, v::Float64) = (ts.hei = v * PT_TO_CM)
set_hei!(o::Ticks, v::Float64)    = set_hei!(o.labels.textstyle, v)
set_hei!(o::Colorbar, v::Float64) = set_hei!(o.ticks, v)
set_hei!(o, v::Float64)           = set_hei!(o.textstyle, v)

####
#### Line related
####

"""
    set_lstyle!(obj, lstyle)

Internal function to set the line style associated with object `obj`. The style
can be described by `lstyle` being a number or a String representing the pattern.
"""
function set_lstyle!(o::LineStyle, v::Int)
    0 < v || throw(OptionValueError("lstyle", v))
    o.lstyle = v
    return nothing
end
function set_lstyle!(o::LineStyle, v::String)
    @assert get_backend() == GLE "lstyle // only GLE backend supported"
    o.lstyle = get(GLE_LSTYLES, v) do
        throw(OptionValueError("lstyle", v))
    end
    return nothing
end
set_lstyle!(o, v::String) = set_lstyle!(o.linestyle, v)

"""
    set_lwidth!(obj, v)

Internal function to set the line width associated with the relevant field of `obj`.
"""
set_lwidth!(o::Union{LineStyle, Axis}, v::Float64) = (o.lwidth = v)
set_lwidth!(o, v::Float64) = set_lwidth!(o.linestyle, v)

"""
    set_smooth!(obj, v)

Internal function to determine whether to use splines for a field of `obj`.
"""
set_smooth!(o::LineStyle, v::Bool) = (o.smooth = v)
set_smooth!(o::Scatter3D, v) = throw(NotImplementedError("smooth for 3d line"))

####
#### Marker related
####

"""
    set_marker!(obj, marker)

Internal function to set the marker associated with object `obj`. The style
can be described by `marker` being a String describing the pattern.
"""
function set_marker!(o::MarkerStyle, v::String)
    @assert get_backend() == GLE "marker // only GLE backend supported"
    o.marker = get(GLE_MARKERS, v) do
        throw(OptionValueError("marker", v))
    end
    return nothing
end
set_marker!(o, v::String) = set_marker!(o.markerstyle, v)

"""
    set_msize!(obj, msize)

Internal function to set the marker size associated with object `obj`. If no marker has been given,
set it to a filled circle.
"""
function set_msize!(o::MarkerStyle, v::Float64)
    isdef(o.marker) || (o.marker = "fcircle")
    o.msize = v
    return nothing
end
set_msize!(o, v::Float64) = set_msize!(o.markerstyle, v)

"""
    set_mcol!(obj, col)

Internal function to set the marker color. If no marker has been given, set it to a filled circle.
"""
function set_mcol!(o::MarkerStyle, c::Color)
    isdef(o.marker) || (o.marker = "fcircle")
    o.color = c
    return nothing
end
set_mcol!(o, v::Color) = set_mcol!(o.markerstyle, v)

# generate functions that take vector inputs for linestyle and markerstyle
for case ∈ ("lstyle", "lwidth", "smooth", "marker", "msize", "mcol")
    f_scalar! = Symbol("set_$(case)!")  # function with scalar input
    f_vector! = Symbol("set_$(case)s!") # e.g. set_markers!
    ex = quote
        # set function for a group of objects (e.g. linestyles, markerstyles)
        function $f_vector!(vs::Vector, v, f::Option{Symbol}=∅)
            v isa Vector || (v = fill(v, length(vs)))
            length(vs) == length(v) || throw(DimensionMismatch($case*"s // dimensions don't match"))
            if !isdef(f)
                for (i, vi) ∈ enumerate(v)
                    $f_scalar!(vs[i], vi)
                end
            else
                for (i, vi) ∈ enumerate(v)
                    $f_scalar!(getfield(vs[i], f), vi)
                end
            end
            return nothing
        end
    end
    eval(ex)
end

####
#### Bar related
####

"""
    set_bwidth!(obj, v)

Internal function to set the bin or box width to value `v`.
"""
set_bwidth!(o::Bar2D, v::Float64) = (o.bwidth = v)

####
#### Boxplot related
####

function set_bp!(b::Boxplot, field::Symbol, v)
    v isa Vector || (v = fill(v, b.nobj))
    length(v) == b.nobj || throw(DimensionMismatch("$(field)s // dimensions don't match"))
    for k ∈ 1:b.nobj
        setfield!(b.boxstyles[k], field, v[k])
    end
    return nothing
end

"""
    set_bwidths!(b, v)

Internal function to set the box widths to values `v`.
"""
set_bwidths!(b::Boxplot, v) = set_bp!(b, :bwidth, v)

"""
    set_wwidth!(b, v)

Internal function to set the whisker width to value `v`.
"""
set_wwidths!(b::Boxplot, v) = set_bp!(b, :wwidth, v)

"""
    set_wrlengths!(b, v)

Internal function to set the whiskers relative lengths to value `v`.
"""
set_wrlengths!(b::Boxplot, v) = set_bp!(b, :wrlength, v)

"""
    set_show!(b, v, f)

Internal function to toggle the display of an element `f` in a boxplot. See [`set_mshow!`](@ref)
and [`set_oshow!`](@ref).
"""
function set_show!(b::Boxplot, v::Union{Bool, Vector{Bool}}, field::Symbol)
    v isa Vector || (v = fill(v, b.nobj))
    length(v) == b.nobj || throw(DimensionMismatch("set_mshow // dimensions don't match"))
    for (k, vk) ∈ enumerate(v)
        setfield!(b.boxstyles[k], field, vk)
    end
    return nothing
end

"""
    set_mshow!(b, v)

Internal function to toggle the display of the mean point in a boxplot. (A vector can
be passed to set this individually on each boxplot shown).
"""
set_mshow!(b::Boxplot, v) = set_show!(b, v, :mshow)

"""
    set_oshow!(b, v)

Internal function to toggle the display of outliers in a boxplot. (A vector can
be passed to set this individually on each boxplot shown).
"""
set_oshow!(b::Boxplot, v) = set_show!(b, v, :oshow)


text:
#
# Example usage of cmocean.jl
#

using PyPlot
include("cmocean.jl")

data = ones(10) * linspace(0,1,101)'

figure(1)
clf()
subplot(3,1,1)
imshow(data, cmap="ice")

subplot(3,1,2)
imshow(data, cmap="speed")

subplot(3,1,3)
imshow(data, cmap="balance")


text:
# ---
# title: 621. Task Scheduler
# id: problem621
# author: AquaIndigo
# date: 2020-12-03
# difficulty: Medium
# categories: Array, Greedy, Queue
# link: <https://leetcode.com/problems/task-scheduler/description/>
# hidden: true
# ---
#
# Given a characters array `tasks`, representing the tasks a CPU needs to do,
# where each letter represents a different task. Tasks could be done in any
# order. Each task is done in one unit of time. For each unit of time, the CPU
# could complete either one task or just be idle.
#
# However, there is a non-negative integer `n` that represents the cooldown
# period between two **same tasks**  (the same letter in the array), that is
# that there must be at least `n` units of time between any two same tasks.
#
# Return _the least number of units of times that the CPU will take to finish
# all the given tasks_.
#
#
#
# **Example 1:**
#
#
#
#     Input: tasks = ["A","A","A","B","B","B"], n = 2
#     Output: 8
#     Explanation:
#     A -> B -> idle -> A -> B -> idle -> A -> B
#     There is at least 2 units of time between any two same tasks.
#
#
# **Example 2:**
#
#
#
#     Input: tasks = ["A","A","A","B","B","B"], n = 0
#     Output: 6
#     Explanation: On this case any permutation of size 6 would work since n = 0.
#     ["A","A","A","B","B","B"]
#     ["A","B","A","B","A","B"]
#     ["B","B","B","A","A","A"]
#     ...
#     And so on.
#
#
# **Example 3:**
#
#
#
#     Input: tasks = ["A","A","A","A","A","A","B","C","D","E","F","G"], n = 2
#     Output: 16
#     Explanation:
#     One possible solution is
#     A -> B -> C -> A -> D -> E -> A -> F -> G -> A -> idle -> idle -> A -> idle -> idle -> A
#
#
#
#
# **Constraints:**
#
#   * `1 <= task.length <= 104`
#   * `tasks[i]` is upper-case English letter.
#   * The integer `n` is in the range `[0, 100]`.
#
#
## @lc code=start
using LeetCode

function least_interval(tasks::Vector{Char}, n::Int)::Int
    counter = fill(0, 26)
    for ch in tasks
        counter[UInt8(ch) - 0x40] += 1
    end
    sort!(counter, rev = true)
    res = (n + 1) * counter[1] - n
    for i in 2:26
        if counter[i] == counter[i - 1]
            res += 1
        else
            break
        end
    end
    return max(length(tasks), res)
end
least_interval(['A','A','A','B','B','B'], 2)
## @lc code=end


text:
using Base.Test
using Compat: readstring
using LegacyStrings: UTF8String, UTF16String, UTF32String
using StringEncodings

for s in ("", "\0", "a", "café crème",
          "a"^(StringEncodings.BUFSIZE-1) * "€ with an incomplete codepoint between two input buffer fills",
          "a string € チャネルパートナーの選択",
          "a string \0€ チャネルパ\0ー\0トナーの選択 with embedded and trailing nuls\0")
    # Test round-trip to Unicode formats, checking against pure-Julia implementation
    for (T, nullen) in ((UTF8String, 0), (UTF16String, 2), (UTF32String, 4))
        enc = StringEncodings.encoding(T)
        a = reinterpret(UInt8, T(s).data)
        # Adjust for explicit \0 only for .data on UTF16String/UTF32String
        a = a[1:end - nullen]
        @test decode(a, enc) == s
        @test decode(UTF16String, a, enc) == s
        @test decode(UTF32String, a, enc) == s
        @test decode(encode(s, enc), enc) == s
    end
end

# Test a few non-Unicode encodings
for (s, enc) in (("noël", "ISO-8859-1"),
                 ("noël €", "ISO-8859-15", "CP1252"),
                 ("Код Обмена Информацией, 8 бит", "KOI8-R"),
                 ("国家标准", "GB18030"))
    @test decode(encode(s, enc), enc) == s
end

# Test that attempt to close stream in the middle of incomplete sequence throws
let s = "a string チャネルパートナーの選択"
    # First, correct version
    p = StringEncoder(IOBuffer(), "UTF-16LE")
    write(p, s.data)
    close(p)
    # Test that closed pipe behaves correctly
    @test_throws ArgumentError write(p, 'a')

    b = IOBuffer()
    p = StringEncoder(b, "UTF-16LE")
    @test string(p) == "StringEncoder{UTF-8, UTF-16LE}($(string(b)))"
    write(p, s.data[1:10])
    @test_throws IncompleteSequenceError close(p)
    # Test that closed pipe behaves correctly even after an error
    @test_throws ArgumentError write(p, 'a')

    # This time, call show
    p = StringEncoder(IOBuffer(), "UTF-16LE")
    write(p, s.data[1:10])
    try
        close(p)
    catch err
        @test isa(err, IncompleteSequenceError)
        io = IOBuffer()
        showerror(io, err)
        @test takebuf_string(io) == "Incomplete byte sequence at end of input"
    end

    b = IOBuffer(encode(s, "UTF-16LE")[1:19])
    p = StringDecoder(b, "UTF-16LE")
    @test string(p) == "StringDecoder{UTF-16LE, UTF-8}($(string(b)))"
    @test readstring(p) == s[1:9]
    @test_throws IncompleteSequenceError close(p)
    # Test that closed pipe behaves correctly even after an error
    @test eof(p)
    @test_throws ArgumentError read(p, UInt8)

    # Test stateful encoding, which output some bytes on final reset
    # with strings containing different scripts
    x = encode(s, "ISO-2022-JP")
    @test decode(x, "ISO-2022-JP") == s

    p = StringDecoder(b, "ISO-2022-JP", "UTF-8")
    b = IOBuffer(x)
    # Test that closed pipe behaves correctly
    close(p)
    @test eof(p)
    @test_throws ArgumentError read(p, UInt8)
    close(p)
end

@test_throws InvalidSequenceError encode("qwertyé€", "ASCII")
try
    encode("qwertyé€", "ASCII")
catch err
     io = IOBuffer()
     showerror(io, err)
     @test takebuf_string(io) ==
        "Byte sequence 0xc3a9e282ac is invalid in source encoding or cannot be represented in target encoding"
end

@test_throws InvalidSequenceError decode(b"qwertyé€", "ASCII")
try
    decode(b"qwertyé€", "ASCII")
catch err
     io = IOBuffer()
     showerror(io, err)
     @test takebuf_string(io) ==
         "Byte sequence 0xc3a9e282ac is invalid in source encoding or cannot be represented in target encoding"
end

let x = encode("ÄÆä", "ISO-8859-1")
    @test_throws InvalidSequenceError decode(x, "UTF-8")
    try
        decode(x, "UTF-8")
    catch err
         io = IOBuffer()
         showerror(io, err)
         @test takebuf_string(io) ==
             "Byte sequence 0xc4c6e4 is invalid in source encoding or cannot be represented in target encoding"
    end
end

mktemp() do p, io
    s = "café crème"
    write(io, encode(s, "CP1252"))
    close(io)
    @test readstring(p, enc"CP1252") == s
end

@test_throws InvalidEncodingError p = StringEncoder(IOBuffer(), "nonexistent_encoding")
@test_throws InvalidEncodingError p = StringDecoder(IOBuffer(), "nonexistent_encoding")

try
    p = StringEncoder(IOBuffer(), "nonexistent_encoding", "absurd_encoding")
catch err
    @test isa(err, InvalidEncodingError)
    io = IOBuffer()
    showerror(io, err)
    @test takebuf_string(io) ==
        "Conversion from absurd_encoding to nonexistent_encoding not supported by iconv implementation, check that specified encodings are correct"
end
try
    p = StringDecoder(IOBuffer(), "nonexistent_encoding", "absurd_encoding")
catch err
    @test isa(err, InvalidEncodingError)
    io = IOBuffer()
    showerror(io, err)
    @test takebuf_string(io) ==
        "Conversion from nonexistent_encoding to absurd_encoding not supported by iconv implementation, check that specified encodings are correct"
end

mktemp() do path, io
    s = "a string \0チャネルパ\0ー\0トナーの選択 with embedded and trailing nuls\0\nand a second line"
    close(io)
    open(path, enc"ISO-2022-JP", "w") do io
        @test iswritable(io) && !isreadable(io)
        write(io, s)
    end

    @test readstring(path, enc"ISO-2022-JP") == s
    @test open(io->readstring(io, enc"ISO-2022-JP"), path) == s
    @test open(readstring, path, enc"ISO-2022-JP") == s

    @test readuntil(path, enc"ISO-2022-JP", '\0') == "a string \0"
    @test open(io->readuntil(io, enc"ISO-2022-JP", '\0'), path) == "a string \0"
    @test readuntil(path, enc"ISO-2022-JP", "チャ") == "a string \0チャ"
    @test open(io->readuntil(io, enc"ISO-2022-JP", "チャ"), path) == "a string \0チャ"

    @test readline(path, enc"ISO-2022-JP") == string(split(s, '\n')[1], '\n')
    @test open(readline, path, enc"ISO-2022-JP") == string(split(s, '\n')[1], '\n')
    a = readlines(path, enc"ISO-2022-JP")
    b = open(readlines, path, enc"ISO-2022-JP")

    c = collect(eachline(path, enc"ISO-2022-JP"))
    d = open(io->collect(eachline(io, enc"ISO-2022-JP")), path)

    @test a[1] == b[1] == c[1] == d[1] == string(split(s, '\n')[1], '\n')
    @test a[2] == b[2] == c[2] == d[2] == split(s, '\n')[2]

    # Test alternative syntaxes for open()
    open(path, enc"ISO-2022-JP", "r") do io
        @test isreadable(io) && !iswritable(io)
        @test readstring(io) == s
    end
    open(path, enc"ISO-2022-JP", true, false, false, false, false) do io
        @test isreadable(io) && !iswritable(io)
        @test readstring(io) == s
    end
    @test_throws ArgumentError open(path, enc"ISO-2022-JP", "r+")
    @test_throws ArgumentError open(path, enc"ISO-2022-JP", "w+")
    @test_throws ArgumentError open(path, enc"ISO-2022-JP", "a+")
    @test_throws ArgumentError open(path, enc"ISO-2022-JP", true, true, false, false, false)
    @test_throws ArgumentError open(path, enc"ISO-2022-JP", true, false, false, false, true)
end

## Test encodings support
b = IOBuffer()
show(b, enc"UTF-8")
@test takebuf_string(b) == "UTF-8 string encoding"
@test string(enc"UTF-8") == "UTF-8"

encodings_list = encodings()
@test "ASCII" in encodings_list
@test "UTF-8" in encodings_list

nothing


text:
##################################################
## evaluate binary operations of symbolic objects
## XXX -- this may prove too narryw with the use of ::Sym

+(x::SymbolicObject, y::SymbolicObject) = x.__add__(y)
*(x::SymbolicObject, y::SymbolicObject) = x.__mul__(y)
-(x::SymbolicObject, y::SymbolicObject) = x.__sub__(y)
-(x::SymbolicObject)                    = x.__neg__()
/(x::SymbolicObject, y::SymbolicObject) = x.__truediv__(y)
^(x::SymbolicObject, y::SymbolicObject) = x.__pow__(y)
^(x::SymbolicObject, y::Rational) = x^convert(Sym,y)
#^(x::SymbolicObject, y::Integer) = x^convert(Sym,y) # no Union{Integer, Rational}, as that has ambiguity
//(x::SymbolicObject, y::Int) = x / Sym(y)
//(x::SymbolicObject, y::Rational) = x / Sym(y)
//(x::SymbolicObject, y::SymbolicObject) = x / y

\(x::SymbolicObject, y::SymbolicObject) = (y'/x')' # ?

inv(x::Sym) = x.__pow__(Sym(-1))

# special case Boolean; issue   351
# promotion for Boolean here is to 0 or  1,  not False,  True
+(x::Bool, y::SymbolicObject) = Sym(Int(x)).__add__(y)
*(x::Bool, y::SymbolicObject) = Sym(Int(x)).__mul__(y)
-(x::Bool, y::SymbolicObject) = Sym(Int(x)).__sub__(y)
/(x::Bool, y::SymbolicObject) = Sym(Int(x)).__truediv__(y)
^(x::Bool, y::SymbolicObject) = Sym(Int(x)).__pow__(y)
+(x::SymbolicObject, y::Bool) = x.__add__(Int(y))
*(x::SymbolicObject, y::Bool) = x.__mul__(Int(y))
-(x::SymbolicObject, y::Bool) = x.__sub__(Int(y))
/(x::SymbolicObject, y::Bool) = x.__truediv__(Int(y))
^(x::SymbolicObject, y::Bool) = x.__pow__(Int(y))


text:
#=
by David Sanders
=#

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
using NCDatasets
using Formatting
using ArgParse
using Statistics
using JSON

include("constants.jl")

function mreplace(x)
    return replace(x, missing=>NaN)
end

function integrate(x, dydx)
    y = copy(dydx) * 0.0

    for i = 2:length(x)
        y[i] = y[i-1] + (dydx[i-1] + dydx[i]) * (x[i] - x[i-1]) / 2.0
    end

    return y
end

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--data-file"
            help = "Ocean data file."
            arg_type = String
            required = true

        "--output-file"
            help = "Output file."
            arg_type = String
            required = true

        "--domain-file"
            help = "Domain file."
            arg_type = String
            required = true

        "--beg-year"
            help = "Year of begin."
            arg_type = Int64
            required = true

        "--end-year"
            help = "Year of end."
            arg_type = Int64
            required = true

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

Dataset(parsed["domain-file"], "r") do ds

    global mask = ds["mask"][:] |> mreplace
    global lat = ds["yc"][1, :] |> mreplace

    global lat_weight = 2π * Re^2.0 * cos.(deg2rad.(lat))
    global ϕ          = deg2rad.(lat)

#    global O_idx = (mask .== 0)

end

Dataset(parsed["data-file"], "r") do ds

    beg_t = (parsed["beg-year"] - 1) * 12 + 1
    end_t = (parsed["end-year"] - 1) * 12 + 12
    rng = (:,:,beg_t:end_t)

    getData = (varname) -> mean( replace(ds[varname][rng...], missing=>0.0, NaN=>0.0), dims=(1, ))[1, :, :]

    global neb = getData("neb")
    global Q_clim = getData("Q_clim")
    global qflx     = getData("qflx")
    global wT       = getData("wT") * ρc
    global nswflx   = getData("nswflx")
    global swflx    = getData("swflx")
    global dHdt    = getData("dHdt")

    global qflx2atm = replace(ds["qflx2atm"][rng...], missing=>0.0, NaN=>0.0)
    qflx2atm[qflx2atm .< 0.0] .= 0.0
    qflx2atm = mean( qflx2atm, dims=(1, ))[1, :, :]

    global (Ny, Nt) = size(neb)

end

IET_neb       = zeros(Float64, Ny, Nt)
IET_Q_clim    = zeros(Float64, Ny, Nt)
IET_qflx      = zeros(Float64, Ny, Nt)
IET_wT        = zeros(Float64, Ny, Nt)
IET_nswflx    = zeros(Float64, Ny, Nt)
IET_swflx     = zeros(Float64, Ny, Nt)
IET_dHdt      = zeros(Float64, Ny, Nt)

for t = 1:Nt
    IET_neb[:, t]    = integrate(ϕ, neb[:, t] .* lat_weight)
    IET_Q_clim[:, t] = integrate(ϕ, Q_clim[:, t] .* lat_weight)
    IET_qflx[:, t]   = integrate(ϕ, qflx[:, t] .* lat_weight)
    IET_wT[:, t]     = integrate(ϕ, wT[:, t] .* lat_weight)
    IET_nswflx[:, t] = integrate(ϕ, nswflx[:, t] .* lat_weight)
    IET_swflx[:, t]  = integrate(ϕ, swflx[:, t] .* lat_weight)
    IET_dHdt[:, t]   = integrate(ϕ, dHdt[:, t] .* lat_weight)
end

Dataset(parsed["output-file"], "c") do ds

    defDim(ds, "Nx",  1)
    defDim(ds, "Ny", Ny)
    defDim(ds, "Nz",  1)
    defDim(ds, "time", Inf)

    for (varname, vardata, vardim, attrib) in [
        ("IET_neb",        reshape(IET_neb,    1, Ny, 1, Nt), ("Nx", "Ny", "Nz", "time"), Dict()),
        ("IET_Q_clim",     reshape(IET_Q_clim, 1, Ny, 1, Nt), ("Nx", "Ny", "Nz", "time"), Dict()),
        ("IET_qflx",       reshape(IET_qflx,   1, Ny, 1, Nt), ("Nx", "Ny", "Nz", "time"), Dict()),
        ("IET_wT",         reshape(IET_wT,     1, Ny, 1, Nt), ("Nx", "Ny", "Nz", "time"), Dict()),
        ("IET_nswflx",     reshape(IET_nswflx, 1, Ny, 1, Nt), ("Nx", "Ny", "Nz", "time"), Dict()),
        ("IET_swflx",      reshape(IET_swflx,  1, Ny, 1, Nt), ("Nx", "Ny", "Nz", "time"), Dict()),
        ("IET_dHdt",       reshape(IET_dHdt,   1, Ny, 1, Nt), ("Nx", "Ny", "Nz", "time"), Dict()),
        ("neb",            reshape(neb,        1, Ny, 1, Nt), ("Nx", "Ny", "Nz", "time"), Dict()),
        ("Q_clim",         reshape(Q_clim,     1, Ny, 1, Nt), ("Nx", "Ny", "Nz", "time"), Dict()),
        ("qflx",           reshape(qflx,       1, Ny, 1, Nt), ("Nx", "Ny", "Nz", "time"), Dict()),
        ("wT",             reshape(wT,         1, Ny, 1, Nt), ("Nx", "Ny", "Nz", "time"), Dict()),
        ("nswflx",         reshape(nswflx,     1, Ny, 1, Nt), ("Nx", "Ny", "Nz", "time"), Dict()),
        ("swflx",          reshape(swflx,      1, Ny, 1, Nt), ("Nx", "Ny", "Nz", "time"), Dict()),
        ("dHdt",           reshape(dHdt,       1, Ny, 1, Nt), ("Nx", "Ny", "Nz", "time"), Dict()),
    ]

        println("Doing var: ", varname)

        var = defVar(ds, varname, Float64, vardim)
        var.attrib["_FillValue"] = 1e20

        var = ds[varname]

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
using LinearAlgebra

struct BarzilaiBorweinIterable{T, F, R}
    w0::T
    f::F
    alpha::R
end

BarzilaiBorweinIterable(w0, f; alpha) = BarzilaiBorweinIterable(w0, f, alpha)

Base.IteratorSize(::Type{<:BarzilaiBorweinIterable}) = Base.IsInfinite()

mutable struct BarzilaiBorweinState{T, G, R}
    w::T
    f_w
    grad_f_w::G
    stepsize::R
end

function Base.iterate(iter::BarzilaiBorweinIterable)
    w = copy(iter.w0)
    grad_f_w, f_w = gradient(iter.f, w)
    state = BarzilaiBorweinState(w, f_w, grad_f_w, iter.alpha)
    return IterationOutput(state.w, state.f_w, state.grad_f_w), state
end

function Base.iterate(iter::BarzilaiBorweinIterable, state::BarzilaiBorweinState)
    state.w .-= state.stepsize .* state.grad_f_w
    y, state.f_w = gradient(iter.f, state.w)
    y .-= state.grad_f_w
    state.stepsize *= -dot(state.grad_f_w, state.grad_f_w) / dot(state.grad_f_w, y)
    state.grad_f_w .+= y
    return IterationOutput(state.w, state.f_w, state.grad_f_w), state
end

struct BarzilaiBorwein
    kwargs
    BarzilaiBorwein(; kwargs...) = new(kwargs)
end

(alg::BarzilaiBorwein)(args...) = BarzilaiBorweinIterable(args...; alg.kwargs...)


text:
# ---
# title: Scale
# cover: scale.gif
# ---

# Relatively resizing image

using Augmentor
using ImageShow, ImageCore
using Random

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

img_in = testpattern(RGB, ratio=0.5)

mosaicview(
    img_in,
    augment(img_in, Scale(0.8)),
    augment(img_in, Scale(0.8, 1));

    fillvalue=colorant"white", nrow=1, npad=10
)

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

Random.seed!(1337)
img_out = [augment(img_in, Scale(0.9:0.05:1.2)) for _ in 1:4]
mosaicview(img_out...; fillvalue=colorant"white", nrow=2)

# ## References

#md # ```@docs
#md # Scale
#md # ```

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