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

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

# ╔═╡ 86e1ee96-f314-11ea-03f6-0f549b79e7c9
begin
	using Pkg
	Pkg.activate(mktempdir())
end

# ╔═╡ a4937996-f314-11ea-2ff9-615c888afaa8
begin
	Pkg.add([
			"Images",
			"ImageMagick",
			"Compose",
			"ImageFiltering",
			"TestImages",
			"Statistics",
			"PlutoUI",
			"BenchmarkTools"
			])

	using Images
	using TestImages
	using ImageFiltering
	using Statistics
	using PlutoUI
	using BenchmarkTools
end

# ╔═╡ e6b6760a-f37f-11ea-3ae1-65443ef5a81a
md"_homework 2, version 2.1_"

# ╔═╡ 85cfbd10-f384-11ea-31dc-b5693630a4c5
md"""

# **Homework 2**: _Dynamic programming_
`18.S191`, fall 2020

This notebook contains _built-in, live answer checks_! In some exercises you will see a coloured box, which runs a test case on your code, and provides feedback based on the result. Simply edit the code, run it, and the check runs again.

_For MIT students:_ there will also be some additional (secret) test cases that will be run as part of the grading process, and we will look at your notebook and write comments.

Feel free to ask questions!
"""

# ╔═╡ 33e43c7c-f381-11ea-3abc-c942327456b1
# edit the code below to set your name and kerberos ID (i.e. email without @mit.edu)

student = (name = "Alison Davey", kerberos_id = "none")

# you might need to wait until all other cells in this notebook have completed running.
# scroll around the page to see what's up

# ╔═╡ ec66314e-f37f-11ea-0af4-31da0584e881
md"""

Submission by: **_$(student.name)_** ($(student.kerberos_id)@mit.edu)
"""

# ╔═╡ 938185ec-f384-11ea-21dc-b56b7469f798
md"_Let's create a package environment:_"

# ╔═╡ 0d144802-f319-11ea-0028-cd97a776a3d0
#img = load(download("https://upload.wikimedia.org/wikipedia/commons/thumb/a/a4/Piet_Mondriaan%2C_1930_-_Mondrian_Composition_II_in_Red%2C_Blue%2C_and_Yellow.jpg/300px-Piet_Mondriaan%2C_1930_-_Mondrian_Composition_II_in_Red%2C_Blue%2C_and_Yellow.jpg"))
#img = load(download("https://upload.wikimedia.org/wikipedia/commons/thumb/5/5b/Hilma_af_Klint_-_Group_IX_SUW%2C_The_Swan_No._1_%2813947%29.jpg/477px-Hilma_af_Klint_-_Group_IX_SUW%2C_The_Swan_No._1_%2813947%29.jpg"))
img = load(download("https://i.imgur.com/4SRnmkj.png"))

# ╔═╡ cc9fcdae-f314-11ea-1b9a-1f68b792f005
md"""
# Arrays: Slices and views

In the lecture (included below) we learned about what array views are. In this exercise we will add to that understanding and look at an important use of `view`s: to reduce the amount of memory allocations when reading sub-sequences within an array.

We will use the `BenchmarkTools` package to emperically understand the effects of using views.
"""

# ╔═╡ b49a21a6-f381-11ea-1a98-7f144c55c9b7
html"""

"""

# ╔═╡ b49e8cc8-f381-11ea-1056-91668ac6ae4e
md"""
## Shrinking an array

Below is a function called `remove_in_each_row(img, pixels)`. It takes a matrix `img` and a vector of integers, `pixels`, and shrinks the image by 1 pixel in width by removing the element `img[i, pixels[i]]` in every row. This function is one of the building blocks of the Image Seam algorithm we saw in the lecture.

Read it and convince yourself that it is correct.
"""

# ╔═╡ e799be82-f317-11ea-3ae4-6d13ece3fe10
function remove_in_each_row(img, column_numbers)
	@assert size(img, 1) == length(column_numbers) # same as the number of rows
	m, n = size(img)
	local img′ = similar(img, m, n-1) # create a similar image with one less column

	# The prime (′) in the variable name is written as \prime<TAB>
    # You cannot use apostrophe for this! (Apostrophe means the transpose of a matrix)

	for (i, j) in enumerate(column_numbers)
		img′[i, :] = vcat(img[i, 1:j-1], img[i, j+1:end])
	end
	img′
end

# ╔═╡ c075a8e6-f382-11ea-2263-cd9507324f4f
md"Let's use it to remove the pixels on the diagonal. These are the image dimensions before and after doing so:"

# ╔═╡ 9cced1a8-f326-11ea-0759-0b2f22e5a1db
(before=size(img), after=size(remove_in_each_row(img, 1:size(img, 1))))

# ╔═╡ 1d893998-f366-11ea-0828-512de0c44915
md"""
## **Exercise 1** - _Making it efficient_

We can use the `@benchmark` macro from the [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl) package to benchmark fragments of Julia code.

`@benchmark` takes an expression and runs it a number of times to obtain statistics about the run time and memory allocation. We generally take the minimum time as the most stable measurement of performance ([for reasons discussed in the paper on BenchmarkTools](http://www-math.mit.edu/~edelman/publications/robust_benchmarking.pdf))
"""

# ╔═╡ 59991872-f366-11ea-1036-afe313fb4ec1
md"""
First, as an example, let's benchmark the `remove_in_each_row` function we defined above by passing in our image and a some indices to remove.
"""

# ╔═╡ e501ea28-f326-11ea-252a-53949fd9ef57
performance_experiment_default = @benchmark remove_in_each_row(img, 1:size(img, 1))

# ╔═╡ f7915918-f366-11ea-2c46-2f4671ae8a22
md"""
#### Exercise 1.1

`vcat(x, y)` is used in julia to concatenate two arrays vertically. This actually creates a new array of size `length(x) + length(y)` and copies `x` and `y` into it.  We use it in `remove_in_each_row` to create rows which have one pixel less.

While using `vcat` might make it easy to write the first version of our function, it's strictly not necessary.

👉 In `remove_in_each_row_no_vcat` below, figure out a way to avoid the use of `vcat` and modify the function to avoid it.
"""

# ╔═╡ 37d4ea5c-f327-11ea-2cc5-e3774c232c2b
function remove_in_each_row_no_vcat(img, column_numbers)
	@assert size(img, 1) == length(column_numbers) # same as the number of rows
	m, n = size(img)
	local img′ = similar(img, m, n-1) # create a similar image with one less column

	for (i, j) in enumerate(column_numbers)
		img′[i, 1:j-1] .= img[i, 1:j-1]
		img′[i, j:end] .= img[i, j+1:end]
	end
	img′
end

# ╔═╡ 67717d02-f327-11ea-0988-bfe661f57f77
performance_experiment_without_vcat = @benchmark remove_in_each_row_no_vcat(img, 1:size(img, 1))

# ╔═╡ 9e149cd2-f367-11ea-28ef-b9533e8a77bb
md"""
If you did it correctly, you should see that this benchmark shows the function running faster! And "memory estimate" should also show a smaller number, and so should "allocs estimate" which is the number of allocations done per call.
"""

# ╔═╡ ba1619d4-f389-11ea-2b3f-fd9ba71cf7e3
md"""
#### Exercise 1.2

👉 How many estimated allocations did this optimization reduce, and how can you explain most of them?
"""

# ╔═╡ e49235a4-f367-11ea-3913-f54a4a6b2d6b
no_vcat_observation = md"""
before 1029 after 687 due to vcat
"""

# ╔═╡ 837c43a4-f368-11ea-00a3-990a45cb0cbd
md"""

#### Exercise 1.3 - `view`-based optimization

👉 In the below `remove_in_each_row_views` function, implement the same optimization to remove `vcat` and use `@view` or `@views` to avoid creating copies or slices of the `img` array.

Pluto will automatically time your change with `@benchmark` below.
"""

# ╔═╡ 90a22cc6-f327-11ea-1484-7fda90283797
function remove_in_each_row_views(img, column_numbers)
	@assert size(img, 1) == length(column_numbers) # same as the number of rows
	m, n = size(img)
	local img′ = similar(img, m, n-1) # create a similar image with one less column

	for (i, j) in enumerate(column_numbers)
		img′[i, 1:j-1] .= @view img[i, 1:j-1]
		img′[i, j:end] .= @view img[i, j+1:end]
	end
	img′
end

# ╔═╡ 3335e07c-f328-11ea-0e6c-8d38c8c0ad5b
performance_experiment_views = @benchmark begin
	remove_in_each_row_views(img, 1:size(img, 1))
end

# ╔═╡ 40d6f562-f329-11ea-2ee4-d7806a16ede3
md"Final tally:"

# ╔═╡ 4f0975d8-f329-11ea-3d10-59a503f8d6b2
(
	default = performance_experiment_default,
	without_vcat = performance_experiment_without_vcat,
	views = performance_experiment_views,
)

# ╔═╡ dc63d32a-f387-11ea-37e2-6f3666a72e03
⧀(a, b) = minimum(a).allocs + size(img, 1) ÷ 2  < minimum(b).allocs;

# ╔═╡ 7eaa57d2-f368-11ea-1a70-c7c7e54bd0b1
md"""

#### Exercise 1.4

Nice! If you did your optimizations right, you should be able to get down the estimated allocations to a single digit number!

👉 How many allocations were avoided by adding the `@view` optimization over the `vcat` optimization? Why is this?
"""

# ╔═╡ fd819dac-f368-11ea-33bb-17148387546a
views_observation = md"""
before 687 after 3 no copying
"""

# ╔═╡ 318a2256-f369-11ea-23a9-2f74c566549b
md"""
## _Brightness and Energy_
"""

# ╔═╡ 7a44ba52-f318-11ea-0406-4731c80c1007
md"""
First, we will define a `brightness` function for a pixel (a color) as the mean of the red, green and blue values.

You should use this function whenever the problem set asks you to deal with _brightness_ of a pixel.
"""

# ╔═╡ 6c7e4b54-f318-11ea-2055-d9f9c0199341
begin
	brightness(c::RGB) = mean((c.r, c.g, c.b))
	brightness(c::RGBA) = mean((c.r, c.g, c.b))
end

# Use a weighted sum of rgb giving more weight to colors we perceive as 'brighter'
# Based on https://www.tutorialspoint.com/dip/grayscale_to_rgb_conversion.htm
#begin
#	brightness(c::AbstractRGB) = 0.3 * c.r + 0.59 * c.g + 0.11 * c.b
#	brightness(c::AbstractRGBA) = 0.3 * c.r + 0.59 * c.g + 0.11 * c.b
#end

# ╔═╡ 74059d04-f319-11ea-29b4-85f5f8f5c610
Gray.(brightness.(img))

# ╔═╡ 0b9ead92-f318-11ea-3744-37150d649d43
md"""We provide you with a convolve function below.
"""

# ╔═╡ d184e9cc-f318-11ea-1a1e-994ab1330c1a
convolve(img, k) = imfilter(img, reflect(k)) # uses ImageFiltering.jl package
# behaves the same way as the `convolve` function used in Lecture 2
# You were asked to implement this in homework 1.

# ╔═╡ cdfb3508-f319-11ea-1486-c5c58a0b9177
float_to_color(x) = RGB(max(0, -x), max(0, x), 0)

# ╔═╡ 5fccc7cc-f369-11ea-3b9e-2f0eca7f0f0e
md"""
finally we define the `energy` function which takes the Sobel gradients along x and y directions and computes the norm of the gradient for each pixel.
"""

# ╔═╡ 6f37b34c-f31a-11ea-2909-4f2079bf66ec
begin
	energy(∇x, ∇y) = sqrt.(∇x.^2 .+ ∇y.^2)
	function energy(img)
		∇y = convolve(brightness.(img), Kernel.sobel()[1])
		∇x = convolve(brightness.(img), Kernel.sobel()[2])
		energy(∇x, ∇y)
	end
end

# ╔═╡ 9fa0cd3a-f3e1-11ea-2f7e-bd73b8e3f302
float_to_color.(energy(img))

# ╔═╡ 87afabf8-f317-11ea-3cb3-29dced8e265a
md"""
## **Exercise 2** - _Building up to dynamic programming_

In this exercise and the following ones, we will use the computational problem of Seam carving. We will think through all the "gut reaction" solutions, and then finally end up with the dynamic programming solution that we saw in the lecture.

In the process we will understand the performance and accuracy of each iteration of our solution.

### How to implement the solutions:

For every variation of the algorithm, your job is to write a function which takes a matrix of energies, and an index for a pixel on the first row, and computes a "seam" starting at that pixel.

The function should return a vector of as many integers as there are rows in the input matrix where each number points out a pixel to delete from the corresponding row. (it acts as the input to `remove_in_each_row`).
"""

# ╔═╡ 8ba9f5fc-f31b-11ea-00fe-79ecece09c25
md"""
#### Exercise 2.1 - _The greedy approach_

The first approach discussed in the lecture (included below) is the _greedy approach_: you start from your top pixel, and at each step you just look at the three neighbors below. The next pixel in the seam is the neighbor with the lowest energy.

"""

# ╔═╡ f5a74dfc-f388-11ea-2577-b543d31576c6
html"""

"""

# ╔═╡ c3543ea4-f393-11ea-39c8-37747f113b96
md"""
👉 Implement the greedy approach.
"""

# ╔═╡ 2f9cbea8-f3a1-11ea-20c6-01fd1464a592
random_seam(m, n, i) = reduce((a, b) -> [a..., clamp(last(a) + rand(-1:1), 1, n)], 1:m-1; init=[i])

# ╔═╡ abf20aa0-f31b-11ea-2548-9bea4fab4c37
function greedy_seam(energies, starting_pixel::Int)
	m, n = size(energies)
	js = fill(1, m)
	js[1] = starting_pixel
	for i=2:m
		j1, j2 = max(1, js[i-1]-1), min(js[i-1]+1, n)
		js[i] = js[i-1]-2+findmin(@view energies[i, j1:j2])[end] + 1*(js[i-1]==1)
	end
	return js
end

# ╔═╡ 5430d772-f397-11ea-2ed8-03ee06d02a22
md"Before we apply your function to our test image, let's try it out on a small matrix of energies (displayed here in grayscale), just like in the lecture snippet above (clicking on the video will take you to the right part of the video). Light pixels have high energy, dark pixels signify low energy."

# ╔═╡ f580527e-f397-11ea-055f-bb9ea8f12015
# try
# 	if length(Set(greedy_seam(greedy_test, 5))) == 1
# 		md"Right now you are seeing the placeholder function. (You haven't done the exercise yet!) This is a straight line from the starting pixel."
# 	end
# catch end

# ╔═╡ 7ddee6fc-f394-11ea-31fc-5bd665a65bef
greedy_test = Gray.(rand(Float64, (8,10)));

# ╔═╡ 6f52c1a2-f395-11ea-0c8a-138a77f03803
md"Starting pixel: $(@bind greedy_starting_pixel Slider(1:size(greedy_test, 2); show_value=true))"

# ╔═╡ 9945ae78-f395-11ea-1d78-cf6ad19606c8
md"_Let's try it on the bigger image!_"

# ╔═╡ 87efe4c2-f38d-11ea-39cc-bdfa11298317
md"Compute shrunk image: $(@bind shrink_greedy CheckBox())"

# ╔═╡ 52452d26-f36c-11ea-01a6-313114b4445d
md"""
#### Exercise 2.2 - _Recursion_

A common pattern in algorithm design is the idea of solving a problem as the combination of solutions to subproblems.

The classic example, is a [Fibonacci number](https://en.wikipedia.org/wiki/Fibonacci_number) generator.

The recursive implementation of Fibonacci looks something like this
"""

# ╔═╡ 2a98f268-f3b6-11ea-1eea-81c28256a19e
function fib(n)
    # base case (basis)
	if n == 0 || n == 1      # `||` means "or"
		return 1
	end

    # recursion (induction)
	return fib(n-1) + fib(n-2)
end

# ╔═╡ e1bb805a-f677-11ea-2372-df15fd3e58ff
fib(1),fib(2),fib(3),fib(4),fib(5),fib(6)

# ╔═╡ 32e9a944-f3b6-11ea-0e82-1dff6c2eef8d
md"""
Notice that you can call a function from within itself which may call itself and so on until a base case is reached. Then the program will combine the result from the base case up to the final result.

In the case of the Fibonacci function, we added the solutions to the subproblems `fib(n-1)`, `fib(n-2)` to produce `fib(n)`.

An analogy can be drawn to the process of mathematical induction in mathematics. And as with mathematical induction there are parts to constructing such a recursive algorithm:

- Defining a base case
- Defining an recursion i.e. finding a solution to the problem as a combination of solutions to smaller problems.

"""

# ╔═╡ 9101d5a0-f371-11ea-1c04-f3f43b96ca4a
md"""
👉 Define a `least_energy` function which returns:
1. the lowest possible total energy for a seam starting at the pixel at $(i, j)$;
2. the column to jump to on the next move (in row $i + 1$),
which is one of $j-1$, $j$ or $j+1$, up to boundary conditions.

Return these two values in a tuple.
"""

# ╔═╡ 8ec27ef8-f320-11ea-2573-c97b7b908cb7
## returns lowest possible sum energy at pixel (i, j), and the column to jump to in row i+1.
function least_energy(energies, i, j)
	# base case
	# if i == something
	#    return energies[...] # no need for recursive computation in the base case!
	# end
	#
	# induction
	# combine results from recursive calls to `least_energy`.

	m, n = size(energies)

	if i==m
		return (energies[i,j],0)
	end

	js = j .+ [-1,0,1]
	js[1], js[3] = max(1, js[1]), min(js[3], n)

	idx,(min_energy,next_loc) = reduce((x,y) -> y[2][1] < x[2][1] ? y : x,
						         enumerate([least_energy(energies, i+1, js[1]),
						  					least_energy(energies, i+1, js[2]),
						  					least_energy(energies, i+1, js[3])
										   ]))

	if next_loc==0
		next_loc = js[idx]
	end

	return (energies[i,j] + min_energy, js[idx])
end

# ╔═╡ ee667a76-f677-11ea-1975-7b8cea947c2f
function least_energy_matrix_(energies)
	m, n = size(energies)
	least_E = zeros(size(energies))
	dirs = zeros(Int, size(energies))
	least_E[end, :] .= @view energies[end, :]
	# the minimum energy on the last row is the energy itself

	# go from the last row up, finding the minimum energy
	for i in m-1:-1:1, j in 1:n
			j1, j2 = max(1, j-1), min(j+1, n)
			e, dir = findmin(@view least_E[i+1, j1:j2])
			least_E[i,j] += e + energies[i,j]
			dirs[i, j] = (-1,0,1)[dir + (j==1)]
	end
	(least_E, dirs)
end

# ╔═╡ 05a47ab0-f678-11ea-3ef4-31cadc21070b
function least_energy_(energies,i,j)
	least_E, dirs = least_energy_matrix_(energies)
	(least_E[i,j],dirs[i,j]+j)
end

# ╔═╡ a7f3d9f8-f3bb-11ea-0c1a-55bbb8408f09
md"""
This is so elegant, correct, but inefficient! If you **check this checkbox** $(@bind compute_access CheckBox()), you will see the number of accesses made to the energies array it took to compute the least energy from the pixel (1,7):
"""

# ╔═╡ 18e0fd8a-f3bc-11ea-0713-fbf74d5fa41a
md"Whoa!"

# ╔═╡ cbf29020-f3ba-11ea-2cb0-b92836f3d04b
begin
	struct AccessTrackerArray{T,N} <: AbstractArray{T, N}
		data::Array{T,N}
		accesses::Ref{Int}
	end
	track_access(x) = AccessTrackerArray(x, Ref(0))

	Base.IndexStyle(::Type{AccessTrackerArray}) = IndexLinear()

	Base.size(x::AccessTrackerArray) = size(x.data)
	Base.getindex(x::AccessTrackerArray, i::Int...) = (x.accesses[] += 1; x.data[i...])
	Base.setindex!(x::AccessTrackerArray, v, i...) = (x.accesses[] += 1; x.data[i...] = v;)
end

# ╔═╡ 8bc930f0-f372-11ea-06cb-79ced2834720
md"""
#### Exercise 2.3 - _Exhaustive search with recursion_

Now use the `least_energy` function you defined above to define the `recursive_seam` function which takes the energies matrix and a starting pixel, and computes the seam with the lowest energy from that starting pixel.

This will give you the method used in the lecture to perform [exhaustive search of all possible paths](https://youtu.be/rpB6zQNsbQU?t=839).
"""

# ╔═╡ 85033040-f372-11ea-2c31-bb3147de3c0d
function recursive_seam(energies, starting_pixel)
	m = size(energies)[1]
	js = zeros(Int,m)
	js[1]=starting_pixel
	for k in 1:m-1
		_, js[k+1] = least_energy(energies,k,js[k])
	end
	js
end

# ╔═╡ 1d55333c-f393-11ea-229a-5b1e9cabea6a
md"Compute shrunk image: $(@bind shrink_recursive CheckBox())"

# ╔═╡ c572f6ce-f372-11ea-3c9a-e3a21384edca
md"""
#### Exercise 2.4

- State clearly why this algorithm does an exhaustive search of all possible paths.
- How does the number of possible seam grow as n increases for a `m×n` image? (Big O notation is fine, or an approximation is fine).
"""

# ╔═╡ 6d993a5c-f373-11ea-0dde-c94e3bbd1552
exhaustive_observation = md"""
It's an exhaustive search because `least_energy` is called at every pixel in the path for each loop, having previously already calculated min energy at the same pixels.
\# of seam will grow by a factor of `n`
"""

# ╔═╡ ea417c2a-f373-11ea-3bb0-b1b5754f2fac
md"""
## **Exercise 3** - _Memoization_

**Memoization** is the name given to the technique of storing results to expensive function calls that will be accessed more than once.

As stated in the video, the function `least_energy` is called repeatedly with the same arguments. In fact, we call it on the order of $3^n$ times, when there are only really $m \times n$ unique ways to call it!

Lets implement memoization on this function with first a [dictionary](https://docs.julialang.org/en/v1/base/collections/#Dictionaries) for storage.
"""

# ╔═╡ 56a7f954-f374-11ea-0391-f79b75195f4d
md"""
#### Exercise 3.1 - _Dictionary as storage_

Let's make a memoized version of least_energy function which takes a dictionary and
first checks to see if the dictionary contains the key (i,j) if it does, returns the value stored in that place, if not, will compute it, and store it in the dictionary at key (i, j) and return the value it computed.

`memoized_least_energy(energies, starting_pixel, memory)`

This function must recursively call itself, and pass the same `memory` object it received as an argument.

You are expected to read and understand the [documentation on dictionaries](https://docs.julialang.org/en/v1/base/collections/#Dictionaries) to find out how to:

1. Create a dictionary
2. Check if a key is stored in the dictionary
3. Access contents of the dictionary by a key.
"""

# ╔═╡ f61bb8e2-f679-11ea-2c07-fd74eabfe504
begin
	function f(i)
		return i*i
	end
	D=Dict(i => f(i) for i = 1:10)
end

# ╔═╡ f87af6ca-f679-11ea-1725-d109a460d06d
D[5] # value of the key 5 if it exists

# ╔═╡ fed7737a-f679-11ea-0c50-91bc97df3a5d
haskey(D, 5) # has a mapping for a given key

# ╔═╡ 087281cc-f67a-11ea-2005-5f06d24e3e2b
haskey(D, 11)

# ╔═╡ 0acd4486-f67a-11ea-0392-b1b58ea189ef
D_ = Dict('a'=>2, 'b'=>3)

# ╔═╡ 14726214-f67a-11ea-02a9-714be6fa228e
typeof(D_)

# ╔═╡ b1d09bc8-f320-11ea-26bb-0101c9a204e2
function memoized_least_energy(energies, i, j, memory)
	m, n = size(energies)
	if !haskey(memory, (i,j))
			memory[i,j]=least_energy(energies,i,j)[1]
	end
	memory[i,j]
end

# ╔═╡ 3e8b0868-f3bd-11ea-0c15-011bbd6ac051
function recursive_memoized_seam(energies, starting_pixel)
	memory = Dict{Tuple{Int,Int}, Float64}() # location => least energy.
	                                         # pass this every time you call memoized_least_energy.
	m, n = size(energies)
	js = zeros(Int,m)
	js[1]=starting_pixel
	for k in 2:m
		j1, j2 = max(1, js[k-1]-1), min(js[k-1]+1, n)
		_, dir = findmin((memoized_least_energy(energies, k, j1, memory),memoized_least_energy(energies, k, j1+1, memory),memoized_least_energy(energies, k, j2, memory)))
		js[k] = js[k-1]+(-1,0,1)[dir + (js[k-1]==1)]
	end
	js
end

# ╔═╡ 4e3bcf88-f3c5-11ea-3ada-2ff9213647b7
md"Compute shrunk image: $(@bind shrink_dict CheckBox())"

# ╔═╡ cf39fa2a-f374-11ea-0680-55817de1b837
md"""
### Exercise 3.2 - _Matrix as storage_

The dictionary-based memoization we tried above works well in general as there is no restriction on what type of keys can be used.

But in our particular case, we can use a matrix as a storage, since a matrix is naturally keyed by two integers.

Write a variation of `matrix_memoized_least_energy` and `matrix_memoized_seam` which use a matrix as storage.
"""

# ╔═╡ c8724b5e-f3bd-11ea-0034-b92af21ca12d
function matrix_memoized_least_energy(energies, i, j, memory)
	m, n = size(energies)
	if memory[i,j] == 0.
			memory[i,j]=least_energy(energies,i,j)[1]
	end
	return memory[i,j]
end

# ╔═╡ be7d40e2-f320-11ea-1b56-dff2a0a16e8d
function matrix_memoized_seam(energies, starting_pixel)
	memory = zeros(size(energies)) # use this as storage -- intially it's all zeros
	m, n = size(energies)
	js = zeros(Int,m)
	js[1]=starting_pixel

	for k in 2:m
		j1, j2 = max(1, js[k-1]-1), min(js[k-1]+1, n)
		_, dir = findmin((matrix_memoized_least_energy(energies, k, j1, memory),matrix_memoized_least_energy(energies, k, j1+1, memory),matrix_memoized_least_energy(energies, k, j2, memory)))
		js[k] = js[k-1]+(-1,0,1)[dir + (js[k-1]==1)]
	end

	js
end

# ╔═╡ 507f3870-f3c5-11ea-11f6-ada3bb087634
md"Compute shrunk image: $(@bind shrink_matrix CheckBox())"

# ╔═╡ 24792456-f37b-11ea-07b2-4f4c8caea633
md"""
## **Exercise 4** - _Dynamic programming without recursion_

Now it's easy to see that the above algorithm is equivalent to one that populates the memory matrix in a for loop.

#### Exercise 4.1

👉 Write a function which takes the energies and returns the least energy matrix which has the least possible seam energy for each pixel. This was shown in the lecture, but attempt to write it on your own.
"""

# ╔═╡ ff055726-f320-11ea-32f6-2bf38d7dd310
function least_energy_matrix(energies)
	least_E = zeros(size(energies))
	m, n = size(energies)
	least_E[end, :] .= @view energies[end, :]
	# the minimum energy on the last row is the energy itself

	# go from the last row up, finding the minimum energy
	for i in m-1:-1:1, j in 1:n
			j1, j2 = max(1, j-1), min(j+1, n)
			least_E[i,j] += findmin(@view least_E[i+1, j1:j2])[1] + energies[i,j]
	end
	least_E
end

# ╔═╡ 92e19f22-f37b-11ea-25f7-e321337e375e
md"""
#### Exercise 4.2

👉 Write a function which, when given the matrix returned by `least_energy_matrix` and a starting pixel (on the first row), computes the least energy seam from that pixel.
"""

# ╔═╡ 795eb2c4-f37b-11ea-01e1-1dbac3c80c13
function seam_from_precomputed_least_energy(energies, starting_pixel::Int)
	E = least_energy_matrix(energies)
	m, n = size(E)
	js = zeros(Int, m)
	js[1] = starting_pixel

	for i = 1:m-1
		j1, j2 = clamp(js[i]-1, 1, n), clamp(js[i]+1, 1, n)
		js[i+1] = js[i] + (-1,0,1)[(js[i]==1) + findmin(@view E[i+1, j1:j2])[2]]
	end
	js
end

# ╔═╡ 51df0c98-f3c5-11ea-25b8-af41dc182bac
md"Compute shrunk image: $(@bind shrink_bottomup CheckBox())"

# ╔═╡ 0fbe2af6-f381-11ea-2f41-23cd1cf930d9
if student.kerberos_id === "jazz"
	md"""
!!! danger "Oops!"
    **Before you submit**, remember to fill in your name and kerberos ID at the top of this notebook!
	"""
end

# ╔═╡ 6b4d6584-f3be-11ea-131d-e5bdefcc791b
md"## Function library

Just some helper functions used in the notebook."

# ╔═╡ ef88c388-f388-11ea-3828-ff4db4d1874e
function mark_path(img, path)
	img′ = copy(img)
	m = size(img, 2)
	for (i, j) in enumerate(path)
		# To make it easier to see, we'll color not just
		# the pixels of the seam, but also those adjacent to it
		for j′ in j-1:j+1
			img′[i, clamp(j′, 1, m)] = RGB(1,0,1)
		end
	end
	img′
end

# ╔═╡ 437ba6ce-f37d-11ea-1010-5f6a6e282f9b
function shrink_n(img, n, min_seam, imgs=[]; show_lightning=true)
	n==0 && return push!(imgs, img)

	e = energy(img)
	seam_energy(seam) = sum(e[i, seam[i]]  for i in 1:size(img, 1))
	_, min_j = findmin(map(j->seam_energy(min_seam(e, j)), 1:size(e, 2)))
	min_seam_vec = min_seam(e, min_j)
	img′ = remove_in_each_row(img, min_seam_vec)
	if show_lightning
		push!(imgs, mark_path(img, min_seam_vec))
	else
		push!(imgs, img′)
	end
	shrink_n(img′, n-1, min_seam, imgs)
end

# ╔═╡ f6571d86-f388-11ea-0390-05592acb9195
if shrink_greedy
	greedy_carved = shrink_n(img, 200, greedy_seam)
	md"Shrink by: $(@bind greedy_n Slider(1:200; show_value=true))"
end

# ╔═╡ f626b222-f388-11ea-0d94-1736759b5f52
if shrink_greedy
	greedy_carved[greedy_n]
end

# ╔═╡ 51e28596-f3c5-11ea-2237-2b72bbfaa001
if shrink_bottomup
	bottomup_carved = shrink_n(img, 200, seam_from_precomputed_least_energy)
	md"Shrink by: $(@bind bottomup_n Slider(1:200, show_value=true))"
end

# time this is run..1. 439 s .1: 171s

# ╔═╡ 0a10acd8-f3c6-11ea-3e2f-7530a0af8c7f
if shrink_bottomup
	bottomup_carved[bottomup_n]
end

# ╔═╡ ef26374a-f388-11ea-0b4e-67314a9a9094
function pencil(X)
	f(x) = RGB(1-x,1-x,1-x)
	map(f, X ./ maximum(X))
end

# ╔═╡ 6bdbcf4c-f321-11ea-0288-fb16ff1ec526
function decimate(img, n)
	img[1:n:end, 1:n:end]
end

# ╔═╡ ddba07dc-f3b7-11ea-353e-0f67713727fc
# Do not make this image bigger, it will be infeasible to compute.
pika = decimate(load(download("https://art.pixilart.com/901d53bcda6b27b.png")),150)

# ╔═╡ 73b52fd6-f3b9-11ea-14ed-ebfcab1ce6aa
size(pika)

# ╔═╡ 1f8d373e-f678-11ea-08ea-c9f349fa3ee0
float_to_color.(energy(pika))

# ╔═╡ fa8e2772-f3b6-11ea-30f7-699717693164
if compute_access
	tracked = track_access(energy(pika))
	least_energy(tracked, 1,7)
	tracked.accesses[]
end

# ╔═╡ d88bc272-f392-11ea-0efd-15e0e2b2cd4e
if shrink_recursive
	recursive_carved = shrink_n(pika, 3, recursive_seam)
	md"Shrink by: $(@bind recursive_n Slider(1:3, show_value=true))"
end

# ╔═╡ e66ef06a-f392-11ea-30ab-7160e7723a17
if shrink_recursive
	recursive_carved[recursive_n]
end

# ╔═╡ 90877dae-f679-11ea-2bd8-8d80ac520234
@benchmark shrink_n(pika, 3, recursive_seam)[2]

# ╔═╡ 4e3ef866-f3c5-11ea-3fb0-27d1ca9a9a3f
if shrink_dict
	dict_carved = shrink_n(pika, 3, recursive_memoized_seam)
	md"Shrink by: $(@bind dict_n Slider(1:3, show_value=true))"
end

# ╔═╡ 6e73b1da-f3c5-11ea-145f-6383effe8a89
if shrink_dict
	dict_carved[dict_n]
end

# ╔═╡ e74abf3a-f682-11ea-2350-d3a9409924a0
matrix_memoized_seam(energy(pika), 4)

# ╔═╡ 50829af6-f3c5-11ea-04a8-0535edd3b0aa
if shrink_matrix
	matrix_carved = shrink_n(pika, 3, matrix_memoized_seam)
	md"Shrink by: $(@bind matrix_n Slider(1:3, show_value=true))"
end

# ╔═╡ 9e56ecfa-f3c5-11ea-2e90-3b1839d12038
if shrink_matrix
	matrix_carved[matrix_n]
end

# ╔═╡ 4eba6b64-f69e-11ea-2afb-cfb15e56648c
least_energy_matrix(energy(pika))

# ╔═╡ e0b5b176-f681-11ea-22f4-4fcb702dd94f
seam_from_precomputed_least_energy(energy(pika), 4)

# ╔═╡ 3a66fcda-f67f-11ea-28e7-932f646316a0
@benchmark seam_from_precomputed_least_energy(energy(pika), 4)

# ╔═╡ ffc17f40-f380-11ea-30ee-0fe8563c0eb1
hint(text) = Markdown.MD(Markdown.Admonition("hint", "Hint", [text]))

# ╔═╡ 9f18efe2-f38e-11ea-0871-6d7760d0b2f6
hint(md"You can call the `least_energy` function recursively within itself to obtain the least energy of the adjacent cells and add the energy at the current cell to get the total energy.")

# ╔═╡ ffc40ab2-f380-11ea-2136-63542ff0f386
almost(text) = Markdown.MD(Markdown.Admonition("warning", "Almost there!", [text]))

# ╔═╡ ffceaed6-f380-11ea-3c63-8132d270b83f
still_missing(text=md"Replace `missing` with your answer.") = Markdown.MD(Markdown.Admonition("warning", "Here we go!", [text]))

# ╔═╡ ffde44ae-f380-11ea-29fb-2dfcc9cda8b4
keep_working(text=md"The answer is not quite right.") = Markdown.MD(Markdown.Admonition("danger", "Keep working on it!", [text]))

# ╔═╡ 980b1104-f394-11ea-0948-21002f26ee25
function visualize_seam_algorithm(algorithm, test_img, starting_pixel)
	seam = algorithm(test_img, starting_pixel)

	display_img = RGB.(test_img)
	for (i, j) in enumerate(seam)
		try
			display_img[i, j] = RGB(0.9, 0.3, 0.6)
		catch ex
			if ex isa BoundsError
				return keep_working("")
			end
			# the solution might give an illegal index
		end
	end
	display_img
end;

# ╔═╡ 2a7e49b8-f395-11ea-0058-013e51baa554
visualize_seam_algorithm(greedy_seam, greedy_test, greedy_starting_pixel)

# ╔═╡ ffe326e0-f380-11ea-3619-61dd0592d409
yays = [md"Great!", md"Yay ❤", md"Great! 🎉", md"Well done!", md"Keep it up!", md"Good job!", md"Awesome!", md"You got the right answer!", md"Let's move on to the next section."]

# ╔═╡ fff5aedc-f380-11ea-2a08-99c230f8fa32
correct(text=rand(yays)) = Markdown.MD(Markdown.Admonition("correct", "Got it!", [text]))

# ╔═╡ e3519118-f387-11ea-0c61-e1c2de1c24c1
if performance_experiment_without_vcat ⧀ performance_experiment_default
	correct()
else
	keep_working(md"We are still using (roughly) the same number of allocations as the default implementation.")
end

# ╔═╡ d4ea4222-f388-11ea-3c8d-db0d651f5282
if performance_experiment_views ⧀ performance_experiment_default
	if minimum(performance_experiment_views).allocs < 10
		correct()
	else
		keep_working(md"We are still using (roughly) the same number of allocations as the implementation without `vcat`.")
	end
else
	keep_working(md"We are still using (roughly) the same number of allocations as the default implementation.")
end

# ╔═╡ 00026442-f381-11ea-2b41-bde1fff66011
not_defined(variable_name) = Markdown.MD(Markdown.Admonition("danger", "Oopsie!", [md"Make sure that you define a variable called **$(Markdown.Code(string(variable_name)))**"]))

# ╔═╡ 145c0f58-f384-11ea-2b71-09ae83f66da2
if !@isdefined(views_observation)
	not_defined(:views_observation)
end

# ╔═╡ d7a9c000-f383-11ea-1516-cf71102d8e94
if !@isdefined(views_observation)
	not_defined(:views_observation)
end

# ╔═╡ e0622780-f3b4-11ea-1f44-59fb9c5d2ebd
if !@isdefined(least_energy_matrix)
	not_defined(:least_energy_matrix)
end

# ╔═╡ 946b69a0-f3a2-11ea-2670-819a5dafe891
if !@isdefined(seam_from_precomputed_least_energy)
	not_defined(:seam_from_precomputed_least_energy)
end

# ╔═╡ fbf6b0fa-f3e0-11ea-2009-573a218e2460
function hbox(x, y, gap=16; sy=size(y), sx=size(x))
	w,h = (max(sx[1], sy[1]),
		   gap + sx[2] + sy[2])

	slate = fill(RGB(1,1,1), w,h)
	slate[1:size(x,1), 1:size(x,2)] .= RGB.(x)
	slate[1:size(y,1), size(x,2) + gap .+ (1:size(y,2))] .= RGB.(y)
	slate
end

# ╔═╡ f010933c-f318-11ea-22c5-4d2e64cd9629
begin
	hbox(
		float_to_color.(convolve(brightness.(img), Kernel.sobel()[1])),
		float_to_color.(convolve(brightness.(img), Kernel.sobel()[2])))
end

# ╔═╡ 256edf66-f3e1-11ea-206e-4f9b4f6d3a3d
vbox(x,y, gap=16) = hbox(x', y')'

# ╔═╡ 00115b6e-f381-11ea-0bc6-61ca119cb628
bigbreak = html"

";

# ╔═╡ c086bd1e-f384-11ea-3b26-2da9e24360ca
bigbreak

# ╔═╡ 8d558c4c-f328-11ea-0055-730ead5d5c34
bigbreak

# ╔═╡ f7eba2b6-f388-11ea-06ad-0b861c764d61
bigbreak

# ╔═╡ 4f48c8b8-f39d-11ea-25d2-1fab031a514f
bigbreak

# ╔═╡ 48089a00-f321-11ea-1479-e74ba71df067
bigbreak

# ╔═╡ Cell order:
# ╟─e6b6760a-f37f-11ea-3ae1-65443ef5a81a
# ╟─ec66314e-f37f-11ea-0af4-31da0584e881
# ╟─85cfbd10-f384-11ea-31dc-b5693630a4c5
# ╠═33e43c7c-f381-11ea-3abc-c942327456b1
# ╟─938185ec-f384-11ea-21dc-b56b7469f798
# ╟─86e1ee96-f314-11ea-03f6-0f549b79e7c9
# ╟─a4937996-f314-11ea-2ff9-615c888afaa8
# ╟─0d144802-f319-11ea-0028-cd97a776a3d0
# ╟─cc9fcdae-f314-11ea-1b9a-1f68b792f005
# ╟─b49a21a6-f381-11ea-1a98-7f144c55c9b7
# ╟─b49e8cc8-f381-11ea-1056-91668ac6ae4e
# ╠═e799be82-f317-11ea-3ae4-6d13ece3fe10
# ╟─c075a8e6-f382-11ea-2263-cd9507324f4f
# ╠═9cced1a8-f326-11ea-0759-0b2f22e5a1db
# ╟─c086bd1e-f384-11ea-3b26-2da9e24360ca
# ╟─1d893998-f366-11ea-0828-512de0c44915
# ╟─59991872-f366-11ea-1036-afe313fb4ec1
# ╠═e501ea28-f326-11ea-252a-53949fd9ef57
# ╟─f7915918-f366-11ea-2c46-2f4671ae8a22
# ╠═37d4ea5c-f327-11ea-2cc5-e3774c232c2b
# ╠═67717d02-f327-11ea-0988-bfe661f57f77
# ╟─9e149cd2-f367-11ea-28ef-b9533e8a77bb
# ╟─e3519118-f387-11ea-0c61-e1c2de1c24c1
# ╟─ba1619d4-f389-11ea-2b3f-fd9ba71cf7e3
# ╠═e49235a4-f367-11ea-3913-f54a4a6b2d6b
# ╟─145c0f58-f384-11ea-2b71-09ae83f66da2
# ╟─837c43a4-f368-11ea-00a3-990a45cb0cbd
# ╠═90a22cc6-f327-11ea-1484-7fda90283797
# ╠═3335e07c-f328-11ea-0e6c-8d38c8c0ad5b
# ╟─d4ea4222-f388-11ea-3c8d-db0d651f5282
# ╟─40d6f562-f329-11ea-2ee4-d7806a16ede3
# ╟─4f0975d8-f329-11ea-3d10-59a503f8d6b2
# ╟─dc63d32a-f387-11ea-37e2-6f3666a72e03
# ╟─7eaa57d2-f368-11ea-1a70-c7c7e54bd0b1
# ╠═fd819dac-f368-11ea-33bb-17148387546a
# ╟─d7a9c000-f383-11ea-1516-cf71102d8e94
# ╟─8d558c4c-f328-11ea-0055-730ead5d5c34
# ╟─318a2256-f369-11ea-23a9-2f74c566549b
# ╟─7a44ba52-f318-11ea-0406-4731c80c1007
# ╠═6c7e4b54-f318-11ea-2055-d9f9c0199341
# ╠═74059d04-f319-11ea-29b4-85f5f8f5c610
# ╟─0b9ead92-f318-11ea-3744-37150d649d43
# ╠═d184e9cc-f318-11ea-1a1e-994ab1330c1a
# ╠═cdfb3508-f319-11ea-1486-c5c58a0b9177
# ╠═f010933c-f318-11ea-22c5-4d2e64cd9629
# ╟─5fccc7cc-f369-11ea-3b9e-2f0eca7f0f0e
# ╠═6f37b34c-f31a-11ea-2909-4f2079bf66ec
# ╠═9fa0cd3a-f3e1-11ea-2f7e-bd73b8e3f302
# ╟─f7eba2b6-f388-11ea-06ad-0b861c764d61
# ╟─87afabf8-f317-11ea-3cb3-29dced8e265a
# ╟─8ba9f5fc-f31b-11ea-00fe-79ecece09c25
# ╟─f5a74dfc-f388-11ea-2577-b543d31576c6
# ╟─c3543ea4-f393-11ea-39c8-37747f113b96
# ╟─2f9cbea8-f3a1-11ea-20c6-01fd1464a592
# ╠═abf20aa0-f31b-11ea-2548-9bea4fab4c37
# ╟─5430d772-f397-11ea-2ed8-03ee06d02a22
# ╟─f580527e-f397-11ea-055f-bb9ea8f12015
# ╟─6f52c1a2-f395-11ea-0c8a-138a77f03803
# ╟─2a7e49b8-f395-11ea-0058-013e51baa554
# ╟─7ddee6fc-f394-11ea-31fc-5bd665a65bef
# ╟─980b1104-f394-11ea-0948-21002f26ee25
# ╟─9945ae78-f395-11ea-1d78-cf6ad19606c8
# ╟─87efe4c2-f38d-11ea-39cc-bdfa11298317
# ╟─f6571d86-f388-11ea-0390-05592acb9195
# ╟─f626b222-f388-11ea-0d94-1736759b5f52
# ╟─52452d26-f36c-11ea-01a6-313114b4445d
# ╠═2a98f268-f3b6-11ea-1eea-81c28256a19e
# ╠═e1bb805a-f677-11ea-2372-df15fd3e58ff
# ╟─32e9a944-f3b6-11ea-0e82-1dff6c2eef8d
# ╟─9101d5a0-f371-11ea-1c04-f3f43b96ca4a
# ╠═ddba07dc-f3b7-11ea-353e-0f67713727fc
# ╠═73b52fd6-f3b9-11ea-14ed-ebfcab1ce6aa
# ╠═1f8d373e-f678-11ea-08ea-c9f349fa3ee0
# ╠═8ec27ef8-f320-11ea-2573-c97b7b908cb7
# ╠═05a47ab0-f678-11ea-3ef4-31cadc21070b
# ╠═ee667a76-f677-11ea-1975-7b8cea947c2f
# ╟─9f18efe2-f38e-11ea-0871-6d7760d0b2f6
# ╟─a7f3d9f8-f3bb-11ea-0c1a-55bbb8408f09
# ╟─fa8e2772-f3b6-11ea-30f7-699717693164
# ╟─18e0fd8a-f3bc-11ea-0713-fbf74d5fa41a
# ╟─cbf29020-f3ba-11ea-2cb0-b92836f3d04b
# ╟─8bc930f0-f372-11ea-06cb-79ced2834720
# ╠═85033040-f372-11ea-2c31-bb3147de3c0d
# ╟─1d55333c-f393-11ea-229a-5b1e9cabea6a
# ╟─d88bc272-f392-11ea-0efd-15e0e2b2cd4e
# ╟─e66ef06a-f392-11ea-30ab-7160e7723a17
# ╠═90877dae-f679-11ea-2bd8-8d80ac520234
# ╟─c572f6ce-f372-11ea-3c9a-e3a21384edca
# ╟─6d993a5c-f373-11ea-0dde-c94e3bbd1552
# ╟─ea417c2a-f373-11ea-3bb0-b1b5754f2fac
# ╟─56a7f954-f374-11ea-0391-f79b75195f4d
# ╠═f61bb8e2-f679-11ea-2c07-fd74eabfe504
# ╠═f87af6ca-f679-11ea-1725-d109a460d06d
# ╠═fed7737a-f679-11ea-0c50-91bc97df3a5d
# ╠═087281cc-f67a-11ea-2005-5f06d24e3e2b
# ╠═0acd4486-f67a-11ea-0392-b1b58ea189ef
# ╠═14726214-f67a-11ea-02a9-714be6fa228e
# ╠═b1d09bc8-f320-11ea-26bb-0101c9a204e2
# ╠═3e8b0868-f3bd-11ea-0c15-011bbd6ac051
# ╟─4e3bcf88-f3c5-11ea-3ada-2ff9213647b7
# ╠═4e3ef866-f3c5-11ea-3fb0-27d1ca9a9a3f
# ╟─6e73b1da-f3c5-11ea-145f-6383effe8a89
# ╟─cf39fa2a-f374-11ea-0680-55817de1b837
# ╠═c8724b5e-f3bd-11ea-0034-b92af21ca12d
# ╠═be7d40e2-f320-11ea-1b56-dff2a0a16e8d
# ╠═e74abf3a-f682-11ea-2350-d3a9409924a0
# ╟─507f3870-f3c5-11ea-11f6-ada3bb087634
# ╠═50829af6-f3c5-11ea-04a8-0535edd3b0aa
# ╟─9e56ecfa-f3c5-11ea-2e90-3b1839d12038
# ╟─4f48c8b8-f39d-11ea-25d2-1fab031a514f
# ╟─24792456-f37b-11ea-07b2-4f4c8caea633
# ╠═ff055726-f320-11ea-32f6-2bf38d7dd310
# ╠═4eba6b64-f69e-11ea-2afb-cfb15e56648c
# ╟─e0622780-f3b4-11ea-1f44-59fb9c5d2ebd
# ╟─92e19f22-f37b-11ea-25f7-e321337e375e
# ╠═795eb2c4-f37b-11ea-01e1-1dbac3c80c13
# ╠═e0b5b176-f681-11ea-22f4-4fcb702dd94f
# ╠═3a66fcda-f67f-11ea-28e7-932f646316a0
# ╟─51df0c98-f3c5-11ea-25b8-af41dc182bac
# ╠═51e28596-f3c5-11ea-2237-2b72bbfaa001
# ╠═0a10acd8-f3c6-11ea-3e2f-7530a0af8c7f
# ╠═946b69a0-f3a2-11ea-2670-819a5dafe891
# ╟─0fbe2af6-f381-11ea-2f41-23cd1cf930d9
# ╟─48089a00-f321-11ea-1479-e74ba71df067
# ╟─6b4d6584-f3be-11ea-131d-e5bdefcc791b
# ╠═437ba6ce-f37d-11ea-1010-5f6a6e282f9b
# ╠═ef88c388-f388-11ea-3828-ff4db4d1874e
# ╟─ef26374a-f388-11ea-0b4e-67314a9a9094
# ╟─6bdbcf4c-f321-11ea-0288-fb16ff1ec526
# ╟─ffc17f40-f380-11ea-30ee-0fe8563c0eb1
# ╟─ffc40ab2-f380-11ea-2136-63542ff0f386
# ╟─ffceaed6-f380-11ea-3c63-8132d270b83f
# ╟─ffde44ae-f380-11ea-29fb-2dfcc9cda8b4
# ╟─ffe326e0-f380-11ea-3619-61dd0592d409
# ╟─fff5aedc-f380-11ea-2a08-99c230f8fa32
# ╟─00026442-f381-11ea-2b41-bde1fff66011
# ╟─fbf6b0fa-f3e0-11ea-2009-573a218e2460
# ╟─256edf66-f3e1-11ea-206e-4f9b4f6d3a3d
# ╟─00115b6e-f381-11ea-0bc6-61ca119cb628


text:
"""
    flatten(a::AbstractArray)

Lorem ipsum dolor sit amet.
"""
function flatten(a::AbstractArray)
  while any(x->typeof(x)<:AbstractArray, a)
    a = collect(Iterators.flatten(a))
  end

  return a
end

ext: jl


text:
#
# TwoDimensional.jl --
#
# Utilities for 2D geometry in Julia.
#
#------------------------------------------------------------------------------
#
# This file if part of the TwoDimensional Julia package licensed under the MIT
# license (https://github.com/emmt/TwoDimensional.jl).
#
# Copyright (C) 2016-2019, Éric Thiébaut.
#
module TwoDimensional

# WeightedPoint is not exported by default because packages may implement a
# different definition.
export
    AbstractPoint,
    AffineTransform,
    BoundingBox,
    Point,
    area,
    center,
    compose,
    distance,
    exterior,
    intercept,
    interior,
    jacobian,
    rotate,
    scale,
    translate

include("basics.jl")

include("xforms.jl")
import .AffineTransforms:
    AffineTransform,
    compose,
    intercept,
    jacobian,
    rotate,
    scale,
    translate

#
# Define some aliases with a 2D suffix to avoid conflicts with other packages.
#

"""

`AffineTransform{T}` is an alias for [`TwoDimensional.AbstractPoint{T}`](@ref
AbstractPoint).

""" AbstractPoint2D
const AbstractPoint2D{T} = AbstractPoint{T}

"""

`Point2D{T}` is an alias for [`TwoDimensional.Point{T}`](@ref Point).

""" Point2D
const Point2D{T} = Point{T}

"""

`WeightedPoint2D{T}` is an alias for [`TwoDimensional.WeightedPoint{T}`](@ref
WeightedPoint).

""" WeightedPoint2D
const WeightedPoint2D{T} = WeightedPoint{T}

"""

`BoundingBox2D{T}` is an alias for [`TwoDimensional.BoundingBox{T}`](@ref
BoundingBox).

""" BoundingBox2D
const BoundingBox2D{T} = BoundingBox{T}

"""

`AffineTransform2D{T}` is an alias for
[`TwoDimensional.AffineTransform{T}`](@ref AffineTransform).

""" AffineTransform2D
const AffineTransform2D{T} = AffineTransform{T}

# The TwoDimensional.Suffixed can be imported/used instead of TwoDimensional to
# have 2D suffix for types and avoid conflicts with other packages.
module Suffixed

export
    AbstractPoint2D,
    AffineTransform2D,
    BoundingBox2D,
    Point2D,
    area,
    center,
    compose,
    distance,
    exterior,
    intercept,
    interior,
    jacobian,
    rotate,
    scale,
    translate

import ..TwoDimensional:
    AbstractPoint2D,
    AffineTransform2D,
    BoundingBox2D,
    Point2D,
    WeightedPoint2D,
    area,
    center,
    compose,
    distance,
    exterior,
    intercept,
    interior,
    jacobian,
    rotate,
    scale,
    translate

end # module Suffixed

end # module TwoDimensional


text:
using ArrayFire, Logging

function nmf(A::AbstractMatrix{T}, k::Integer, n::Int=10) where {T<:Real}
    W = rand(T, size(A,1), k)
    H = rand(T, k, size(A,2))
    for j = 1:n
        H .*= (W'*A) ./ ((W'*W)*H)
        W .*= (A*H') ./ (W*(H*H'))
    end
    return W, H
end

function nmf_gpu(A::AFArray{T}, k::Integer, n::Int=10) where {T<:Real}
    W = rand(AFArray{T}, size(A,1), k)
    H = rand(AFArray{T}, k, size(A,2))
    for j = 1:n
        H = H .* (W'*A) ./ ((W'*W)*H)
        W = W .* (A*H') ./ (W*(H*H'))
    end
    return sync(W), sync(H)
end

@info("Warmup")
a = rand(Float32,10,10)
ad = AFArray(a)
nmf(a,2)
nmf_gpu(ad,2)
for i = 1:10
    dim = i*1000
    a = rand(Float32,dim,dim)
    ad = AFArray(a)
    @info("Start Timing Dim = $dim")
    t1 = @elapsed nmf(a,2)
    t2 = @elapsed nmf_gpu(ad,2)
    println("tcpu = $t1")
    println("tgpu = $t2")
    println("Speedup = $(t1/t2)")
    GC.gc()
end
#=a = rand(1000,1000)
ad = AFArray(a)
W = rand(AFArray{Float64}, size(a,1), 2)
H = rand(AFArray{Float64}, 2, size(a,2))
W1,H1 = nmf(a, Array(W), Array(H), 2)
W2,H2 = nmf_gpu(ad, W, H, 2)
@show sum(W1-Array(W2))
@show sum(H1-Array(H2))=#


text:
export relu

doc"""
    relu(x::Var)

Rectified Linear Unit.

```math
f(x) = \max(0, x)
```
"""
relu(x::Var) = Var(relu(x.data), ∇relu!, (x,))
relu(x::T) where T<:AbstractFloat = max(x, zero(T))
relu(x::Array) = relu.(x)
relu(x::CuArray) = CUDNN.relu(x)
relu(x::Node) = Node(relu, (x,))

function ∇relu!(y::Var, x::Var)
    isnothing(x.grad) && return
    ∇relu!(y.data, y.grad, x.data, x.grad)
end

function ∇relu!(y::Array{T}, gy::Array{T}, x::Array{T}, gx::Array{T}) where T
    @inbounds for i = 1:length(x)
        gx[i] += x[i] > T(0) ? gy[i] : T(0)
    end
end

∇relu!(y::CuArray, gy::CuArray, x::CuArray, gx::CuArray) = CUDNN.∇relu!(y, gy, x, gx)


text:
include("header.jl")

@testset "rosenbrock" begin
    # info("Test rosenbrock with map...")
    b0(x) = sum(map((i, j) -> (1 - j)^2 + 10*(i - j^2)^2, x[2:end], x[1:end-1]))
    b1 = grad(b0)
    b1sum(x)=sum(b1(x))
    @test gradcheck(b0,randn(2))
    @test gradcheck(b0,(randn(2)...,))
    @test gradcheck(b0,randn(3))
    @test gradcheck(b0,(randn(3)...,))
    @test gradcheck(b1sum,randn(2))
    @test gradcheck(b1sum,(randn(2)...,))
    @test gradcheck(b1sum,randn(3))
    @test gradcheck(b1sum,(randn(3)...,))
    @time b1sum(randn(10000))
end

nothing


text:
#! format: off
get_variable_multiplier(_, ::Type{<:PSY.RenewableGen}, ::AbstractRenewableFormulation) = 1.0
get_expression_type_for_reserve(::ActivePowerReserveVariable, ::Type{<:PSY.RenewableGen}, ::Type{<:PSY.Reserve{PSY.ReserveUp}}) = ActivePowerRangeExpressionUB
get_expression_type_for_reserve(::ActivePowerReserveVariable, ::Type{<:PSY.RenewableGen}, ::Type{<:PSY.Reserve{PSY.ReserveDown}}) = ActivePowerRangeExpressionLB
########################### ActivePowerVariable, RenewableGen #################################

get_variable_binary(::ActivePowerVariable, ::Type{<:PSY.RenewableGen}, ::AbstractRenewableFormulation) = false

get_variable_lower_bound(::ActivePowerVariable, d::PSY.RenewableGen, ::AbstractRenewableFormulation) = 0.0
get_variable_upper_bound(::ActivePowerVariable, d::PSY.RenewableGen, ::AbstractRenewableFormulation) = PSY.get_max_active_power(d)

########################### ReactivePowerVariable, RenewableGen #################################

get_variable_binary(::ReactivePowerVariable, ::Type{<:PSY.RenewableGen}, ::AbstractRenewableFormulation) = false

get_multiplier_value(::TimeSeriesParameter, d::PSY.RenewableGen, ::FixedOutput) = PSY.get_max_active_power(d)
get_multiplier_value(::TimeSeriesParameter, d::PSY.RenewableGen, ::AbstractRenewableFormulation) = PSY.get_max_active_power(d)

########################Objective Function##################################################
objective_function_multiplier(::ActivePowerVariable, ::AbstractRenewableDispatchFormulation)=OBJECTIVE_FUNCTION_NEGATIVE

variable_cost(::Nothing, ::ActivePowerVariable, ::PSY.RenewableDispatch, ::AbstractRenewableDispatchFormulation)=0.0
variable_cost(cost::PSY.OperationalCost, ::ActivePowerVariable, ::PSY.RenewableDispatch, ::AbstractRenewableDispatchFormulation)=PSY.get_variable(cost)

#! format: on

get_initial_conditions_device_model(
    ::OperationModel,
    ::DeviceModel{T, <:AbstractRenewableFormulation},
) where {T <: PSY.RenewableGen} = DeviceModel(T, RenewableFullDispatch)

get_initial_conditions_device_model(
    ::OperationModel,
    ::DeviceModel{T, FixedOutput},
) where {T <: PSY.RenewableGen} = DeviceModel(T, FixedOutput)

function get_min_max_limits(
    device,
    ::Type{ReactivePowerVariableLimitsConstraint},
    ::Type{<:AbstractRenewableFormulation},
)
    return PSY.get_reactive_power_limits(device)
end

function get_default_time_series_names(
    ::Type{<:PSY.RenewableGen},
    ::Type{<:Union{FixedOutput, AbstractRenewableFormulation}},
)
    return Dict{Type{<:TimeSeriesParameter}, String}(
        ActivePowerTimeSeriesParameter => "max_active_power",
        ReactivePowerTimeSeriesParameter => "max_active_power",
    )
end

function get_default_attributes(
    ::Type{<:PSY.RenewableGen},
    ::Type{<:Union{FixedOutput, AbstractRenewableFormulation}},
)
    return Dict{String, Any}()
end

####################################### Reactive Power constraint_infos #########################

function add_constraints!(
    container::OptimizationContainer,
    T::Type{<:ReactivePowerVariableLimitsConstraint},
    U::Type{<:ReactivePowerVariable},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    X::Type{<:PM.AbstractPowerModel},
) where {V <: PSY.RenewableGen, W <: AbstractDeviceFormulation}
    add_range_constraints!(container, T, U, devices, model, X)
    return
end

"""
Reactive Power Constraints on Renewable Gen Constant power_factor
"""
function add_constraints!(
    container::OptimizationContainer,
    ::Type{<:ReactivePowerVariableLimitsConstraint},
    ::Type{<:ReactivePowerVariable},
    devices::IS.FlattenIteratorWrapper{V},
    ::DeviceModel{V, W},
    ::Type{<:PM.AbstractPowerModel},
) where {V <: PSY.RenewableGen, W <: RenewableConstantPowerFactor}
    names = [PSY.get_name(d) for d in devices]
    time_steps = get_time_steps(container)
    p_var = get_variable(container, ActivePowerVariable(), V)
    q_var = get_variable(container, ReactivePowerVariable(), V)
    jump_model = get_jump_model(container)
    constraint =
        add_constraints_container!(container, EqualityConstraint(), V, names, time_steps)
    for t in time_steps, d in devices
        name = PSY.get_name(d)
        pf = sin(acos(PSY.get_power_factor(d)))
        constraint[name, t] =
            JuMP.@constraint(jump_model, q_var[name, t] == p_var[name, t] * pf)
    end
end

function add_constraints!(
    container::OptimizationContainer,
    T::Type{ActivePowerVariableLimitsConstraint},
    U::Type{<:Union{VariableType, ExpressionType}},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    X::Type{<:PM.AbstractPowerModel},
) where {V <: PSY.RenewableGen, W <: AbstractRenewableDispatchFormulation}
    add_parameterized_upper_bound_range_constraints(
        container,
        ActivePowerVariableTimeSeriesLimitsConstraint,
        U,
        ActivePowerTimeSeriesParameter,
        devices,
        model,
        X,
    )
    return
end

##################################### renewable generation cost ############################
function objective_function!(
    container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::DeviceModel{T, U},
    ::Type{<:PM.AbstractPowerModel},
) where {T <: PSY.RenewableGen, U <: AbstractRenewableDispatchFormulation}
    add_variable_cost!(container, ActivePowerVariable(), devices, U())
    return
end


text:
@time @testset "Bars" begin

  @testset "BarInEnglish" begin

    foo = EnglishFoo()

    bar = Bar(foo,"John")

    @test greet(bar) == "Hi everybody! I am John"

  end

  @testset "BarInGerman" begin

    foo = GermanFoo()

    bar = Bar(foo,"Johanes")

    @test greet(bar) == "Guten Tag! Ich bin Johanes"

  end

end


text:
using Test
include("module-3-optional-exercise.jl")

#quantum number = 1
area_under_curve_zero_to_zero = one_dimensional_PIB_numerical_integration(1, 4, 100000, 0, 0)
area_under_curve_zero_to_L = one_dimensional_PIB_numerical_integration(1, 4, 100000, 0, 4)
area_under_curve_zero_to_half_L = one_dimensional_PIB_numerical_integration(1, 10, 100000, 0, 5)
area_under_curve_one_quarter_L_to_L = one_dimensional_PIB_numerical_integration(1, 50.0, 10000, 50.0 / 4, 50.0)

#quantum number = 2
area_under_curve_zero_to_zero_q2 = one_dimensional_PIB_numerical_integration(2, 4, 100000, 0, 0)
area_under_curve_zero_to_half_L_q2 = one_dimensional_PIB_numerical_integration(2, 10, 100000, 0, 5)
area_under_curve_zero_to_L_q2 = one_dimensional_PIB_numerical_integration(2, 4, 100000, 0, 4)

#quantum number = 3
area_under_curve_one_quarter_L_to_L_q3 = one_dimensional_PIB_numerical_integration(3, 50.0, 10000, 50.0 / 4, 50.0)
area_under_curve_one_quarter_L_to_L_q3 = one_dimensional_PIB_numerical_integration(3, 50.0, 10000, 50.0 / 4, 50.0)

@testset "$(get_module_2_excersise_3_author_name()): numerical integration tests" begin
    @testset " quantum number = 1" begin
        @test isapprox(0.0, area_under_curve_zero_to_zero; atol = 0.05)
        @test isapprox(1.0, area_under_curve_zero_to_L; atol = 0.05)
        @test isapprox(0.5, area_under_curve_zero_to_half_L; atol = 0.05)
        @test isapprox(0.90915, area_under_curve_one_quarter_L_to_L; atol = 0.05)
    end
    @testset " quantum number = 2" begin
        #quantum number = 2
        @test isapprox(0.0, area_under_curve_zero_to_zero_q2; atol = 0.05)
        @test isapprox(0.5, area_under_curve_zero_to_half_L_q2; atol = 0.05)
        @test isapprox(1.0, area_under_curve_zero_to_L_q2; atol = 0.05)
    end
    @testset " quantum number = 3" begin
        #quantum number = 3
        @test isapprox(0.69695, area_under_curve_one_quarter_L_to_L_q3; atol = 0.05)
    end
end


text:
using LinearAlgebra
using Pathfinder: WoodburyPDMat
using PDMats
using Test

include("test_utils.jl")

function decompose(W)
    n, m = size(W.B)
    k = min(n, m)
    UA = W.UA
    Q = W.Q
    UC = W.UC
    X = zeros(eltype(UC), n, n)
    X[diagind(X)] .= true
    X[1:k, 1:k] .= UC'
    T = UA' * Q * X
    return T
end

function test_decompositions(W)
    UA = cholesky(W.A).U
    Q, R = qr(UA' \ W.B)
    C = I + R * W.D * R'
    UC = cholesky(Symmetric(C)).U
    @test W.UA ≈ UA
    @test W.Q ≈ Q
    @test W.UC ≈ UC
end

@testset "WoodburyPDMat" begin
    @testset "A $Atype, D $Dtype" for T in (Float64, Float32),
        Atype in (:dense, :diag),
        Dtype in (:dense, :diag),
        n in (5, 10),
        k in (5, 10)

        A = Atype === :diag ? rand_pd_diag_mat(T, n) : rand_pd_mat(T, n)
        B = randn(T, n, k)
        D = Dtype === :diag ? rand_pd_diag_mat(T, k) : rand_pd_mat(T, k)
        W = @inferred WoodburyPDMat{T} WoodburyPDMat(A, B, D)
        Wmat = A + B * D * B'

        @testset "basic" begin
            @test eltype(W) === T
            @test eltype(Wmat) === T
            @test size(W) == (n, n)
            @test Matrix(W) ≈ Wmat
            @test W[3, 5] ≈ Wmat[3, 5]
            @test W ≈ Wmat
            test_decompositions(W)
            M = decompose(W)
            @test M * M' ≈ Wmat
            @test WoodburyPDMat(A, B, big.(D)) isa WoodburyPDMat{BigFloat}
            @test Matrix(WoodburyPDMat(A, B, big.(D))) ≈ Wmat
            Wbig = convert(AbstractMatrix{BigFloat}, W)
            @test Wbig isa WoodburyPDMat{BigFloat}
            @test Wbig ≈ Wmat
            test_decompositions(W)
            @test convert(AbstractMatrix{T}, W) === W
        end

        @testset "inv" begin
            invW = @inferred WoodburyPDMat{T} inv(W)
            @test eltype(invW) === T
            invWmat = inv(Matrix(W))
            @test invW isa WoodburyPDMat
            @test invW ≈ invWmat
            test_decompositions(invW)
        end

        @testset "determinant" begin
            @test @inferred(det(W)) ≈ det(Wmat)
            @test @inferred(logdet(W)) ≈ logdet(Wmat)
            @test @inferred(logabsdet(W))[1] ≈ logabsdet(Wmat)[1]
            @test logabsdet(W)[2] == logabsdet(Wmat)[2]
        end

        @testset "diag" begin
            @test @inferred(diag(W)) ≈ diag(Wmat)
        end

        @testset "adjoint/transpose" begin
            @test W' === W
            @test transpose(W) === W
        end

        @testset "lmul!" begin
            x = randn(T, n)
            y = Wmat * x
            @test lmul!(W, x) === x
            @test x ≈ y

            X = randn(T, n, 5)
            Y = Wmat * X
            @test lmul!(W, X) === X
            @test X ≈ Y
        end

        @testset "mul!" begin
            x = randn(T, n)
            y = similar(x)
            @test mul!(y, W, x) === y
            @test y ≈ Wmat * x

            X = randn(T, n, 5)
            Y = similar(X)
            @test mul!(Y, W, X) === Y
            @test Y ≈ Wmat * X
        end

        @testset "*" begin
            @inferred Union{WoodburyPDMat{Float64},Matrix{Float64}} W * 5.0
            @test W * 5.0 isa WoodburyPDMat
            @test W * 5.0 ≈ Wmat * 5
            test_decompositions(W * 5.0)
            test_decompositions(W * 3)
            @test W * -2 isa Matrix
            @test W * -2 ≈ Wmat * -2

            x = randn(T, n)
            @test W * x ≈ Wmat * x

            X = randn(T, n)
            @test W * X ≈ Wmat * X
        end

        @testset "PDMats.dim" begin
            @test PDMats.dim(W) == n
        end

        @testset "unwhiten" begin
            M = decompose(W)

            x = randn(T, n)
            @test @inferred(unwhiten(W, x)) ≈ M * x

            X = randn(T, n, 100)
            @test @inferred(unwhiten(W, X)) ≈ M * X
        end

        @testset "invunwhiten!" begin
            M = decompose(W)

            x = randn(T, n)
            @test @inferred(Pathfinder.invunwhiten!(similar(x), W, x)) ≈ M' \ x

            X = randn(T, n, 100)
            @test @inferred(Pathfinder.invunwhiten!(similar(X), W, X)) ≈ M' \ X
        end

        @testset "PDMats.quad" begin
            x = randn(T, n)
            @test @inferred(quad(W, x)) ≈ dot(x, Wmat, x)

            u = randn(T, n)
            @test quad(W, Pathfinder.invunwhiten!(similar(u), W, u)) ≈ dot(u, u)

            X = randn(T, n, 100)
            @test @inferred(quad(W, X)) ≈ quad(PDMats.PDMat(Symmetric(Wmat)), X)

            U = randn(T, n, 100)
            @test quad(W, Pathfinder.invunwhiten!(similar(U), W, U)) ≈
                vec(sum(abs2, U; dims=1))
        end

        @testset "PDMats.invquad" begin
            x = randn(T, n)
            @test @inferred(invquad(W, x)) ≈ dot(x, inv(Wmat), x)

            u = randn(T, n)
            @test invquad(W, unwhiten(W, u)) ≈ dot(u, u)

            X = randn(T, n, 100)
            @test @inferred(invquad(W, X)) ≈ invquad(PDMats.PDMat(Symmetric(Wmat)), X)

            U = randn(T, n, 100)
            @test invquad(W, unwhiten(W, U)) ≈ vec(sum(abs2, U; dims=1))
        end
    end
end


text:
function chopzero!{T<:Real,D}(arr ::Array{Complex{T},D};
                              tol::Real=eps(T))
  for i in eachindex(arr)
    if isapprox(real(arr[i]), 0; atol=tol)
      arr[i] = 1im * imag(arr[i])
    end
    if isapprox(imag(arr[i]), 0; atol=tol)
      arr[i] = real(arr[i])
    end
  end
  arr
end

function chopzero!{T<:Real,D}(arr ::Array{T, D};
                              tol::Real=eps(T))
  for i in eachindex(arr)
    if isapprox(arr[i], 0; atol=tol)
      arr[i] = 0
    end
  end
  arr
end

function chopzero{T<:Real,D}(arr_in ::Array{Complex{T},D};
                             tol::Real=eps(T))
    arr = copy(arr_in)
    for i in eachindex(arr)
        if isapprox(real(arr[i]), 0; atol=tol)
            arr[i] = 1im * imag(arr[i])
        end
        if isapprox(imag(arr[i]), 0; atol=tol)
            arr[i] = real(arr[i])
        end
    end
    arr
end

function chopzero{T<:Real,D}(arr_in :: Array{T, D};
                             tol::Real=eps(T))
  arr = copy(arr_in)
  for i in eachindex(arr)
    if isapprox(arr[i], 0; atol=tol)
      arr[i] = 0
    end
  end
  arr
end