text:
### A Pluto.jl notebook ###
# v0.14.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 el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 7dea2baa-8f2c-11eb-184b-cb353ca6ca79
using PlutoUI

# ╔═╡ 9fd099fc-8f2c-11eb-286e-a548737f4572
using Plots

# ╔═╡ 9633d3f8-9000-11eb-3d5a-b3dd2a013890
using LinearAlgebra

# ╔═╡ 2c140db8-9013-11eb-21d9-f7718dc9cf57
using Optim

# ╔═╡ 56018eb0-90e1-11eb-18c0-e1f1aecf307c
using Vec

# ╔═╡ 5fead5d4-8f2d-11eb-2d3c-a34658637923
∇fline(x) = 0.75 * exp(-0.5 * (x - 7)^2) * (x - 7) + exp(-0.25 * (x - 3)^2) * 0.5 *(x - 3);

# ╔═╡ bf810ea8-8f2c-11eb-2d6e-3dcb7961b1fa
fline(x) = 1.5 - exp(-0.25 * (x - 3)^2)  - 0.75 * exp(-0.5 * (x - 7)^2);

# ╔═╡ 99e23b7a-8ffd-11eb-1b86-9940c74de5fc
rosenbrock(x) = (1 - x[1])^2 + 5 * (x[2] - x[1]^2)^2;

# ╔═╡ 9d46db42-8ffe-11eb-1f0f-1d6553bcd608
rosenbrock(x, y) = (1 - x)^2 + 5 * (y - x^2)^2;

# ╔═╡ 9a62773e-8ffd-11eb-315a-611790d3d945
∇rosenbrock(x) = [2 * (10 * x[1]^3 - 10 * x[1] * x[2] + x[1] - 1), 10 * (x[2] - x[1]^2)];

# ╔═╡ fb5587cc-9012-11eb-313a-512bd0f9651b
Hrosenbrock(x) = [-20 * (x[2] - x[1]^2) + 40 * x[1]^2 + 2 -20*x[1];
	              -20*x[1] 10];

# ╔═╡ fb7ea5ef-a2be-43a5-aae6-1bacc3a65f1e
rosenbrock2(x) = (1 - x[1])^2 + 100 * (4x[2] - x[1]^2)^2;

# ╔═╡ db569392-347f-45b1-b7c6-e2a0c235cff9
rosenbrock2(x, y) = (1-x)^2 + 100 * (4y - x^2)^2;

# ╔═╡ 1d09ca7a-8b25-4f3a-9f7d-ef586dd80be8
rosenbrock3(x, y) = (1-x)^2 + 100 * (y - x^2)^2;

# ╔═╡ cb524eeb-d78d-4125-8e5a-5da5e14644bb
∇rosenbrock2(x) = [2 * (200x[1]^3 - 800x[1] * x[2] + x[1] - 1), -800 * (x[1]^2 - 4x[2])];

# ╔═╡ 5e9fcab4-d2fe-445d-9b8e-2fe9fe50444d
md"""
# Descent Direction Iteration

1. Determine if out current x satisfies some termination condition. If so, terminate. If not, proceed.
2. Determine the **descent direction**.
3. Determine the **step size**.
4. Compute the next design point according to $\mathbf{x}^{(k+1)} \gets \mathbf{x}^{(k)} + \alpha^{(k)}\mathbf{d}^{(k)}$.
"""

# ╔═╡ 70855f33-c00d-49c3-82ec-f891ae4b40c5
md"""
nsteps = $(@bind nsteps1d NumberField(0:1:10, default = 0))

α = $(@bind α1d NumberField(0.1:0.1:3, default = 0.1))
"""

# ╔═╡ a3ae828c-3354-43e4-bf49-2284e465660d
begin
	p10 = plot(x -> x^2, legend = false, xlabel = "x", ylabel = "f(x)", color = :black, lw = 3, xlims = (-4,4), ylims = (0, 15))

	function run_gd(α, nsteps)
		gd_xs = [-3.5]
		gd_ys = [gd_xs[end]^2]

		for i = 1:nsteps
			next_x = gd_xs[end] - α * 2 * gd_xs[end]
			push!(gd_xs, next_x)
			push!(gd_ys, next_x^2)
		end

		return gd_xs, gd_ys
	end

	gd_xs, gd_ys = run_gd(α1d, nsteps1d)

	scatter!(gd_xs, gd_ys, color = :aqua, markersize = 5)
	scatter!(gd_xs, zeros(length(gd_xs)), color = :black, markersize = 5)
	for i = 1:length(gd_xs)
		plot!([gd_xs[i], gd_xs[i]], [0.0, gd_ys[i]], color = :black)
	end
	p10
end

# ╔═╡ 0259e325-9558-4dfc-9ca8-06daade248b9
md"""
# Multiple Dimensions
"""

# ╔═╡ e54c3874-667d-42a8-ba9f-5376f3ae3bd6
md"""
nsteps = $(@bind nsteps2d NumberField(0:1:10, default = 0))

α = $(@bind α2d NumberField(0.1:0.1:3, default = 0.1))
"""

# ╔═╡ 33f938e6-8229-4baa-9e1c-3c0215c53197
begin
	xvec = -4:0.1:4
	yvec = -4:0.1:4
	fbowl(x, y) = (x^2 + y^2) / 2

	function run_gd_2d(α, nsteps)
		gd_xs = [-3.5]
		gd_ys = [3.5]
		gd_zs = [fbowl(gd_xs[end], gd_ys[end])]

		for i = 1:nsteps
			next_x = gd_xs[end] - α * 2 * gd_xs[end]
			next_y = gd_ys[end] - α * 2 * gd_ys[end]
			push!(gd_xs, next_x)
			push!(gd_ys, next_y)
			push!(gd_zs, fbowl(gd_xs[end], gd_ys[end]))
		end

		return gd_xs, gd_ys, gd_zs
	end

	gd_xs2, gd_ys2, gd_zs2 = run_gd_2d(α2d, nsteps2d)

	sp1 = plot(xvec, yvec, fbowl, st = :surface, c = cgrad(:viridis, rev = true), colorbar = false, size = (600, 300), alpha = 0.5, legend = false, xlabel = "x₁", ylabel = "x₂", zlabel = "f(x)")
	#for i = 1:length(gd_xs)
	scatter!(sp1, gd_xs2, gd_ys2, gd_zs2, color = :red)
	plot!(sp1, gd_xs2, gd_ys2, gd_zs2, color = :black)
		#plot!(sp1, [-3.5, -3.5], [3.5, 3.5], [0.0, 9.0])

	sp2 = contour(xvec, yvec, fbowl, c = cgrad(:viridis, rev = true), colorbar = false, aspectratio = :equal, size = (600, 300), xlabel = "x₁", ylabel = "x₂", legend = false)
	scatter!(sp2, gd_xs2, gd_ys2, color = :red)
	plot!(sp2, gd_xs2, gd_ys2, color = :black)

	plot(sp1, sp2)
end

# ╔═╡ 9402e2b6-757e-4632-aac2-7251b28a768a
# begin
# # 	# plot(x -> x^2, legend = false, xlabel = "x", ylabel = "f(x)", color = :black, lw = 3, xlims = (-4,4), ylims = (0, 15))
# 	sp10 = plot(xvec, yvec, fbowl, st = :surface, c = cgrad(:viridis, rev = true), colorbar = false, size = (600, 400), alpha = 0.5, legend = false, xlabel = "x₁", ylabel = "x₂", zlabel = "f(x)")
# 	#scatter!(sp10, [-3.5], [3.5], [fbowl(-3.5, 3.5)], color = :red, markersize = 6)
# 	xsline = collect(-4:0.1:4)
# 	ysline = -xsline
# 	fsline = [fbowl(xsline[i], ysline[i]) for i = 1:length(xsline)]

# 	plot!(xsline, ysline, fsline, color = :white, lw = 6)
# 	plot!(xsline, ysline, zeros(length(xsline)), color = :black, lw = 6)
# end

# ╔═╡ a364ef52-8f2c-11eb-26a6-5d47e19a7230
md"""
# Line Search
"""

# ╔═╡ 24b2fe6c-9004-11eb-1e5f-0317dd2e85aa
md"""
### First Wolfe Condition
"""

# ╔═╡ 4794ee64-8f34-11eb-1ef2-5981f0e059df
β = 0.2

# ╔═╡ a00a4ec4-8f34-11eb-3db7-49be2cca0e5f
x₀ = 1.5

# ╔═╡ 80e1c23c-8f2e-11eb-3beb-519be0766468
md"""
α = $(@bind αlinedemo Slider(0.05:0.05:5.0, default = 5.0))
"""

# ╔═╡ ed8b99da-8f2c-11eb-14c5-0393c6d531c4
begin
	p1 = plot(fline, xlims = (0, 7.0), ylims = (0.45, 1.65), color = :black, xlabel = "x", ylabel = "f(x)", label = "")
	first_wolf_bound = fline(x₀) + β * αlinedemo * ∇fline(x₀)
	color = fline(x₀ + αlinedemo) < first_wolf_bound ? :green : :red
	plot!(Shape([0.0, 7.0, 7.0, 0.0, 0.0], [0.0, 0.0, first_wolf_bound, first_wolf_bound, 0.0]), color = color, alpha = 0.3, label = "")
	scatter!([x₀], [fline(x₀)], color = :black, label = "x₀")
	scatter!([x₀ + αlinedemo], [fline(x₀ + αlinedemo)], color = :gray, label = "f(x₀ + αd)")
	if fline(x₀ + αlinedemo) < first_wolf_bound
		xs = collect(x₀ + αlinedemo:0.05:4.05)
		ys = fline.(xs)
		plot!(xs, ys, color = :green, linewidth = 6, label = "")
	end
	p1
end

# ╔═╡ 4c0b5eac-8f35-11eb-06b4-a75f220c7f78
md"""
**f(x₀ + α):** $(fline(x₀ + αlinedemo))

**f(x₀) + βα∇f(x₀):** $first_wolf_bound
"""

# ╔═╡ 3954970e-9004-11eb-3be9-5f7adea26bc2
md"""
### Backtracking Line Search
"""

# ╔═╡ 5479f7d8-8fee-11eb-134b-a781d0affff4
function backtracking_line_search(f, ∇f, x, d, α; p = 0.5, β = 1e-4)
	y, g = f(x), ∇f(x)
	while f(x + α*d) > y + β*α*(g'd)
		α *= p
	end
	α
end

# ╔═╡ f87f64da-8fee-11eb-15fe-8f8090d4d226
md"""
p = $(@bind p NumberField(0.3:0.1:0.9, default = 0.8))
"""

# ╔═╡ ecaf9fe4-8fee-11eb-1c05-d76ef0f7ee1e
begin
	plot(fline, xlims = (0, 7.0), ylims = (0.45, 1.65), color = :black, xlabel = "x", ylabel = "f(x)", label = "", legend = false)

	xs1 = collect(x₀:0.05:4.05)
	ys1 = fline.(xs1)

	plot!(xs1, ys1, color = :green, linewidth = 6)

	function get_trace()
		αbt = 5.0
		attempted_xs = [x₀ + αbt]
		attempted_fs = [fline(x₀ + αbt)]

		while fline(x₀ + αbt) > fline(x₀) + β * αbt * ∇fline(x₀)
			αbt *= p
			push!(attempted_xs, x₀ + αbt)
			push!(attempted_fs, fline(x₀ + αbt))
		end
		return attempted_xs, attempted_fs
	end

	attempted_xs, attempted_fs = get_trace()

	scatter!(attempted_xs[1:end-1], attempted_fs[1:end-1], color = :red, markersize = 6)
	scatter!([attempted_xs[end]], [attempted_fs[end]], color = :green, markersize = 6)
	scatter!([x₀], [fline(x₀)], color = :black, label = "x₀", markersize = 6)
end

# ╔═╡ 827f2e30-8ffd-11eb-2cd6-9b1edfa701f9
md"""
### Backtracking line search on the Rosenbrock function
"""

# ╔═╡ 3f536d2c-9003-11eb-23c7-7b0b0a92d51a
md"""
nsteps = $(@bind nstepsbls NumberField(0:1:7, default = 0))
"""

# ╔═╡ e0fea120-8ffd-11eb-2cac-4f6e4547cc26
begin
	xr = -2:0.1:2
	yr = -2:0.1:2

	# Contour of rosenbrock
	c = contour(xr, yr, rosenbrock, levels=[0.5,1,2,3,5,10,20,50,100], colorbar = false, c = cgrad(:viridis, rev = true), legend = false, xlims = (-2, 2), ylims = (-2, 2), xlabel = "x₁", ylabel = "x₂", aspectratio = :equal)
	xbls = [-1.75, -1.75]
	scatter!([xbls[1]], [xbls[2]], color = :black)

	pbls = 0.5
	βbls = 1e-4
	αmaxbls = 6.0
	#nstepsbl = 7
	for k in 1:nstepsbls
		d = -∇rosenbrock(xbls)
		α = αmaxbls
		y, g = rosenbrock(xbls), ∇rosenbrock(xbls)

		dir = normalize(d, 2)
		plot!([xbls[1], xbls[1] + 6dir[1]], [xbls[2], xbls[2] + 6dir[2]], color = :red)

		while rosenbrock(xbls + α*d) ≥ y + βbls*α*(g⋅d)
			x_push = xbls + d*α
			if -2 ≤ x_push[1] ≤ 2 && -2 ≤ x_push[2] ≤ 2
				scatter!([x_push[1]], [x_push[2]], color = :red)
			end
			α *= pbls
		end
		scatter!([xbls[1] + d[1]*α], [xbls[2] + d[2]*α], color = :black)
		plot!([xbls[1], xbls[1] + d[1]*α], [xbls[2], xbls[2] + d[2]*α], color = :black)
		xbls += d*α
	end
	c
end

# ╔═╡ 4848aa34-9004-11eb-1408-5598d9e145f3
md"""
### Curvature Condition
"""

# ╔═╡ 3ae1a64c-9005-11eb-32cd-adbd7cb6651a
md"""
σ = $(@bind σcurve NumberField(0.1:0.1:1, default = 0.4))
"""

# ╔═╡ 55f1ca30-9004-11eb-2658-972069811aa5
begin
	p2 = plot(fline, xlims = (0, 7.0), ylims = (0.45, 1.65), color = :black, xlabel = "x", ylabel = "f(x)", label = "")
	scatter!([x₀], [fline(x₀)], color = :black, label = "x₀")

	gradx₀ = ∇fline(x₀)
	points_curv = collect(x₀:0.05:7)
	gradpoints_curv = ∇fline.(points_curv)
	inds_curv = findall(gradpoints_curv .≥ σcurve * gradx₀)

	diffs_inds_curv = [inds_curv[i+1] - inds_curv[i] for i = 1:length(inds_curv)-1]
	breaks_curv = findall(diffs_inds_curv .> 1)

	if breaks_curv != nothing
		for i = 1:length(breaks_curv)
			if i == 1
				start_ind = 1
			else
				start_ind = breaks_curv[i-1] + 1
			end
			plot_inds = inds_curv[start_ind:breaks_curv[i]]
			plot!(points_curv[plot_inds], fline.(points_curv[plot_inds]), color = :green, linewidth = 6, label = "")
		end
		plot!(points_curv[inds_curv[breaks_curv[end]+1:end]], fline.(points_curv[inds_curv[breaks_curv[end]+1:end]]), color = :green, linewidth = 6, label = "")
	else
		plot!(points_curv[inds_curv], fline.(points_curv[inds_curv]), color = :green, linewidth = 6, label = "")
	end
	p2
end

# ╔═╡ 529ef5e6-9007-11eb-05b2-a70d8ad51908
md"""
### Strong Curvature Condition (Second Wolfe Condition)
"""

# ╔═╡ 9afa6fa8-9007-11eb-0d4e-01881e229cc7
md"""
σ = $(@bind σscurve NumberField(0.1:0.1:1, default = 0.4))
"""

# ╔═╡ 63249a1a-9007-11eb-2630-0523c391ee59
begin
	p3 = plot(fline, xlims = (0, 7.0), ylims = (0.45, 1.65), color = :black, xlabel = "x", ylabel = "f(x)", label = "")
	scatter!([x₀], [fline(x₀)], color = :black, label = "x₀")

	points_scurv = collect(x₀:0.05:7)
	gradpoints_scurv = ∇fline.(points_scurv)
	inds_scurv = findall(abs.(gradpoints_scurv) .≤ -σscurve * gradx₀)

	diffs_inds_scurv = [inds_scurv[i+1] - inds_scurv[i] for i = 1:length(inds_scurv)-1]
	breaks_scurv = findall(diffs_inds_scurv .> 1)

	if breaks_scurv != nothing
		for i = 1:length(breaks_scurv)
			if i == 1
				start_ind = 1
			else
				start_ind = breaks_scurv[i-1] + 1
			end
			plot_inds = inds_scurv[start_ind:breaks_scurv[i]]
			plot!(points_scurv[plot_inds], fline.(points_scurv[plot_inds]), color = :green, linewidth = 6, label = "")
		end
		plot!(points_scurv[inds_scurv[breaks_scurv[end]+1:end]], fline.(points_scurv[inds_scurv[breaks_scurv[end]+1:end]]), color = :green, linewidth = 6, label = "")
	else
		plot!(points_scurv[inds_scurv], fline.(points_scurv[inds_scurv]), color = :green, linewidth = 6, label = "")
	end
	p3
end

# ╔═╡ e7ca1b5c-9008-11eb-2e5b-972ffbcdad4b
md"""
### Wolfe Conditions
"""

# ╔═╡ 0201e428-9009-11eb-3b2f-d7f3bb5b1831
begin
	p4 = plot(fline, xlims = (0, 7.0), ylims = (0.45, 1.65), color = :black, xlabel = "x", ylabel = "f(x)", label = "")
	scatter!([x₀], [fline(x₀)], color = :black, label = "x₀")

	points_wolfe = collect(x₀:0.05:7)
	gradpoints_wolfe = ∇fline.(points_wolfe)
	inds_wolfe1 = findall(fline.(points_wolfe) .≤ fline(x₀) .+ (β .* (points_wolfe .- x₀) .* ∇fline(x₀)))
	inds_wolfe2 = findall(abs.(gradpoints_wolfe) .≤ -σscurve * gradx₀)

	diffs_inds_wolfe2 = [inds_wolfe2[i+1] - inds_wolfe2[i] for i = 1:length(inds_wolfe2)-1]
	breaks_wolfe2 = findall(diffs_inds_wolfe2 .> 1)

	if breaks_wolfe2 != nothing
		for i = 1:length(breaks_wolfe2)
			if i == 1
				start_ind = 1
			else
				start_ind = breaks_wolfe2[i-1] + 1
			end
			plot_inds = inds_wolfe2[start_ind:breaks_wolfe2[i]]
			plot!(points_wolfe[plot_inds], fline.(points_wolfe[plot_inds]), color = :red, linewidth = 6, label = "")
		end
		plot!(points_wolfe[inds_wolfe2[breaks_wolfe2[end]+1:end]], fline.(points_wolfe[inds_wolfe2[breaks_wolfe2[end]+1:end]]), color = :red, linewidth = 6, label = "2nd Wolfe")
	else
		plot!(points_wolfe[inds_wolfe2], fline.(points_wolfe[inds_wolfe2]), color = :red, linewidth = 6, label = "2nd Wolfe")
	end

	plot!(points_wolfe[inds_wolfe1], fline.(points_wolfe[inds_wolfe1]), color = :blue, linewidth = 6, label = "1st Wolfe")

	inds_both_wolfe = intersect(inds_wolfe1, inds_wolfe2)

	plot!(points_wolfe[inds_both_wolfe], fline.(points_wolfe[inds_both_wolfe]), color = :green, linewidth = 6, label = "Both")

	p4
end

# ╔═╡ 631fc672-900b-11eb-1fe2-a531dac443b1
md"""
### Strong Backtracking Line Search
"""

# ╔═╡ a738a7fb-0ebb-42c5-bc2d-af1a4ea72bfc
md"""
If one of the following conditions hold, an interval is guaranteed to contain step lengths satisfying the Wolfe conditions:

$f(\mathbf{x} + \alpha^{(k)}\mathbf{d}) \geq f(\mathbf{x})$

$f(\mathbf{x^{(k)}} + \alpha^{(k)}\mathbf{d^{(k)}}) \geq f(\mathbf{x^{(k)}}) + \beta\alpha\nabla_{\mathbf{d}^{(k)}}f(\mathbf{x}^{(k)})$

$\nabla f(\mathbf{x} + \alpha^{(k)}\mathbf{d}) \geq \mathbf{0}$
"""

# ╔═╡ 272c43fc-9012-11eb-07f6-3d0625cafbab
md"""
nsteps = $(@bind nstepssbls NumberField(0:1:7, default = 0))
"""

# ╔═╡ 75b076d0-900b-11eb-2875-b9b74ab9a943
begin
	x_init = 0.0
	cond1(α_hi) = ∇fline(x_init + α_hi) >= 0
	cond2(α_hi) = fline(x_init + α_hi) > fline(x_init) + β*α_hi*∇fline(x_init)
	cond3(α_hi) = fline(x_init + α_hi) >= fline(x_init)

	p5 = plot(fline, xlims = (0, 7.0), ylims = (0.45, 1.65), color = :black, xlabel = "x", ylabel = "f(x)", label = "")
	scatter!([x_init], [fline(x_init)], color = :black, label = "x₀", markersize = 6)

	points_cd = collect(0.5:0.05:7)
	cd1 = cond1.(points_cd)
	cd2 = cond2.(points_cd)
	cd3 = cond3.(points_cd)
	inds_cd = findall(cd1 .| cd2 .| cd3)

	diffs_inds_cd = [inds_cd[i+1] - inds_cd[i] for i = 1:length(inds_cd)-1]
	breaks_cd = findall(diffs_inds_cd .> 1)

	if breaks_cd != nothing
		for i = 1:length(breaks_cd)
			if i == 1
				start_ind = 1
			else
				start_ind = breaks_cd[i-1] + 1
			end
			plot_inds = inds_cd[start_ind:breaks_cd[i]]
			plot!(points_cd[plot_inds], fline.(points_cd[plot_inds]), color = :blue, linewidth = 6, label = "")
		end
		plot!(points_cd[inds_cd[breaks_cd[end]+1:end]], fline.(points_cd[inds_cd[breaks_cd[end]+1:end]]), color = :blue, linewidth = 6, label = "Interval Conditions")
	else
		plot!(points_cd[inds_cd], fline.(points_cd[inds_cd]), color = :blue, linewidth = 6, label = "Interval Conditions")
	end

	function sbls()
		α = 0.6
		y0, g0, y_prev, α_prev = fline(x_init), ∇fline(x_init), NaN, 0

		for k = 1:nstepssbls
			y = fline(x_init + α)

			scatter!([α], [y], color = :black, label = "", markersize = 6)

			if y > y0 + β*α*g0 || (!isnan(y_prev) && y >= y_prev)

				α = 0.75α # bisection
				y = fline(α)
				scatter!([α], [y], color = :red, label = "", markersize = 6)
			end
			g = ∇fline(x_init + α)
			if abs(g) <= -σscurve * g0
				break
			elseif g >= 0
				break
			end
			α_prev, α = α, 2α
			y_prev = y
		end
	end

	sbls()
	p5
end

# ╔═╡ 81f2c770-9012-11eb-099f-dda2f9c4246f
md"""
# Trust Region Methods
"""

# ╔═╡ e9e89ffc-9018-11eb-1028-856604a9f26b
md"""
nsteps = $(@bind nstepstrust NumberField(0:1:10, default = 0))
"""

# ╔═╡ 447b2884-9013-11eb-1457-0b057903165f
begin
	function solve_trust_region_subproblem(∇f, H, x0, δ)

		Hx0 = H(x0)
		∇fx0 = ∇f(x0)

		# Unconstrained
		f = x -> ∇fx0⋅(x-x0) + ((x-x0)'*Hx0)⋅(x-x0)/2

		# Constrained with huge penalty
		g = x -> norm(x-x0) <= δ ? f(x) : 9999999.0

		result = optimize(g, x0, NelderMead())

		return (result.minimizer, result.minimum)
	end

	# Contour of rosenbrock
	cr = contour(xr, yr, rosenbrock, levels=[0.5,1,2,3,5,10,20,50,100], colorbar = false, c = cgrad(:viridis, rev = true), legend = false, xlims = (-2, 2), ylims = (-2, 2), xlabel = "x₁", ylabel = "x₂", aspectratio = :equal)

	#scatter!([xbls[1]], [xbls[2]], color = :blue)

	function trust_region_rosenbrock(k_max)
		x₀ = [-1.75,-1.75]
		x = deepcopy(x₀)
		y = rosenbrock(x)
		δ = 1.0
		θs = range(0,stop=2π,length=101)
		η1=0.25; η2=0.5; γ1=0.5; γ2=2.0; δ=1.0

		for k = 1:k_max
			scatter!([x[1]], [x[2]], color = :black, alpha = k/k_max)
			plot!([x[1] + δ*cos(θ) for θ in θs], [x[2] + δ*sin(θ) for θ in θs], color = :black, alpha = k/k_max)

			x′, y′ = solve_trust_region_subproblem(∇rosenbrock, Hrosenbrock, x, δ)
			r = (y - rosenbrock(x′)) / (y - y′)
			if r < η1
				δ *= γ1
			else
				x, y = x′, rosenbrock(x′)
				if r > η2
					δ *= γ2
				end
			end
		end
	end

	trust_region_rosenbrock(nstepstrust)

	cr
end

# ╔═╡ 82a6b55a-90df-11eb-119b-9d94c9420aae
md"""
# First-Order Methods
"""

# ╔═╡ 9c3268ac-90df-11eb-3a09-cd5b1e01213d
md"""
### Gradient Descent
"""

# ╔═╡ 9397f2be-90e1-11eb-2a0c-559077244bd9
md"""
nsteps = $(@bind nstepsgd NumberField(0:1:10, default = 0))
"""

# ╔═╡ a9d3e5b2-90df-11eb-388f-eb7a9433a6dd
begin
	# Contour of rosenbrock
	cr1 = contour(xr, yr, rosenbrock, levels=[0.5,1,2,3,5,10,20,50,100], colorbar = false, c = cgrad(:viridis, rev = true), legend = false, xlims = (-2, 2), ylims = (-2, 2), xlabel = "x₁", ylabel = "x₂", aspectratio = :equal)
	scatter!([-1.0], [-1.0], color = :black)

	# For solving for α
	function secant_method(df, x1, x2, ϵ)
	    df1 = df(x1)

	    delta = Inf
	    while abs(delta) > ϵ
	    	df2 = df(x2)
	        delta = (x2 - x1)/(df2 - df1)*df2
	        x1, x2, df1 = x2, x2 - delta, df2
	    end
	    x2
	end

	# Run
	df(x,y) = [2*(10*x^3-10*x*y+x-1), 10*(y-x^2)]

	function grad_descent(nsteps)
		p0 = (-1, -1)
		pts_grad = Tuple{Float64,Float64}[p0]
		for i in 1 : nsteps
			x, y = pts_grad[end]
			dp = normalize(-VecE2{Float64}(df(x, y)...))

			f1d = a -> begin
				x2 = x + a*dp.x
				y2 = y + a*dp.y

				da = df(x2, y2)
				pa = VecE2{Float64}(da[1], da[2])
				proj(pa, dp, Float64)
			end
			alpha = secant_method(f1d, 0.0, 1.0, 0.0001)
			push!(pts_grad, (x + alpha*dp.x, y + alpha*dp.y))
		end
		return pts_grad
	end

	pts_grad = grad_descent(10)

	plot!([pts_grad[i][1] for i = 1:nstepsgd+1], [pts_grad[i][2] for i = 1:nstepsgd+1], color = :black, linewidth = 2)
end

# ╔═╡ d59fd8f2-90e1-11eb-102a-f36334c851b4
md"""
### Conjugate Gradient Descent
"""

# ╔═╡ 0263d884-90e7-11eb-2d19-37a0ebcf527b
md"""
nsteps = $(@bind nstepsgd_conj NumberField(0:1:10, default = 0))
"""

# ╔═╡ c589570a-90e5-11eb-2100-45f865a96a02
begin
	# Contour of rosenbrock
	cr2 = contour(xr, yr, rosenbrock, levels=[0.5,1,2,3,5,10,20,50,100], colorbar = false, c = cgrad(:viridis, rev = true), legend = false, xlims = (-2, 2), ylims = (-2, 2), xlabel = "x₁", ylabel = "x₂", aspectratio = :equal)
	scatter!([-1.0], [-1.0], color = :black)

	function conj_grad(nsteps)
		x0 = Float64[-1, -1]
		pts_conj = Vector{Float64}[x0]
		g0 = ∇rosenbrock(x0)
		d0 = -g0
		f1d = a -> begin
			da = ∇rosenbrock(x0 + a*d0)
			proj(VecE2{Float64}(da[1], da[2]), VecE2{Float64}(d0[1], d0[2]), Float64)
		end
		alpha = secant_method(f1d, 0.0, 0.1, 0.0001) # NOTE: is sensitive to the value of 0.1, does not work with 1.0!
		push!(pts_conj, x0 + alpha*d0)

		for i in 1 : nsteps
			x = pts_conj[end]
			g1 = ∇rosenbrock(x)
			β = max(0.0, dot(g1, g1-g0)/dot(g0,g0))
			d1 = -g1 + β*d0

			f1d = a -> begin
				da = ∇rosenbrock(x + a*d1)
				proj(VecE2{Float64}(da[1], da[2]), VecE2{Float64}(d1[1], d1[2]), Float64)
			end
			alpha = secant_method(f1d, 0.0, 0.1, 0.0001)
			push!(pts_conj, x + alpha*d1)

			d0 = d1
			g0 = g1
		end
		return pts_conj
	end

	pts_conj = conj_grad(10)

	plot!([pts_grad[i][1] for i = 1:nstepsgd_conj+1], [pts_grad[i][2] for i = 1:nstepsgd_conj+1], color = :black, linewidth = 2, alpha = 0.5)
	plot!([pts_conj[i][1] for i = 1:nstepsgd_conj+1], [pts_conj[i][2] for i = 1:nstepsgd_conj+1], color = :black, linewidth = 2)
end

# ╔═╡ 04834139-2837-4891-9ba1-4f16c32df04f
md"""
### Momentum
"""

# ╔═╡ 50688a3b-c1bf-4cde-a820-874158384fae
md"""
nsteps = $(@bind nstepsgd_mom NumberField(0:1:40, default = 0))
"""

# ╔═╡ 7967346e-bf2f-4791-8e5a-11fee51ba9de
begin
	xr2 = -3:0.1:2
	yr2 = -0.5:0.1:2
	# Contour of rosenbrock 2
	cr3 = contour(xr2, yr2, rosenbrock2, levels=[25,50,100,200,250,300], colorbar = false, c = cgrad(:viridis, rev = true), legend = false, xlims = (-3, 2), ylims = (-0.5, 2), xlabel = "x₁", ylabel = "x₂", aspectratio = :equal, clim = (2, 500))
	scatter!([-2.0], [1.5], color = :black)

	abstract type DescentMethod end

	function this_step!(M::DescentMethod, v::VecE2{Float64})
	    x = Float64[v.x, v.y]
	    return VecE2{Float64}(step!(M, rosenbrock2, ∇rosenbrock2, x)...)
	end

	function run_descent_method(M::DescentMethod, x₀::VecE2{Float64}, N::Int)
	    pts = [x₀]
	    init!(M, rosenbrock2, ∇rosenbrock2, Float64[x₀.x, x₀.y])
	    for i in 1 : N
	        push!(pts, this_step!(M, pts[end]))
	    end
	    return pts
	end

	struct GradientDescent <: DescentMethod
		α
	end

	function init!(M::GradientDescent, f, ∇f, x)
			return M
	end

	function step!(M::GradientDescent, f, ∇f, x)
		α, g = M.α, ∇f(x)
		return x - α*g
	end

	mutable struct Momentum <: DescentMethod
		α # learning rate
		β # momentum decay
		v # momentum
	end

	function init!(M::Momentum, f, ∇f, x)
		M.v = zeros(length(x))
		return M
	end

	function step!(M::Momentum, f, ∇f, x)
		α, β, v, g = M.α, M.β, M.v, ∇f(x)
		v[:] = β*v - α*g
		return x + v
	end

	mutable struct NesterovMomentum <: DescentMethod
		α # learning rate
		β # momentum decay
		v # momentum
	end

	function init!(M::NesterovMomentum, f, ∇f, x)
		M.v = zeros(length(x))
		return M
	end

	function step!(M::NesterovMomentum, f, ∇f, x)
		α, β, v = M.α, M.β, M.v
		v[:] = β*v - α*∇f(x + β*v)
		return x + v
	end

	mutable struct HyperGradientDescent <: DescentMethod
		α0 # initial learning rate
		μ # learning rate of the learning rate
		α # current learning rate
		g_prev # previous gradient
	end

	function init!(M::HyperGradientDescent, f, ∇f, x)
		M.α = M.α0
		M.g_prev = zeros(length(x))
		return M
	end

	function step!(M::HyperGradientDescent, f, ∇f, x)
		α, μ, g, g_prev = M.α, M.μ, ∇f(x), M.g_prev
		α = α + μ*(g⋅g_prev)
		M.g_prev, M.α = g, α
		return x - α*g
	end

	mutable struct HyperNesterovMomentum <: DescentMethod
		α0 # initial learning rate
		μ # learning rate of the learning rate
		β # momentum decay
		v # momentum
		α # current learning rate
		g_prev # previous gradient
	end

	function init!(M::HyperNesterovMomentum, f, ∇f, x)
		M.α = M.α0
		M.v = zeros(length(x))
		M.g_prev = zeros(length(x))
		return M
	end

	function step!(M::HyperNesterovMomentum, f, ∇f, x)
		α, β, μ = M.α, M.β, M.μ
		v, g, g_prev = M.v, ∇f(x), M.g_prev
		α = α - μ*(g⋅(-g_prev - β*v))
		v[:] = β*v + g
		M.g_prev, M.α = g, α
		return x - α*(g + β*v)
	end

	N = 40

	pts_gd = run_descent_method(GradientDescent(0.0003), VecE2{Float64}(-2,1.5), N)
	pts_mom = run_descent_method(Momentum(0.0003, 0.9, zeros(2)), VecE2{Float64}(-2,1.5), N)

	plot!([pts_gd[i][1] for i = 1:nstepsgd_mom + 1], [pts_gd[i][2] for i = 1:nstepsgd_mom + 1], color = :black, linewidth = 2, alpha = 0.5)
	plot!([pts_mom[i][1] for i = 1:nstepsgd_mom + 1], [pts_mom[i][2] for i = 1:nstepsgd_mom + 1], color = :black, linewidth = 2)
end

# ╔═╡ bf301087-6624-4c56-a8d2-dc7a9856d33e
md"""
### Nesterov Momentum
"""

# ╔═╡ 5c4506ab-c281-4f94-a60a-ef5d727fe125
md"""
nsteps = $(@bind nstepsnest_mom NumberField(0:1:40, default = 0))
"""

# ╔═╡ b638d908-35a9-4203-945a-1086f4dba559
begin
	# Contour of rosenbrock 2
	cr4 = contour(xr2, yr2, rosenbrock2, levels=[25,50,100,200,250,300], colorbar = false, c = cgrad(:viridis, rev = true), legend = false, xlims = (-3, 2), ylims = (-0.5, 2), xlabel = "x₁", ylabel = "x₂", aspectratio = :equal, clim = (2, 500))
	scatter!([-2.0], [1.5], color = :black)

	pts_nest_mom = run_descent_method(NesterovMomentum(0.0002, 0.92, zeros(2)), VecE2{Float64}(-2,1.5), N)

	plot!([pts_mom[i][1] for i = 1:nstepsnest_mom + 1], [pts_mom[i][2] for i = 1:nstepsnest_mom + 1], color = :black, linewidth = 2, alpha = 0.5)
	plot!([pts_nest_mom[i][1] for i = 1:nstepsnest_mom + 1], [pts_nest_mom[i][2] for i = 1:nstepsnest_mom + 1], color = :black, linewidth = 2)
end

# ╔═╡ d0a4ca13-8ea4-4451-ba88-b46d28a1ba89
md"""
### Hypergradient Descent and Hyper-Nesterov
"""

# ╔═╡ b3759a50-11d4-4d96-a80b-bcb6a76ba944
md"""
nsteps = $(@bind nstepshyper NumberField(0:1:40, default = 0))
"""

# ╔═╡ 36099212-4d39-4e22-a142-0b9b88263de2
begin
	# Contour of rosenbrock 2
	cr5 = contour(xr2, yr2, rosenbrock2, levels=[25,50,100,200,250,300], colorbar = false, c = cgrad(:viridis, rev = true), legend = false, xlims = (-3, 2), ylims = (-0.5, 2), xlabel = "x₁", ylabel = "x₂", aspectratio = :equal, clim = (2, 500))
	scatter!([-2.0], [1.5], color = :black)

	pts_hyp_gd = run_descent_method(HyperGradientDescent(0.0004, 8e-13, NaN, zeros(2)), VecE2{Float64}(-2,1.5), N)
	pts_hyp_nest = run_descent_method(HyperNesterovMomentum(0.00023, 1e-12, 0.93, zeros(2), NaN, zeros(2)), VecE2{Float64}(-2,1.5), N)

	plot!([pts_hyp_gd[i][1] for i = 1:nstepshyper + 1], [pts_hyp_gd[i][2] for i = 1:nstepshyper + 1], color = :black, linewidth = 2, alpha = 0.5)
	plot!([pts_hyp_nest[i][1] for i = 1:nstepshyper + 1], [pts_hyp_nest[i][2] for i = 1:nstepshyper+ 1], color = :black, linewidth = 2)
end

# ╔═╡ Cell order:
# ╠═7dea2baa-8f2c-11eb-184b-cb353ca6ca79
# ╠═9fd099fc-8f2c-11eb-286e-a548737f4572
# ╠═9633d3f8-9000-11eb-3d5a-b3dd2a013890
# ╠═2c140db8-9013-11eb-21d9-f7718dc9cf57
# ╠═56018eb0-90e1-11eb-18c0-e1f1aecf307c
# ╠═5fead5d4-8f2d-11eb-2d3c-a34658637923
# ╠═bf810ea8-8f2c-11eb-2d6e-3dcb7961b1fa
# ╠═99e23b7a-8ffd-11eb-1b86-9940c74de5fc
# ╠═9d46db42-8ffe-11eb-1f0f-1d6553bcd608
# ╠═9a62773e-8ffd-11eb-315a-611790d3d945
# ╠═fb5587cc-9012-11eb-313a-512bd0f9651b
# ╠═fb7ea5ef-a2be-43a5-aae6-1bacc3a65f1e
# ╠═db569392-347f-45b1-b7c6-e2a0c235cff9
# ╠═1d09ca7a-8b25-4f3a-9f7d-ef586dd80be8
# ╠═cb524eeb-d78d-4125-8e5a-5da5e14644bb
# ╟─5e9fcab4-d2fe-445d-9b8e-2fe9fe50444d
# ╟─70855f33-c00d-49c3-82ec-f891ae4b40c5
# ╟─a3ae828c-3354-43e4-bf49-2284e465660d
# ╟─0259e325-9558-4dfc-9ca8-06daade248b9
# ╟─e54c3874-667d-42a8-ba9f-5376f3ae3bd6
# ╟─33f938e6-8229-4baa-9e1c-3c0215c53197
# ╟─9402e2b6-757e-4632-aac2-7251b28a768a
# ╟─a364ef52-8f2c-11eb-26a6-5d47e19a7230
# ╟─24b2fe6c-9004-11eb-1e5f-0317dd2e85aa
# ╟─4794ee64-8f34-11eb-1ef2-5981f0e059df
# ╟─a00a4ec4-8f34-11eb-3db7-49be2cca0e5f
# ╟─80e1c23c-8f2e-11eb-3beb-519be0766468
# ╟─4c0b5eac-8f35-11eb-06b4-a75f220c7f78
# ╟─ed8b99da-8f2c-11eb-14c5-0393c6d531c4
# ╟─3954970e-9004-11eb-3be9-5f7adea26bc2
# ╠═5479f7d8-8fee-11eb-134b-a781d0affff4
# ╟─f87f64da-8fee-11eb-15fe-8f8090d4d226
# ╟─ecaf9fe4-8fee-11eb-1c05-d76ef0f7ee1e
# ╟─827f2e30-8ffd-11eb-2cd6-9b1edfa701f9
# ╟─3f536d2c-9003-11eb-23c7-7b0b0a92d51a
# ╟─e0fea120-8ffd-11eb-2cac-4f6e4547cc26
# ╟─4848aa34-9004-11eb-1408-5598d9e145f3
# ╟─3ae1a64c-9005-11eb-32cd-adbd7cb6651a
# ╟─55f1ca30-9004-11eb-2658-972069811aa5
# ╟─529ef5e6-9007-11eb-05b2-a70d8ad51908
# ╟─9afa6fa8-9007-11eb-0d4e-01881e229cc7
# ╟─63249a1a-9007-11eb-2630-0523c391ee59
# ╟─e7ca1b5c-9008-11eb-2e5b-972ffbcdad4b
# ╟─0201e428-9009-11eb-3b2f-d7f3bb5b1831
# ╟─631fc672-900b-11eb-1fe2-a531dac443b1
# ╟─a738a7fb-0ebb-42c5-bc2d-af1a4ea72bfc
# ╟─272c43fc-9012-11eb-07f6-3d0625cafbab
# ╟─75b076d0-900b-11eb-2875-b9b74ab9a943
# ╟─81f2c770-9012-11eb-099f-dda2f9c4246f
# ╟─e9e89ffc-9018-11eb-1028-856604a9f26b
# ╟─447b2884-9013-11eb-1457-0b057903165f
# ╟─82a6b55a-90df-11eb-119b-9d94c9420aae
# ╟─9c3268ac-90df-11eb-3a09-cd5b1e01213d
# ╟─9397f2be-90e1-11eb-2a0c-559077244bd9
# ╟─a9d3e5b2-90df-11eb-388f-eb7a9433a6dd
# ╟─d59fd8f2-90e1-11eb-102a-f36334c851b4
# ╟─0263d884-90e7-11eb-2d19-37a0ebcf527b
# ╟─c589570a-90e5-11eb-2100-45f865a96a02
# ╟─04834139-2837-4891-9ba1-4f16c32df04f
# ╟─50688a3b-c1bf-4cde-a820-874158384fae
# ╟─7967346e-bf2f-4791-8e5a-11fee51ba9de
# ╟─bf301087-6624-4c56-a8d2-dc7a9856d33e
# ╟─5c4506ab-c281-4f94-a60a-ef5d727fe125
# ╟─b638d908-35a9-4203-945a-1086f4dba559
# ╟─d0a4ca13-8ea4-4451-ba88-b46d28a1ba89
# ╟─b3759a50-11d4-4d96-a80b-bcb6a76ba944
# ╟─36099212-4d39-4e22-a142-0b9b88263de2


text:
using Plots;gr()

f(x) = x*x - 2
df(x) = 2x

function newton(x_0, max_iter)
    x = x_0
    error = zeros(max_iter)
    for i in 1:max_iter
        x = x - f(x) / df(x)
        error[i] = abs(x - sqrt(2))

    end

    (x, error)
end

x_0 = 1.0
max_iter = 20

x, error = newton(x_0, max_iter)
print(error)
t = 1:max_iter

plot(t, error, color=:cornflowerblue, label="newton",xlabel="t",ylabel="Error")
savefig("newton.png")

ext: jl


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ATKmm_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ATKmm")
JLLWrappers.@generate_main_file("ATKmm", UUID("ae83496e-abbd-5b27-8e4c-ced0103e1cfe"))
end  # module ATKmm_jll


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Test
using JuliaFEM: get_polygon_clip, calculate_polygon_area

## polygon clipping

S = [[0.375, 0.000, 0.500], [0.600, 0.000, 0.500], [0.500, 0.250, 0.500]]
M = [[0.500, 0.000, 0.500], [0.250, 0.000, 0.500], [0.375, 0.250, 0.500]]
n0 = [0.0, 0.0, 1.0]
P = get_polygon_clip(S, M, n0)
@test length(P) == 3
@test isapprox(calculate_polygon_area(P), 1/128)

S = [[0.250, 0.000, 0.500], [0.750, 0.000, 0.500], [0.500, 0.250, 0.500]]
M = [[0.500, 0.000, 0.500], [0.250, 0.000, 0.500], [0.375, 0.250, 0.500]]
n0 = [0.0, 0.0, 1.0]
P = get_polygon_clip(S, M, n0)
@test length(P) == 3
@test isapprox(calculate_polygon_area(P), 1/48)

# visually inspected
Xs = [[ 0.00, 0.00, 0.50], [1.00,  0.00, 0.50], [0.00, 1.00, 0.50]]
Xm = [[-0.25, 0.50, 0.50], [0.50, -0.25, 0.50], [0.75, 0.75, 0.50]]
P_ = [[0.65, 0.35, 0.00], [0.5625, 0.0,    0.0], [0.25, 0.00, 0.0],
      [0.00, 0.25, 0.00], [0.0,    0.5625, 0.0], [0.35, 0.65, 0.0]]
P = get_polygon_clip(Xs, Xm, [0.0, 0.0, 1.0])
@test length(P) == length(P_)


text:
type RHS_IE_Scalar{F,uType,tType} <: Function
  f::F
  u_old::uType
  t::tType
  dt::tType
end

function (p::RHS_IE_Scalar)(u,resid)
  resid[1] = u[1] - p.u_old[1] - p.dt*p.f(p.t+p.dt,u[1])[1]
end

@inline function initialize!(integrator,cache::ImplicitEulerConstantCache,f=integrator.f)
  cache.uhold[1] = integrator.uprev; cache.u_old[1] = integrator.uprev
  integrator.kshortsize = 2
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  integrator.fsalfirst = f(integrator.t,integrator.uprev)
end

@inline function perform_step!(integrator,cache::ImplicitEulerConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  @unpack uhold,u_old,rhs,adf = cache
  u_old[1] = uhold[1]
  if integrator.iter > 1 && !integrator.u_modified
    uhold[1] = current_extrapolant(t+dt,integrator)
  end # else uhold is previous value.
  rhs.t = t
  rhs.dt = dt
  if alg_autodiff(integrator.alg)
    nlres = NLsolve.nlsolve(adf,uhold)
  else
    nlres = NLsolve.nlsolve(rhs,uhold,autodiff=alg_autodiff(integrator.alg))
  end
  uhold[1] = nlres.zero[1]
  k = f(t+dt,uhold[1])
  integrator.fsallast = k
  u = uhold[1]
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = k
  @pack integrator = t,dt,u
end

type RHS_IE{F,uType,tType,DiffCacheType,SizeType,uidxType} <: Function
  f::F
  u_old::uType
  t::tType
  dt::tType
  dual_cache::DiffCacheType
  sizeu::SizeType
  uidx::uidxType
end
function (p::RHS_IE)(uprev,resid)
  du = get_du(p.dual_cache, eltype(uprev))
  p.f(p.t+p.dt,reshape(uprev,p.sizeu),du)
  for i in p.uidx
    resid[i] = uprev[i] - p.u_old[i] - p.dt*du[i]
  end
end

@inline function initialize!(integrator,cache::ImplicitEulerCache,f=integrator.f)
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  f(integrator.t,integrator.uprev,integrator.fsalfirst)
  integrator.kshortsize = 2
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@inline function perform_step!(integrator,cache::ImplicitEulerCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  uidx = eachindex(integrator.uprev)
  @unpack u_old,dual_cache,k,adf,rhs,uhold = cache
  copy!(u_old,uhold)
  if integrator.iter > 1 && !integrator.u_modified
    current_extrapolant!(uhold,t+dt,integrator)
  end # else uhold is previous value.
  rhs.t = t
  rhs.dt = dt
  rhs.uidx = uidx
  rhs.sizeu = size(u)
  if alg_autodiff(integrator.alg)
    nlres = NLsolve.nlsolve(adf,uhold)
  else
    nlres = NLsolve.nlsolve(rhs,uhold,autodiff=alg_autodiff(integrator.alg))
  end
  copy!(uhold,nlres.zero)
  f(t+dt,u,k)
  @pack integrator = t,dt,u
end

type RHS_Trap{F,uType,rateType,tType,SizeType,DiffCacheType,uidxType} <: Function
  f::F
  u_old::uType
  f_old::rateType
  t::tType
  dt::tType
  sizeu::SizeType
  dual_cache::DiffCacheType
  uidx::uidxType
end

function (p::RHS_Trap)(uprev,resid)
  du1 = get_du(p.dual_cache, eltype(uprev))
  p.f(p.t+p.dt,reshape(uprev,p.sizeu),du1)
  for i in p.uidx
    resid[i] = @muladd uprev[i] - p.u_old[i] - (p.dt/2)*(du1[i]+p.f_old[i])
  end
end

@inline function initialize!(integrator,cache::TrapezoidCache,f=integrator.f)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = cache.k
  f(integrator.t,integrator.uprev,integrator.fsalfirst)
  integrator.kshortsize = 2
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@inline function perform_step!(integrator,cache::TrapezoidCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  uidx = eachindex(integrator.uprev)
  @unpack u_old,dual_cache,k,rhs,adf,uhold = cache
  copy!(u_old,uhold)
  if integrator.iter > 1 && !integrator.u_modified
    current_extrapolant!(uhold,t+dt,integrator)
  end # else uhold is previous value.
  # copy!(rhs.fsalfirst,fsalfirst) Implicitly done by pointers: fsalfirst === fsalfirst == rhs.fsalfirst
  rhs.t = t
  rhs.dt = dt
  rhs.uidx = uidx
  rhs.sizeu = size(u)
  if alg_autodiff(integrator.alg)
    nlres = NLsolve.nlsolve(adf,uhold)
  else
    nlres = NLsolve.nlsolve(rhs,uhold,autodiff=alg_autodiff(integrator.alg))
  end
  copy!(uhold,nlres.zero)
  f(t+dt,u,k)
  @pack integrator = t,dt,u
end

type RHS_Trap_Scalar{F,uType,rateType,tType} <: Function
  f::F
  u_old::uType
  f_old::rateType
  t::tType
  dt::tType
end

function (p::RHS_Trap_Scalar)(uprev,resid)
  resid[1] = @muladd uprev[1] - p.u_old[1] - (p.dt/2)*(p.f_old + p.f(p.t+p.dt,uprev[1])[1])
end

@inline function initialize!(integrator,cache::TrapezoidConstantCache,f=integrator.f)
  cache.uhold[1] = integrator.uprev; cache.u_old[1] = integrator.uprev
  integrator.fsalfirst = f(integrator.t,integrator.uprev)
  integrator.kshortsize = 2
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
end

@inline function perform_step!(integrator,cache::TrapezoidConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  @unpack uhold,u_old,rhs,adf = cache
  u_old[1] = uhold[1]
  if integrator.iter > 1 && !integrator.u_modified
    uhold[1] = current_extrapolant(t+dt,integrator)
  end # else uhold is previous value.
  rhs.t = t
  rhs.dt = dt
  rhs.f_old = integrator.fsalfirst
  if alg_autodiff(integrator.alg)
    nlres = NLsolve.nlsolve(adf,uhold)
  else
    nlres = NLsolve.nlsolve(rhs,uhold,autodiff=alg_autodiff(integrator.alg))
  end
  uhold[1] = nlres.zero[1]
  k = f(t+dt,uhold[1])
  integrator.fsallast = k
  u = uhold[1]
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  @pack integrator = t,dt,u
end


text:
# ---
# title: 309. Best Time to Buy and Sell Stock with Cooldown
# id: problem309
# author: Indigo
# date: 2021-06-10
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/best-time-to-buy-and-sell-stock-with-cooldown/description/>
# hidden: true
# ---
#
# Say you have an array for which the _i_ th element is the price of a given
# stock on day _i_.
#
# Design an algorithm to find the maximum profit. You may complete as many
# transactions as you like (ie, buy one and sell one share of the stock multiple
# times) with the following restrictions:
#
#   * You may not engage in multiple transactions at the same time (ie, you must sell the stock before you buy again).
#   * After you sell your stock, you cannot buy stock on next day. (ie, cooldown 1 day)
#
# **Example:**
#
#
#
#     Input: [1,2,3,0,2]
#     Output: 3
#     Explanation: transactions = [buy, sell, cooldown, buy, sell]
#
#
#
## @lc code=start
using LeetCode

function max_profit_309(prices)
    isempty(prices) && return 0
    ## f0: this day I have a stock
    ## f1: this day I sell a stock
    ## f2: this day I can buy the stock
    f0, f1, f2 = -first(prices), 0, 0
    for price in @view(prices[2:end])
        f0, f1, f2 = max(f2 - price, f0), f0 + price, max(f1, f2)
    end
    return max(f1, f2)
end
## @lc code=end


text:
using PyPlotly
using Test

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

max_stars_repo_name: isentropic/PyPlotly

max_issues_repo_name: isentropic/PyPlotly

max_forks_repo_name: isentropic/PyPlotly


text:
abstract type AbstractParserConfiguration{D} end

deptype(::Type{<:AbstractParserConfiguration{D}}) where D = D
deptype(c::AbstractParserConfiguration) = deptype(typeof(c))

for f in (:leftdeps, :rightdeps, :leftmostdep, :rightmostdep)
    @eval $f(cfg::AbstractParserConfiguration, args...) = $f(tokens(cfg), args...)
end

"""
    stacktoken(cfg, i)

Return the token at stack index `i` (starting at 1).
"""
function stacktoken(cfg, i=1)
    stk = stack(cfg)
    stklen = length(stk)
    s_index = stacklength(cfg) - i + 1
    if 1 <= s_index <= stacklength(cfg)
        a_index = stk[s_index]
        a_index == 0 ? root(deptype(cfg)) : token(cfg, a_index)
    else
        noval(deptype(cfg))
    end
end

"""
    buffertoken(cfg, i)

Return the token at buffer index `i` (starting at 1).
"""
buffertoken(cfg, i) =
    1 <= i <= bufferlength(cfg) ? token(cfg, buffer(cfg)[i]) : noval(deptype(cfg))

struct StackBufferConfiguration{T <: Dependency}
    stack::Vector{Int}
    buffer::Vector{Int}
    A::Vector{T}
end

function StackBufferConfiguration{T}(words) where T
    σ = [0]
    β = collect(1:length(words))
    A = [unk(T, id, w) for (id,w) in enumerate(words)]
    StackBufferConfiguration{T}(σ, β, A)
end

function StackBufferConfiguration{T}(gold::DependencyTree) where T
    σ = [0]
    β = collect(1:length(gold))
    A = [dep(word, head=-1) for word in gold]
    StackBufferConfiguration{T}(σ, β, A)
end

StackBufferConfiguration(gold::DependencyTree) =
    StackBufferConfiguration{eltype(gold)}(gold)

deptype(cfg::StackBufferConfiguration) = eltype(cfg.A)

stacklength(cfg) = length(cfg.stack)
bufferlength(cfg) = length(cfg.stack)

token(cfg::StackBufferConfiguration, i) = iszero(i) ? root(deptype(cfg)) :
                                i == -1   ? noval(deptype(cfg)) :
                                cfg.A[i]

tokens(cfg::StackBufferConfiguration) = cfg.A

tokens(cfg::StackBufferConfiguration, is) =
    [token(cfg, i) for i in is if 0 <= i <= length(cfg.A)]

function popstack(cfg, n=1)
    sh = [s for s in cfg.stack[end-n+1:end]]
    st = length(cfg.stack) >= n ?
        [s for s in cfg.stack[1:end-n]] :
        Int[]
    return (st, sh...)
end

function shiftbuffer(cfg, n=1)
    @assert length(cfg.buffer) >= n
    bh = [b for b in cfg.buffer[1:n]]
    be = length(cfg.buffer) > n ?
        [b for b in cfg.buffer[n+1:end]] :
        Int[]
    return (bh..., be)
end

macro stackbufconfig(T, f=:c)
    @eval begin
        $T{T}(words::AbstractVector) where T =
            $T{T}(StackBufferConfiguration{T}(words))

        $T{T}(gold::DependencyTree) where T =
            $T{T}(StackBufferConfiguration{T}(gold))

        $T(gold::DependencyTree) = $T{eltype(gold)}(gold)

        stack(cfg::$T)  = cfg.$f.stack
        buffer(cfg::$T) = cfg.$f.buffer

        stacklength(cfg::$T) = length(cfg.$f.stack)
        bufferlength(cfg::$T) = length(cfg.$f.buffer)

        popstack(cfg::$T, args...)    = popstack(cfg.$f, args...)
        shiftbuffer(cfg::$T, args...) = shiftbuffer(cfg.$f, args...)

        token(cfg::$T, args...)  = token(cfg.$f, args...)
        tokens(cfg::$T, args...) = tokens(cfg.$f, args...)

        DependencyTree(cfg::$T, args...) = DependencyTree(cfg.$f.A, check=false)

        function Base.show(io::IO, cfg::$T)
            println(io, "$(typeof(cfg))($(stack(cfg)),$(buffer(cfg)))")
            for t in tokens(cfg)
                println(join((id(t), form(t), head(t)), "\t"))
            end
        end
    end
end

# LeftArc in ArcStandard
function leftarc_reduce2(cfg, args...; kwargs...)
    # assert a head-dependent relation between the word at the top of
    # the stack and the word directly beneath it; remove the lower
    # word from the stack
    @assert length(cfg.stack) >= 2
    σ, s1, s0 = cfg.stack[1:end-2], cfg.stack[end-1], cfg.stack[end]
    A = copy(cfg.A)
    A[s1] = dep(A[s1], args...; head=s0, kwargs...)
    StackBufferConfiguration([σ ; s0], cfg.buffer, A)
end

# LeftArc in ArcEager
# LeftArc in ArcHybrid
function leftarc_reduce(cfg, args...; kwargs...)
    # Assert a head-dependent relation between the word at the front
    # of the input buffer and the word at the top of the stack; pop
    # the stack.
    s, b, A = cfg.stack[end], cfg.buffer[1], copy(cfg.A)
    if s > 0
        A[s] = dep(A[s], args...; head=b, kwargs...)
    end
    StackBufferConfiguration(cfg.stack[1:end-1], cfg.buffer, A)
end

# RightArc in ArcHybrid
# RightArc in ArcStandard
function rightarc_reduce(cfg, args...; kwargs...)
    # assert a head-dependent relation btwn the 2nd word on the stack
    # and the word on top; remove the word at the top of the stack
    # σ, s1, s0 = σs1s0(cfg)
    σ, s1, s0 = popstack(cfg, 2)
    A = copy(cfg.A)
    if s0 > 0
        A[s0] = dep(A[s0], args...; head=s1, kwargs...)
    end
    StackBufferConfiguration([σ ; s1], cfg.buffer, A)
end

# RightArc in ArcEager
function rightarc_shift(cfg, args...; kwargs...)
    # Assert a head-dependent relation between the word on the top of
    # the σ and the word at front of the input buffer; shift the
    # word at the front of the input buffer to the stack.
    (σ, s), (b, β), A = popstack(cfg), shiftbuffer(cfg), copy(cfg.A)
    A[b] = dep(A[b], args...; head=s, kwargs...)
    StackBufferConfiguration([cfg.stack ; b], β, A)
end

# Reduce in ArcEager
function Base.reduce(cfg)
    # Pop the stack.
    StackBufferConfiguration(cfg.stack[1:end-1], cfg.buffer, cfg.A)
end

# Shift in ArcEager
# Shift in ArcHybrid
# Shift in ArcStandard
function shift(cfg::StackBufferConfiguration)
    # Remove the word from the front of the input buffer and push it
    # onto the stack.
    b, β = shiftbuffer(cfg)
    StackBufferConfiguration([cfg.stack ; b], β, cfg.A)
end

==(cfg1::StackBufferConfiguration, cfg2::StackBufferConfiguration) =
    cfg1.stack == cfg2.stack && cfg1.buffer == cfg2.buffer && cfg1.A == cfg2.A


text:
module SatisfiabilityInterface

using Symbolics
using Symbolics: Variable, Sym
using Symbolics: Assignment, get_variables, operation, arguments, value, istree

# using ReversePropagation

export DiscreteVariable, ConstraintSatisfactionProblem, DiscreteCSP, SymbolicSATProblem, SATProblem
export read_cnf, SATProblem
export satisfies
export encode, solve, all_solutions

# functions to be extend:
import Base: push!
import CommonSolve: solve

include("symbolics_interface.jl")
include("discrete_variables.jl")

include("sat_problem.jl")
include("symbolic_problem.jl")
include("MOI_wrapper.jl")
include("read_cnf.jl")
include("solver.jl")
# include("encode.jl")

include("relations.jl")

end


text:
"""

    make_julia_model(problem::VLJuliaModelObject;
        intermediate_representation_dictionary::Union{Dict{String,Any},Nothing}=nothing,
        logger::Union{Nothing,SimpleLogger}=nothing) -> VLResult

This documentation is going to be awesome. Magical.
"""
function make_julia_model(problem::VLJuliaModelObject;
    intermediate_representation_dictionary::Union{Dict{String,Any},Nothing}=nothing,
    logger::Union{Nothing,SimpleLogger}=nothing)::VLResult

    # initialize -
    src_component_set = Set{NamedTuple}()
    network_component_set = Set{NamedTuple}()
    root_component_set = Set{NamedTuple}()

    try

        # check: do we have an intermediate representation?
        ir_dictionary = intermediate_representation_dictionary
        if (isnothing(intermediate_representation_dictionary) == true)

            # parse the vff document -
            result = read_model_document(problem)
            if (isa(result.value, Exception) == true)
                throw(result.value) # TODO: what is the diff between re-throw -vs- throw?
            end
            ir_dictionary = result.value
        end

        # Step 0: load the Defaults.toml file -
        path_to_model_file = problem.path_to_model_file
        path_to_defaults_file = joinpath(dirname(path_to_model_file), problem.defaults_file_name)
        default_dictionary = TOML.parsefile(path_to_defaults_file)
        ir_dictionary["default_dictionary"] = default_dictionary

        # Step 1: Generate the "distribution" files, copy to thier location -
        path_to_output_dir = problem.path_to_output_dir

        # build the problem dictionary -
        data_dictionary_component = generate_data_dictionary_program_component(problem, ir_dictionary)
        push!(src_component_set, data_dictionary_component)

        # build Kinetics.jl program component -
        kinetics_program_component = generate_kinetics_program_component(problem, ir_dictionary)
        push!(src_component_set, kinetics_program_component)

        # build Control.jl program component -
        control_program_component = generate_control_program_component(problem, ir_dictionary)
        push!(src_component_set, control_program_component)

        # build Balances.jl program component -
        balances_program_component = generate_balances_program_component(problem, ir_dictionary)
        push!(src_component_set, balances_program_component)

        # build the include program component -
        include_program_component = generate_include_program_component(problem, ir_dictionary)
        push!(root_component_set, include_program_component)

        # build the include program component -
        init_program_component = generate_init_program_component(problem, ir_dictionary)
        push!(root_component_set, init_program_component)

        # build the include program component -
        driver_program_component = generate_driver_program_component(problem, ir_dictionary)
        push!(root_component_set, driver_program_component)

        # build the stoichiometric_matrix -
        stoichiometric_matrix_component = generate_stochiometric_matrix_component(ir_dictionary)
        push!(network_component_set, stoichiometric_matrix_component)

        # build the degradation_matrix -
        degradation_matrix_component = generate_dilution_degradation_matrix_component(ir_dictionary)
        push!(network_component_set, degradation_matrix_component)

        # dump src components to disk -
        _output_path_to_src_distribution_files = joinpath(path_to_output_dir, "src")
        write_program_components_to_disk(_output_path_to_src_distribution_files, src_component_set)

        # dump root components to disk -
        write_program_components_to_disk(path_to_output_dir, root_component_set)

        # dump network components to disk -
        _output_path_to_network_distribution_files = joinpath(_output_path_to_src_distribution_files, "network")
        write_program_components_to_disk(_output_path_to_network_distribution_files, network_component_set)

        # return -
        return VLResult(nothing)
    catch error
        return VLResult(error)
    end
end


text:
module InformationSensitivityFunctions

# Write your package code here.

end