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

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

# ╔═╡ 83eb9ca0-ed68-11ea-0bc5-99a09c68f867
md"_homework 1, version 4_"

# ╔═╡ ac8ff080-ed61-11ea-3650-d9df06123e1f
md"""

# **Homework 1** - _convolutions_
`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!
"""

# ╔═╡ 911ccbce-ed68-11ea-3606-0384e7580d7c
# edit the code below to set your name and kerberos ID (i.e. email without @mit.edu)

student = (name = "Debaditya Bhattacharya", kerberos_id = "debbh")

# press the ▶ button in the bottom right of this cell to run your edits
# or use Shift+Enter

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

# ╔═╡ 8ef13896-ed68-11ea-160b-3550eeabbd7d
md"""

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

# ╔═╡ 5f95e01a-ee0a-11ea-030c-9dba276aba92
md"_Let's create a package environment:_"

# ╔═╡ 65780f00-ed6b-11ea-1ecf-8b35523a7ac0
begin
	import Pkg
	Pkg.activate(mktempdir())
end

# ╔═╡ 74b008f6-ed6b-11ea-291f-b3791d6d1b35
begin
	Pkg.add(["Images", "ImageMagick"])
	using Images
end

# ╔═╡ 6b30dc38-ed6b-11ea-10f3-ab3f121bf4b8
begin
	Pkg.add("PlutoUI")
	using PlutoUI
end

# ╔═╡ 67461396-ee0a-11ea-3679-f31d46baa9b4
md"_We set up Images.jl again:_"

# ╔═╡ 540ccfcc-ee0a-11ea-15dc-4f8120063397
md"""
## **Exercise 1** - _Manipulating vectors (1D images)_

A `Vector` is a 1D array. We can think of that as a 1D image.

"""

# ╔═╡ 467856dc-eded-11ea-0f83-13d939021ef3
example_vector = [0.5, 0.4, 0.3, 0.2, 0.1, 0.0, 0.7, 0.0, 0.7, 0.9]

# ╔═╡ ad6a33b0-eded-11ea-324c-cfabfd658b56
md"#### Exerise 1.1
👉 Make a random vector `random_vect` of length 10 using the `rand` function.
"

# ╔═╡ f51333a6-eded-11ea-34e6-bfbb3a69bcb0
random_vect = rand(10) # replace this with your code!

# ╔═╡ cf738088-eded-11ea-2915-61735c2aa990
md"👉 Make a function `mean` using a `for` loop, which computes the mean/average of a vector of numbers."

# ╔═╡ 0ffa8354-edee-11ea-2883-9d5bfea4a236
function mean(x)
	val = sum(x)/ length(x)
	return val
end

# ╔═╡ 1f104ce4-ee0e-11ea-2029-1d9c817175af
mean([1, 2, 3])

# ╔═╡ 1f229ca4-edee-11ea-2c56-bb00cc6ea53c
md"👉 Define `m` to be the mean of `random_vect`."

# ╔═╡ 2a391708-edee-11ea-124e-d14698171b68
m = mean(random_vect)

# ╔═╡ e2863d4c-edef-11ea-1d67-332ddca03cc4
md"""👉 Write a function `demean`, which takes a vector `x` and subtracts the mean from each value in `x`."""

# ╔═╡ ec5efe8c-edef-11ea-2c6f-afaaeb5bc50c
function demean(x)
	return x .- mean(x)
end

# ╔═╡ 29e10640-edf0-11ea-0398-17dbf4242de3
md"Let's check that the mean of the `demean(random_vect)` is 0:

_Due to floating-point round-off error it may *not* be *exactly* 0._"

# ╔═╡ 6f67657e-ee1a-11ea-0c2f-3d567bcfa6ea
if ismissing(random_vect)
	md"""
	!!! info
	    The following cells error because `random_vect` is not yet defined. Have you done the first exercise?
	"""
end

# ╔═╡ 73ef1d50-edf0-11ea-343c-d71706874c82
copy_of_random_vect = copy(random_vect); # in case demean modifies `x`

# ╔═╡ 38155b5a-edf0-11ea-3e3f-7163da7433fb
mean(demean(copy_of_random_vect))

# ╔═╡ a5f8bafe-edf0-11ea-0da3-3330861ae43a
md"""
#### Exercise 1.2

👉 Generate a vector of 100 zeros. Change the center 20 elements to 1.
"""

# ╔═╡ b6b65b94-edf0-11ea-3686-fbff0ff53d08
function create_bar()
	x = zeros(100)
	x[40:60] .= 1
	return x
end

# ╔═╡ 22f28dae-edf2-11ea-25b5-11c369ae1253
md"""
#### Exercise 1.3

👉 Write a function that turns a `Vector` of `Vector`s into a `Matrix`.
"""

# ╔═╡ 8c19fb72-ed6c-11ea-2728-3fa9219eddc4
function vecvec_to_matrix(vecvec)
	mat = vecvec[1]
	flag = 1
	for vec in vecvec
		if flag == 1
			flag = 0
			continue
		end
		mat = [mat vec]
	end
	return mat
end

# ╔═╡ c4761a7e-edf2-11ea-1e75-118e73dadbed
vecvec_to_matrix([[1,2], [3,4]])

# ╔═╡ 40881910-4865-11eb-1ff7-3798bccfbc7b
typeof([[1,2], [3,4]])

# ╔═╡ 393667ca-edf2-11ea-09c5-c5d292d5e896
md"""

👉 Write a function that turns a `Matrix` into a`Vector` of `Vector`s .
"""

# ╔═╡ 9f1c6d04-ed6c-11ea-007b-75e7e780703d
function matrix_to_vecvec(matrix)
	vecvec = [matrix[:,1]]
	for i in 1:size(matrix)[2]
		if i == 1
			continue
		end
		vecvec = [vecvec; [matrix[:,i]]]
	end
	return vecvec
end

# ╔═╡ 70955aca-ed6e-11ea-2330-89b4d20b1795
matrix_to_vecvec([6 7; 8 9])

# ╔═╡ ab4fc2c0-4865-11eb-0a6c-db550806559a
matrix_to_vecvec([6 7 8; 8 9 10])

# ╔═╡ 5da8cbe8-eded-11ea-2e43-c5b7cc71e133
begin
	colored_line(x::Vector{<:Real}) = Gray.(Float64.((hcat(x)')))
	colored_line(x::Any) = nothing
end

# ╔═╡ 56ced344-eded-11ea-3e81-3936e9ad5777
colored_line(example_vector)

# ╔═╡ b18e2c54-edf1-11ea-0cbf-85946d64b6a2
colored_line(random_vect)

# ╔═╡ d862fb16-edf1-11ea-36ec-615d521e6bc0
colored_line(create_bar())

# ╔═╡ e083b3e8-ed61-11ea-2ec9-217820b0a1b4
md"""
## **Exercise 2** - _Manipulating images_

In this exercise we will get familiar with matrices (2D arrays) in Julia, by manipulating images.
Recall that in Julia images are matrices of `RGB` color objects.

Let's load a picture of Philip again.
"""

# ╔═╡ c5484572-ee05-11ea-0424-f37295c3072d
philip_file = download("https://i.imgur.com/VGPeJ6s.jpg")

# ╔═╡ e86ed944-ee05-11ea-3e0f-d70fc73b789c
md"_Hi there Philip_"

# ╔═╡ c54ccdea-ee05-11ea-0365-23aaf053b7d7
md"""
#### Exercise 2.1
👉 Write a function **`mean_colors`** that accepts an object called `image`. It should calculate the mean (average) amounts of red, green and blue in the image and return a tuple `(r, g, b)` of those means.
"""

# ╔═╡ f6898df6-ee07-11ea-2838-fde9bc739c11
function mean_colors(image)
	mean(image)
	return mean(image)
end

# ╔═╡ f68d4a36-ee07-11ea-0832-0360530f102e
md"""
#### Exercise 2.2
👉 Look up the documentation on the `floor` function. Use it to write a function `quantize(x::Number)` that takes in a value $x$ (which you can assume is between 0 and 1) and "quantizes" it into bins of width 0.1. For example, check that 0.267 gets mapped to 0.2.
"""

# ╔═╡ f6991a50-ee07-11ea-0bc4-1d68eb028e6a
begin
	function quantize(x::Number)
		x*=10
		return floor(x)/10
	end

	function quantize(color::AbstractRGB)
		# you will write me in a later exercise!
		return RGB(quantize(color.r),quantize(color.g),quantize(color.b))
	end

	function quantize(image::AbstractMatrix)
		# you will write me in a later exercise!

		return quantize.(image)
	end
end

# ╔═╡ f6a655f8-ee07-11ea-13b6-43ca404ddfc7
quantize(0.267), quantize(0.91)

# ╔═╡ 1fc7dd0e-4869-11eb-1761-552ef1069515
color = RGB(0.2,0.45,0.13)

# ╔═╡ ce65f880-4868-11eb-288c-ad505b2bacf8
typeof(RGB(quantize(color.r),quantize(color.g),quantize(color.b)))

# ╔═╡ f6b218c0-ee07-11ea-2adb-1968c4fd473a
md"""
#### Exercise 2.3
👉 Write the second **method** of the function `quantize`, i.e. a new *version* of the function with the *same* name. This method will accept a color object called `color`, of the type `AbstractRGB`.

_Write the function in the same cell as `quantize(x::Number)` from the last exercise. 👆_

Here, `::AbstractRGB` is a **type annotation**. This ensures that this version of the function will be chosen when passing in an object whose type is a **subtype** of the `AbstractRGB` abstract type. For example, both the `RGB` and `RGBX` types satisfy this.

The method you write should return a new `RGB` object, in which each component ($r$, $g$ and $b$) are quantized.
"""

# ╔═╡ f6bf64da-ee07-11ea-3efb-05af01b14f67
md"""
#### Exercise 2.4
👉 Write a method `quantize(image::AbstractMatrix)` that quantizes an image by quantizing each pixel in the image. (You may assume that the matrix is a matrix of color objects.)

_Write the function in the same cell as `quantize(x::Number)` from the last exercise. 👆_
"""

# ╔═╡ 25dad7ce-ee0b-11ea-3e20-5f3019dd7fa3
md"Let's apply your method!"

# ╔═╡ f6cc03a0-ee07-11ea-17d8-013991514d42
md"""
#### Exercise 2.5
👉 Write a function `invert` that inverts a color, i.e. sends $(r, g, b)$ to $(1 - r, 1-g, 1-b)$.
"""

# ╔═╡ 63e8d636-ee0b-11ea-173d-bd3327347d55
function invert(color::AbstractRGB)

	return RGB(1,1,1)-color
end

# ╔═╡ 2cc2f84e-ee0d-11ea-373b-e7ad3204bb00
md"Let's invert some colors:"

# ╔═╡ b8f26960-ee0a-11ea-05b9-3f4bc1099050
black = RGB(0.0, 0.0, 0.0)

# ╔═╡ 5de3a22e-ee0b-11ea-230f-35df4ca3c96d
invert(black)

# ╔═╡ 4e21e0c4-ee0b-11ea-3d65-b311ae3f98e9
red = RGB(0.8, 0.1, 0.1)

# ╔═╡ 6dbf67ce-ee0b-11ea-3b71-abc05a64dc43
invert(red)

# ╔═╡ 846b1330-ee0b-11ea-3579-7d90fafd7290
md"Can you invert the picture of Philip?"

# ╔═╡ f6d6c71a-ee07-11ea-2b63-d759af80707b
md"""
#### Exercise 2.6
👉 Write a function `noisify(x::Number, s)` to add randomness of intensity $s$ to a value $x$, i.e. to add a random value between $-s$ and $+s$ to $x$. If the result falls outside the range $(0, 1)$ you should "clamp" it to that range. (Note that Julia has a `clamp` function, but you should write your own function `myclamp(x)`.)
"""

# ╔═╡ 70ac1b40-486b-11eb-2b38-1920b4dcf32b
function myclamp(x::Number)
	if x >= 1
		return 1
	elseif x <= 0
		return 0
	else
		return x
	end
end

# ╔═╡ f6e2cb2a-ee07-11ea-06ee-1b77e34c1e91
begin
	function noisify(x::Number, s)
		return myclamp(x + rand((-1,1))*s)
	end

	function noisify(color::AbstractRGB, s)
		# you will write me in a later exercise!
		return RGB(noisify(color.r,s),noisify(color.g,s),noisify(color.b,s))
	end

	function noisify(image::AbstractMatrix, s)
		# you will write me in a later exercise!
		return noisify.(image,s)
	end
end

# ╔═╡ f6fc1312-ee07-11ea-39a0-299b67aee3d8
md"""
👉  Write the second method `noisify(c::AbstractRGB, s)` to add random noise of intensity $s$ to each of the $(r, g, b)$ values in a colour.

_Write the function in the same cell as `noisify(x::Number)` from the last exercise. 👆_
"""

# ╔═╡ 774b4ce6-ee1b-11ea-2b48-e38ee25fc89b
@bind color_noise Slider(0:0.01:1, show_value=true)

# ╔═╡ 7e4aeb70-ee1b-11ea-100f-1952ba66f80f
noisify(red, color_noise)

# ╔═╡ 6a05f568-ee1b-11ea-3b6c-83b6ada3680f

# ╔═╡ f70823d2-ee07-11ea-2bb3-01425212aaf9
md"""
👉 Write the third method `noisify(image::AbstractMatrix, s)` to noisify each pixel of an image.

_Write the function in the same cell as `noisify(x::Number)` from the last exercise. 👆_
"""

# ╔═╡ e70a84d4-ee0c-11ea-0640-bf78653ba102
@bind philip_noise Slider(0:0.01:2, show_value=true)

# ╔═╡ 9604bc44-ee1b-11ea-28f8-7f7af8d0cbb2

# ╔═╡ f714699e-ee07-11ea-08b6-5f5169861b57
md"""
👉 For which noise intensity does it become unrecognisable?

You may need noise intensities larger than 1. Why?

"""

# ╔═╡ bdc2df7c-ee0c-11ea-2e9f-7d2c085617c1
answer_about_noise_intensity = md"""
The image is unrecognisable with intensity ...
"""

# ╔═╡ 81510a30-ee0e-11ea-0062-8b3327428f9d

# ╔═╡ e3b03628-ee05-11ea-23b6-27c7b0210532
decimate(image, ratio=5) = image[1:ratio:end, 1:ratio:end]

# ╔═╡ c8ecfe5c-ee05-11ea-322b-4b2714898831
philip = let
	original = Images.load(philip_file)
	decimate(original, 8)
end

# ╔═╡ 5be9b144-ee0d-11ea-2a8d-8775de265a1d
mean_colors(philip)

# ╔═╡ 9751586e-ee0c-11ea-0cbb-b7eda92977c9
quantize(philip)

# ╔═╡ 943103e2-ee0b-11ea-33aa-75a8a1529931
philip_inverted = invert.(philip)

# ╔═╡ ac15e0d0-ee0c-11ea-1eaf-d7f88b5df1d7
noisify(philip, philip_noise)

# ╔═╡ e08781fa-ed61-11ea-13ae-91a49b5eb74a
md"""

## **Exercise 3** - _Convolutions_

As we have seen in the videos, we can produce cool effects using the mathematical technique of **convolutions**. We input one image $M$ and get a new image $M'$ back.

Conceptually we think of $M$ as a matrix. In practice, in Julia it will be a `Matrix` of color objects, and we may need to take that into account. Ideally, however, we should write a **generic** function that will work for any type of data contained in the matrix.

A convolution works on a small **window** of an image, i.e. a region centered around a given point $(i, j)$. We will suppose that the window is a square region with odd side length $2\ell + 1$, running from $-\ell, \ldots, 0, \ldots, \ell$.

The result of the convolution over a given window, centred at the point $(i, j)$ is a *single number*; this number is the value that we will use for $M'_{i, j}$.
(Note that neighbouring windows overlap.)

To get started let's restrict ourselves to convolutions in 1D.
So a window is just a 1D region from $-\ell$ to $\ell$.

"""

# ╔═╡ 7fc8ee1c-ee09-11ea-1382-ad21d5373308
md"""
---

Let's create a vector `v` of random numbers of length `n=100`.
"""

# ╔═╡ 7fcd6230-ee09-11ea-314f-a542d00d582e
n = 100

# ╔═╡ 7fdb34dc-ee09-11ea-366b-ffe10d1aa845
v = rand(n)

# ╔═╡ 7fe9153e-ee09-11ea-15b3-6f24fcc20734
md"_Feel free to experiment with different values!_"

# ╔═╡ 80108d80-ee09-11ea-0368-31546eb0d3cc
md"""
#### Exercise 3.1
You've seen some colored lines in this notebook to visualize arrays. Can you make another one?

👉 Try plotting our vector `v` using `colored_line(v)`.
"""

# ╔═╡ 01070e28-ee0f-11ea-1928-a7919d452bdd
colored_line(v)

# ╔═╡ 7522f81e-ee1c-11ea-35af-a17eb257ff1a
md"Try changing `n` and `v` around. Notice that you can run the cell `v = rand(n)` again to regenerate new random values."

# ╔═╡ 801d90c0-ee09-11ea-28d6-61b806de26dc
md"""
#### Exercise 3.2
We need to decide how to handle the **boundary conditions**, i.e. what happens if we try to access a position in the vector `v` beyond `1:n`.  The simplest solution is to assume that $v_{i}$ is 0 outside the original vector; however, this may lead to strange boundary effects.

A better solution is to use the *closest* value that is inside the vector. Effectively we are extending the vector and copying the extreme values into the extended positions. (Indeed, this is one way we could implement this; these extra positions are called **ghost cells**.)

👉 Write a function `extend(v, i)` that checks whether the position $i$ is inside `1:n`. If so, return the $i$th component of `v`; otherwise, return the nearest end value.
"""

# ╔═╡ 802bec56-ee09-11ea-043e-51cf1db02a34
function extend(v, i)
	if i < 1
		x = v[1]
	elseif i > length(v)
		x = v[length(v)]
	else
		x = v[i]
	end
	return x
end

# ╔═╡ b7f3994c-ee1b-11ea-211a-d144db8eafc2
md"_Some test cases:_"

# ╔═╡ 803905b2-ee09-11ea-2d52-e77ff79693b0
extend(v, 1)

# ╔═╡ 80479d98-ee09-11ea-169e-d166eef65874
extend(v, -8)

# ╔═╡ 805691ce-ee09-11ea-053d-6d2e299ee123
extend(v, n + 10)

# ╔═╡ 806e5766-ee0f-11ea-1efc-d753cd83d086
md"Extended with 0:"

# ╔═╡ 38da843a-ee0f-11ea-01df-bfa8b1317d36
colored_line([0, 0, example_vector..., 0, 0])

# ╔═╡ 9bde9f92-ee0f-11ea-27f8-ffef5fce2b3c
md"Extended with your `extend`:"

# ╔═╡ 45c4da9a-ee0f-11ea-2c5b-1f6704559137
if extend(v,1) === missing
	missing
else
	colored_line([extend(example_vector, i) for i in -1:12])
end

# ╔═╡ 80664e8c-ee09-11ea-0702-711bce271315
md"""
#### Exercise 3.3
👉 Write a function `blur_1D(v, l)` that blurs a vector `v` with a window of length `l` by averaging the elements within a window from $-\ell$ to $\ell$. This is called a **box blur**.
"""

# ╔═╡ 807e5662-ee09-11ea-3005-21fdcc36b023
function blur_1D(v, l)
	s = length(v)
	ev = [extend(v, i) for i in (1-l):(s+l)]
	return [mean(ev[i:(i+2*l)]) for i in 1:s]

	# v_blur = zeros(length(v))
	# for i in 1:length(v)
	# 	v_blur[i] = mean([extend(v,i+j) for j in -l:l])
	# end
	# return v_blur
end

# ╔═╡ 808deca8-ee09-11ea-0ee3-1586fa1ce282
let
	try
		test_v = rand(n)
		original = copy(test_v)
		blur_1D(test_v, 5)
		if test_v != original
			md"""
			!!! danger "Oopsie!"
			    It looks like your function _modifies_ `v`. Can you write it without doing so? Maybe you can use `copy`.
			"""
		end
	catch
	end
end

# ╔═╡ 809f5330-ee09-11ea-0e5b-415044b6ac1f
md"""
#### Exercise 3.4
👉 Apply the box blur to your vector `v`. Show the original and the new vector by creating two cells that call `colored_line`. Make the parameter $\ell$ interactive, and call it `l_box` instead of just `l` to avoid a variable naming conflict.
"""

# ╔═╡ 4c5cdff0-48ce-11eb-10b8-691432d0012a
@bind l_box Slider(0:1:100, show_value=true)

# ╔═╡ ca1ac5f4-ee1c-11ea-3d00-ff5268866f87
colored_line(v)

# ╔═╡ 5cde0890-48ce-11eb-03a2-fda7bcffcfec
colored_line(blur_1D(v,l_box))

# ╔═╡ 80ab64f4-ee09-11ea-29b4-498112ed0799
md"""
#### Exercise 3.5
The box blur is a simple example of a **convolution**, i.e. a linear function of a window around each point, given by

$$v'_{i} = \sum_{n}  \, v_{i - n} \, k_{n},$$

where $k$ is a vector called a **kernel**.

Again, we need to take care about what happens if $v_{i -n }$ falls off the end of the vector.

👉 Write a function `convolve_vector(v, k)` that performs this convolution. You need to think of the vector $k$ as being *centred* on the position $i$. So $n$ in the above formula runs between $-\ell$ and $\ell$, where $2\ell + 1$ is the length of the vector $k$. You will need to do the necessary manipulation of indices.
"""

# ╔═╡ 28e20950-ee0c-11ea-0e0a-b5f2e570b56e
function convolve_vector(v, k)
	s = length(v)
	l = (length(k)-1)÷2
	ev = [extend(v, i) for i in (1-l):(s+l)]
	return [sum(ev[i:(i+2*l)].*k) for i in 1:s]
end

# ╔═╡ 93284f92-ee12-11ea-0342-833b1a30625c
test_convolution = let
	v = [1, 10, 100, 1000, 10000]
	k = [0, 1, 0]
	convolve_vector(v, k)
end

# ╔═╡ 5eea882c-ee13-11ea-0d56-af81ecd30a4a
colored_line(test_convolution)

# ╔═╡ cf73f9f8-ee12-11ea-39ae-0107e9107ef5
md"_Edit the cell above, or create a new cell with your own test cases!_"

# ╔═╡ 80b7566a-ee09-11ea-3939-6fab470f9ec8
md"""
#### Exercise 3.6
👉 Write a function `gaussian_kernel`.

The definition of a Gaussian in 1D is

$$G(x) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp \left( \frac{-x^2}{2\sigma^2} \right)$$

We need to **sample** (i.e. evaluate) this at each pixel in a region of size $n^2$,
and then **normalize** so that the sum of the resulting kernel is 1.

For simplicity you can take $\sigma=1$.
"""

# ╔═╡ 1c8b4658-ee0c-11ea-2ede-9b9ed7d3125e
function gaussian_kernel(n)
	g(x) = exp(-(x^2)/2)/sqrt(2π)
	a = [g(i) for i in -(n÷2):n÷2]
	return a./sum(a)
end

# ╔═╡ 664b6270-48e1-11eb-1c33-3b4d2fa7d90d
gaussian_kernel(7)

# ╔═╡ bae5f940-48e0-11eb-2762-a38b8a33bf0d
colored_line(gaussian_kernel(7))

# ╔═╡ f8bd22b8-ee14-11ea-04aa-ab16fd01826e
md"Let's test your kernel function!"

# ╔═╡ 576f53c0-48e0-11eb-0d9e-d9e5f1e794eb
@bind kernal_size Slider(1:2:50, show_value=true)

# ╔═╡ 2a9dd06a-ee13-11ea-3f84-67bb309c77a8
gaussian_kernel_size_1D = kernal_size # change this value, or turn me into a slider!

# ╔═╡ 38eb92f6-ee13-11ea-14d7-a503ac04302e
test_gauss_1D_a = let
	v = random_vect
	k = gaussian_kernel(gaussian_kernel_size_1D)

	if k !== missing
		convolve_vector(v, k)
	end
end

# ╔═╡ b424e2aa-ee14-11ea-33fa-35491e0b9c9d
colored_line(test_gauss_1D_a)

# ╔═╡ 24c21c7c-ee14-11ea-1512-677980db1288
test_gauss_1D_b = let
	v = create_bar()
	k = gaussian_kernel(gaussian_kernel_size_1D)

	if k !== missing
		convolve_vector(v, k)
	end
end

# ╔═╡ bc1c20a4-ee14-11ea-3525-63c9fa78f089
colored_line(test_gauss_1D_b)

# ╔═╡ b01858b6-edf3-11ea-0826-938d33c19a43
md"""

## **Exercise 4** - _Convolutions of images_

Now let's move to 2D images. The convolution is then given by a **kernel** matrix $K$:

$$M'_{i, j} = \sum_{k, l}  \, M_{i- k, j - l} \, K_{k, l},$$

where the sum is over the possible values of $k$ and $l$ in the window. Again we think of the window as being *centered* at $(i, j)$.

A common notation for this operation is $*$:

$$M' = M * K.$$
"""

# ╔═╡ 7c1bc062-ee15-11ea-30b1-1b1e76520f13
md"""
#### Exercise 4.1
👉 Write a function `extend_mat` that takes a matrix `M` and indices `i` and `j`, and returns the closest element of the matrix.
"""

# ╔═╡ 7c2ec6c6-ee15-11ea-2d7d-0d9401a5e5d1
function extend_mat(M::AbstractMatrix, i, j)
	(nr, nc) = size(M)
	t_i, t_j = 0,0
	if i < 1
		t_i = 1
	elseif i > nr
		t_i = nr
	else
		t_i = i
	end

	if j < 1
		t_j = 1
	elseif j > nc
		t_j = nc
	else
		t_j = j
	end
	return M[t_i,t_j]
end

# ╔═╡ 9afc4dca-ee16-11ea-354f-1d827aaa61d2
md"_Let's test it!_"

# ╔═╡ cf6b05e2-ee16-11ea-3317-8919565cb56e
small_image = Gray.(rand(5,5))

# ╔═╡ e3616062-ee27-11ea-04a9-b9ec60842a64
md"Extended with `0`:"

# ╔═╡ e5b6cd34-ee27-11ea-0d60-bd4796540b18
[get(small_image, (i, j), Gray(0)) for (i,j) in Iterators.product(-1:7,-1:7)]

# ╔═╡ d06ea762-ee27-11ea-2e9c-1bcff86a3fe0
md"Extended with your `extend`:"

# ╔═╡ e1dc0622-ee16-11ea-274a-3b6ec9e15ab5
[extend_mat(small_image, i, j) for (i,j) in Iterators.product(-1:7,-1:7)]

# ╔═╡ 3cd535e4-ee26-11ea-2482-fb4ad43dda19
let
	philip_head = philip[250:430,110:230]
	[extend_mat(philip_head, i, j) for (i,j) in Iterators.product(-50:size(philip_head,1)+51, (-50:size(philip_head,2)+51))]
end

# ╔═╡ 7c41f0ca-ee15-11ea-05fb-d97a836659af
md"""
#### Exercise 4.2
👉 Implement a function `convolve_image(M, K)`.
"""

# ╔═╡ 8b96e0bc-ee15-11ea-11cd-cfecea7075a0
function convolve_image(M::AbstractMatrix, K::AbstractMatrix)
	nc_K, nr_K = size(K)
	nc_M, nr_M = size(M)
	M_pad = [extend_mat(M,i,j) for (i,j) in Iterators.product(1-nc_K÷2:nc_M+nc_K÷2,1-nr_K÷2:nr_M+nr_K÷2)]
	return [sum(M_pad[i:i+nc_K-1,j:j+nr_K-1].*K) for (i,j) in Iterators.product(1:nc_M, 1:nr_M)]
end

# ╔═╡ 5a5135c6-ee1e-11ea-05dc-eb0c683c2ce5
md"_Let's test it out! 🎃_"

# ╔═╡ 577c6daa-ee1e-11ea-1275-b7abc7a27d73
test_image_with_border = [get(small_image, (i, j), Gray(0)) for (i,j) in Iterators.product(-1:7,-1:7)]

# ╔═╡ 09e8fde0-48ef-11eb-2dfa-11cad1f2ea84
M = test_image_with_border

# ╔═╡ b3b85330-48ee-11eb-349a-cd5bc8f200ad
begin
	K = K_test
	nc_K, nr_K = size(K)
	nc_M, nr_M = size(M)
	M_pad = [extend_mat(M,i,j) for (i,j) in Iterators.product(1-nc_K÷2:nc_M+nc_K÷2,1-nr_K÷2:nr_M+nr_K÷2)]
	[sum(M_pad[i:i+nc_K -1,j:j+nr_K -1]*K) for (i,j) in Iterators.product(1:nc_M, 1:nr_M)]
end

# ╔═╡ 275a99c8-ee1e-11ea-0a76-93e3618c9588
K_test = [
	0  	0 	0
	0 	1 	0
	0   0 	0
]

# ╔═╡ 42dfa206-ee1e-11ea-1fcd-21671042064c
convolve_image(test_image_with_border, K_test)

# ╔═╡ 6e53c2e6-ee1e-11ea-21bd-c9c05381be07
md"_Edit_ `K_test` _to create your own test case!_"

# ╔═╡ e7f8b41a-ee25-11ea-287a-e75d33fbd98b
convolve_image(philip, K_test)

# ╔═╡ 7c50ea80-ee15-11ea-328f-6b4e4ff20b7e
md"""
#### Exercise 4.3
👉 Apply a **Gaussian blur** to an image.

Here, the 2D Gaussian kernel will be defined as

$$G(x,y)=\frac{1}{2\pi \sigma^2}e^{\frac{-(x^2+y^2)}{2\sigma^2}}$$
"""

# ╔═╡ 8a335044-ee19-11ea-0255-b9391246d231
md"""
---

You can create all sorts of effects by choosing the kernel in a smart way. Today, we will implement two special kernels, to produce a **Gaussian blur** and a **Sobel edge detect** filter.

Make sure that you have watched [the lecture](https://www.youtube.com/watch?v=8rrHTtUzyZA) about convolutions!
"""

# ╔═╡ e918f550-48f0-11eb-0240-9bdc1b3ad4a6
function normalize(K)
	return K./sum(K)
end

# ╔═╡ aad67fd0-ee15-11ea-00d4-274ec3cda3a3
function with_gaussian_blur(image, n=7)
	g(x,y) = exp(-(x^2+y^2)/2)/(2π)

	K = normalize([g(i-(n+1)/2,j-(n+1)/2) for (i,j) in Iterators.product(1:n,1:n)])
	# @show K
	return convolve_image(image,K)
end

# ╔═╡ 8ae59674-ee18-11ea-3815-f50713d0fa08
md"_Let's make it interactive. 💫_"

# ╔═╡ d4a85820-48f2-11eb-1040-99b046126b5e
K, img=with_gaussian_blur(philip,7)

# ╔═╡ b41383d0-48f4-11eb-2c7f-cfa803500f77
K

# ╔═╡ 7c6642a6-ee15-11ea-0526-a1aac4286cdd
md"""
#### Exercise 4.4
👉 Create a **Sobel edge detection filter**.

Here, we will need to create two separate filters that separately detect edges in the horizontal and vertical directions:

```math
\begin{align}

G_x &= \left(\begin{bmatrix}
1 \\
2 \\
1 \\
\end{bmatrix} \otimes [1~0~-1]
\right) * A = \begin{bmatrix}
1 & 0 & -1 \\
2 & 0 & -2 \\
1 & 0 & -1 \\
\end{bmatrix}*A\\
G_y &= \left(
\begin{bmatrix}
1 \\
0 \\
-1 \\
\end{bmatrix} \otimes [1~2~1]
\right) * A = \begin{bmatrix}
1 & 2 & 1 \\
0 & 0 & 0 \\
-1 & -2 & -1 \\
\end{bmatrix}*A
\end{align}
```
Here $A$ is the array corresponding to your image.
We can think of these as derivatives in the $x$ and $y$ directions.

Then we combine them by finding the magnitude of the **gradient** (in the sense of multivariate calculus) by defining

$$G_\text{total} = \sqrt{G_x^2 + G_y^2}.$$

For simplicity you can choose one of the "channels" (colours) in the image to apply this to.
"""

# ╔═╡ 9eeb876c-ee15-11ea-1794-d3ea79f47b75
function with_sobel_edge_detect(image)
	a = [1 2 1]
	b = [1 0 -1]
	K_1 = a'*b
	K_2 = b'*a
	G_x = convolve_image(image,K_1)
	G_y = convolve_image(image,K_2)
	G_t = (G_x.*G_x + G_y.*G_y).^.5
	@show size(K_1)
	return G_t
end

# ╔═╡ 1b85ee76-ee10-11ea-36d7-978340ef61e6
md"""
## **Exercise 5** - _Lecture transcript_
_(MIT students only)_

Please see the Canvas post for transcript document for week 1 [here](https://canvas.mit.edu/courses/5637/discussion_topics/27880).

We need each of you to correct about 100 lines (see instructions in the beginning of the document.)

👉 Please mention the name of the video and the line ranges you edited:
"""

# ╔═╡ 477d0a3c-ee10-11ea-11cf-07b0e0ce6818
lines_i_edited = md"""
Convolution, lines 100-0 (_for example_)
"""

# ╔═╡ 8ffe16ce-ee20-11ea-18bd-15640f94b839
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

# ╔═╡ 5516c800-edee-11ea-12cf-3f8c082ef0ef
hint(text) = Markdown.MD(Markdown.Admonition("hint", "Hint", [text]))

# ╔═╡ b1d5ca28-edf6-11ea-269e-75a9fb549f1d
hint(md"You can find out more about any function (like `rand`) by creating a new cell and typing:

```
?rand
```

Once the Live Docs are open, you can select any code to learn more about it. It might be useful to leave it open all the time, and get documentation while you type code.")

# ╔═╡ f6ef2c2e-ee07-11ea-13a8-2512e7d94426
hint(md"The `rand` function generates (uniform) random floating-point numbers between $0$ and $1$.")

# ╔═╡ ea435e58-ee11-11ea-3785-01af8dd72360
hint(md"Have a look at Exercise 2 to see an example of adding interactivity with a slider. You can read the [Interactivity](./sample/Interactivity.jl) and the [PlutoUI](./sample/PlutoUI.jl) sample notebooks _(right click -> Open in new tab)_ to learn more.")

# ╔═╡ e9aadeee-ee1d-11ea-3525-95f6ba5fda31
hint(md"`l = (length(k) - 1) ÷ 2`")

# ╔═╡ 649df270-ee24-11ea-397e-79c4355e38db
hint(md"`num_rows, num_columns = size(M)`")

# ╔═╡ 0cabed84-ee1e-11ea-11c1-7d8a4b4ad1af
hint(md"`num_rows, num_columns = size(K)`")

# ╔═╡ 57360a7a-edee-11ea-0c28-91463ece500d
almost(text) = Markdown.MD(Markdown.Admonition("warning", "Almost there!", [text]))

# ╔═╡ dcb8324c-edee-11ea-17ff-375ff5078f43
still_missing(text=md"Replace `missing` with your answer.") = Markdown.MD(Markdown.Admonition("warning", "Here we go!", [text]))

# ╔═╡ 58af703c-edee-11ea-2963-f52e78fc2412
keep_working(text=md"The answer is not quite right.") = Markdown.MD(Markdown.Admonition("danger", "Keep working on it!", [text]))

# ╔═╡ f3d00a9a-edf3-11ea-07b3-1db5c6d0b3cf
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."]

# ╔═╡ 5aa9dfb2-edee-11ea-3754-c368fb40637c
correct(text=rand(yays)) = Markdown.MD(Markdown.Admonition("correct", "Got it!", [text]))

# ╔═╡ 74d44e22-edee-11ea-09a0-69aa0aba3281
not_defined(variable_name) = Markdown.MD(Markdown.Admonition("danger", "Oopsie!", [md"Make sure that you define a variable called **$(Markdown.Code(string(variable_name)))**"]))

# ╔═╡ 397941fc-edee-11ea-33f2-5d46c759fbf7
if !@isdefined(random_vect)
	not_defined(:random_vect)
elseif ismissing(random_vect)
	still_missing()
elseif !(random_vect isa Vector)
	keep_working(md"`random_vect` should be a `Vector`.")
elseif length(random_vect) != 10
	keep_working(md"`random_vect` does not have the correct size.")
else
	correct()
end

# ╔═╡ 38dc80a0-edef-11ea-10e9-615255a4588c
if !@isdefined(mean)
	not_defined(:mean)
else
	let
		result = mean([1,2,3])
		if ismissing(result)
			still_missing()
		elseif isnothing(result)
			keep_working(md"Did you forget to write `return`?")
		elseif result != 2
			keep_working()
		else
			correct()
		end
	end
end

# ╔═╡ 2b1ccaca-edee-11ea-34b0-c51659f844d0
if !@isdefined(m)
	not_defined(:m)
elseif ismissing(m)
	still_missing()
elseif !(m isa Number)
	keep_working(md"`m` should be a number.")
elseif m != mean(random_vect)
	keep_working()
else
	correct()
end

# ╔═╡ e3394c8a-edf0-11ea-1bb8-619f7abb6881
if !@isdefined(create_bar)
	not_defined(:create_bar)
else
	let
		result = create_bar()
		if ismissing(result)
			still_missing()
		elseif isnothing(result)
			keep_working(md"Did you forget to write `return`?")
		elseif !(result isa Vector) || length(result) != 100
			keep_working(md"The result should be a `Vector` with 100 elements.")
		elseif result[[1,50,100]] != [0,1,0]
			keep_working()
		else
			correct()
		end
	end
end

# ╔═╡ adfbe9b2-ed6c-11ea-09ac-675262f420df
if !@isdefined(vecvec_to_matrix)
	not_defined(:vecvec_to_matrix)
else
	let
		input = [[6,7],[8,9]]

		result = vecvec_to_matrix(input)
		shouldbe = [6 7; 8 9]

		if ismissing(result)
			still_missing()
		elseif isnothing(result)
			keep_working(md"Did you forget to write `return`?")
		elseif !(result isa Matrix)
			keep_working(md"The result should be a `Matrix`")
		elseif result != shouldbe && result != shouldbe'
			keep_working()
		else
			correct()
		end
	end
end

# ╔═╡ e06b7fbc-edf2-11ea-1708-fb32599dded3
if !@isdefined(matrix_to_vecvec)
	not_defined(:matrix_to_vecvec)
else
	let
		input = [6 7 8; 8 9 10]
		result = matrix_to_vecvec(input)
		shouldbe = [[6,7,8],[8,9,10]]
		shouldbe2 = [[6,8], [7,9], [8,10]]

		if ismissing(result)
			still_missing()
		elseif isnothing(result)
			keep_working(md"Did you forget to write `return`?")
		elseif result != shouldbe && result != shouldbe2
			keep_working()
		else
			correct()
		end
	end
end

# ╔═╡ 4d0158d0-ee0d-11ea-17c3-c169d4284acb
if !@isdefined(mean_colors)
	not_defined(:mean_colors)
else
	let
		input = reshape([RGB(1.0, 1.0, 1.0), RGB(1.0, 1.0, 0.0)], (2,1))

		result = mean_colors(input)
		shouldbe = (1.0, 1.0, 0.5)
		shouldbe2 = RGB(shouldbe...)

		if ismissing(result)
			still_missing()
		elseif isnothing(result)
			keep_working(md"Did you forget to write `return`?")
		elseif !(result == shouldbe) && !(result == shouldbe2)
			keep_working()
		else
			correct()
		end
	end
end

# ╔═╡ c905b73e-ee1a-11ea-2e36-23b8e73bfdb6
if !@isdefined(quantize)
	not_defined(:quantize)
else
	let
		result = quantize(.3)

		if ismissing(result)
			still_missing()
		elseif isnothing(result)
			keep_working(md"Did you forget to write `return`?")
		elseif result != .3
			if quantize(0.35) == .3
				almost(md"What should quantize(`0.2`) be?")
			else
				keep_working()
			end
		else
			correct()
		end
	end
end

# ╔═╡ bcf98dfc-ee1b-11ea-21d0-c14439500971
if !@isdefined(extend)
	not_defined(:extend)
else
	let
		result = extend([6,7],-10)

		if ismissing(result)
			still_missing()
		elseif isnothing(result)
			keep_working(md"Did you forget to write `return`?")
		elseif result != 6 || extend([6,7],10) != 7
			keep_working()
		else
			correct()
		end
	end
end

# ╔═╡ 7ffd14f8-ee1d-11ea-0343-b54fb0333aea
if !@isdefined(convolve_vector)
	not_defined(:convolve_vector)
else
	let
		x = [1, 10, 100]
		result = convolve_vector(x, [0, 1, 1])
		shouldbe = [11, 110, 200]
		shouldbe2 = [2, 11, 110]

		if ismissing(result)
			still_missing()
		elseif isnothing(result)
			keep_working(md"Did you forget to write `return`?")
		elseif !(result isa AbstractVector)
			keep_working(md"The returned object is not a `Vector`.")
		elseif size(result) != size(x)
			keep_working(md"The returned vector has the wrong dimensions.")
		elseif result != shouldbe && result != shouldbe2
			keep_working()
		else
			correct()
		end
	end
end

# ╔═╡ efd1ceb4-ee1c-11ea-350e-f7e3ea059024
if !@isdefined(extend_mat)
	not_defined(:extend_mat)
else
	let
		input = [42 37; 1 0]
		result = extend_mat(input, -2, -2)

		if ismissing(result)
			still_missing()
		elseif isnothing(result)
			keep_working(md"Did you forget to write `return`?")
		elseif result != 42 || extend_mat(input, -1, 3) != 37
			keep_working()
		else
			correct()
		end
	end
end

# ╔═╡ 115ded8c-ee0a-11ea-3493-89487315feb7
bigbreak = html"<br><br><br><br><br>";

# ╔═╡ 54056a02-ee0a-11ea-101f-47feb6623bec
bigbreak

# ╔═╡ 45815734-ee0a-11ea-2982-595e1fc0e7b1
bigbreak

# ╔═╡ 4139ee66-ee0a-11ea-2282-15d63bcca8b8
bigbreak

# ╔═╡ 27847dc4-ee0a-11ea-0651-ebbbb3cfd58c
bigbreak

# ╔═╡ 0001f782-ee0e-11ea-1fb4-2b5ef3d241e2
bigbreak

# ╔═╡ 91f4778e-ee20-11ea-1b7e-2b0892bd3c0f
bigbreak

# ╔═╡ 5842895a-ee10-11ea-119d-81e4c4c8c53b
bigbreak

# ╔═╡ dfb7c6be-ee0d-11ea-194e-9758857f7b20
function camera_input(;max_size=200, default_url="https://i.imgur.com/SUmi94P.png")
"""
<span class="pl-image waiting-for-permission">
<style>

	.pl-image.popped-out {
		position: fixed;
		top: 0;
		right: 0;
		z-index: 5;
	}

	.pl-image #video-container {
		width: 250px;
	}

	.pl-image video {
		border-radius: 1rem 1rem 0 0;
	}
	.pl-image.waiting-for-permission #video-container {
		display: none;
	}
	.pl-image #prompt {
		display: none;
	}
	.pl-image.waiting-for-permission #prompt {
		width: 250px;
		height: 200px;
		display: grid;
		place-items: center;
		font-family: monospace;
		font-weight: bold;
		text-decoration: underline;
		cursor: pointer;
		border: 5px dashed rgba(0,0,0,.5);
	}

	.pl-image video {
		display: block;
	}
	.pl-image .bar {
		width: inherit;
		display: flex;
		z-index: 6;
	}
	.pl-image .bar#top {
		position: absolute;
		flex-direction: column;
	}

	.pl-image .bar#bottom {
		background: black;
		border-radius: 0 0 1rem 1rem;
	}
	.pl-image .bar button {
		flex: 0 0 auto;
		background: rgba(255,255,255,.8);
		border: none;
		width: 2rem;
		height: 2rem;
		border-radius: 100%;
		cursor: pointer;
		z-index: 7;
	}
	.pl-image .bar button#shutter {
		width: 3rem;
		height: 3rem;
		margin: -1.5rem auto .2rem auto;
	}

	.pl-image video.takepicture {
		animation: pictureflash 200ms linear;
	}

	@keyframes pictureflash {
		0% {
			filter: grayscale(1.0) contrast(2.0);
		}

		100% {
			filter: grayscale(0.0) contrast(1.0);
		}
	}
</style>

	<div id="video-container">
		<div id="top" class="bar">
			<button id="stop" title="Stop video">✖</button>
			<button id="pop-out" title="Pop out/pop in">⏏</button>
		</div>
		<video playsinline autoplay></video>
		<div id="bottom" class="bar">
		<button id="shutter" title="Click to take a picture">📷</button>
		</div>
	</div>

	<div id="prompt">
		<span>
		Enable webcam
		</span>
	</div>

<script>
	// based on https://github.com/fonsp/printi-static (by the same author)

	const span = currentScript.parentElement
	const video = span.querySelector("video")
	const popout = span.querySelector("button#pop-out")
	const stop = span.querySelector("button#stop")
	const shutter = span.querySelector("button#shutter")
	const prompt = span.querySelector(".pl-image #prompt")

	const maxsize = $(max_size)

	const send_source = (source, src_width, src_height) => {
		const scale = Math.min(1.0, maxsize / src_width, maxsize / src_height)

		const width = Math.floor(src_width * scale)
		const height = Math.floor(src_height * scale)

		const canvas = html`<canvas width=\${width} height=\${height}>`
		const ctx = canvas.getContext("2d")
		ctx.drawImage(source, 0, 0, width, height)

		span.value = {
			width: width,
			height: height,
			data: ctx.getImageData(0, 0, width, height).data,
		}
		span.dispatchEvent(new CustomEvent("input"))
	}

	const clear_camera = () => {
		window.stream.getTracks().forEach(s => s.stop());
		video.srcObject = null;

		span.classList.add("waiting-for-permission");
	}

	prompt.onclick = () => {
		navigator.mediaDevices.getUserMedia({
			audio: false,
			video: {
				facingMode: "environment",
			},
		}).then(function(stream) {

			stream.onend = console.log

			window.stream = stream
			video.srcObject = stream
			window.cameraConnected = true
			video.controls = false
			video.play()
			video.controls = false

			span.classList.remove("waiting-for-permission");

		}).catch(function(error) {
			console.log(error)
		});
	}
	stop.onclick = () => {
		clear_camera()
	}
	popout.onclick = () => {
		span.classList.toggle("popped-out")
	}

	shutter.onclick = () => {
		const cl = video.classList
		cl.remove("takepicture")
		void video.offsetHeight
		cl.add("takepicture")
		video.play()
		video.controls = false
		console.log(video)
		send_source(video, video.videoWidth, video.videoHeight)
	}

	document.addEventListener("visibilitychange", () => {
		if (document.visibilityState != "visible") {
			clear_camera()
		}
	})

	// Set a default image

	const img = html`<img crossOrigin="anonymous">`

	img.onload = () => {
	console.log("helloo")
		send_source(img, img.width, img.height)
	}
	img.src = "$(default_url)"
	console.log(img)
</script>
</span>
""" |> HTML
end

# ╔═╡ 94c0798e-ee18-11ea-3212-1533753eabb6
@bind gauss_raw_camera_data camera_input(;max_size=100)

# ╔═╡ 1a0324de-ee19-11ea-1d4d-db37f4136ad3
@bind sobel_raw_camera_data camera_input(;max_size=100)

# ╔═╡ e15ad330-ee0d-11ea-25b6-1b1b3f3d7888

function process_raw_camera_data(raw_camera_data)
	# the raw image data is a long byte array, we need to transform it into something
	# more "Julian" - something with more _structure_.

	# The encoding of the raw byte stream is:
	# every 4 bytes is a single pixel
	# every pixel has 4 values: Red, Green, Blue, Alpha
	# (we ignore alpha for this notebook)

	# So to get the red values for each pixel, we take every 4th value, starting at
	# the 1st:
	reds_flat = UInt8.(raw_camera_data["data"][1:4:end])
	greens_flat = UInt8.(raw_camera_data["data"][2:4:end])
	blues_flat = UInt8.(raw_camera_data["data"][3:4:end])

	# but these are still 1-dimensional arrays, nicknamed 'flat' arrays
	# We will 'reshape' this into 2D arrays:

	width = raw_camera_data["width"]
	height = raw_camera_data["height"]

	# shuffle and flip to get it in the right shape
	reds = reshape(reds_flat, (width, height))' / 255.0
	greens = reshape(greens_flat, (width, height))' / 255.0
	blues = reshape(blues_flat, (width, height))' / 255.0

	# we have our 2D array for each color
	# Let's create a single 2D array, where each value contains the R, G and B value of
	# that pixel

	RGB.(reds, greens, blues)
end

# ╔═╡ f461f5f2-ee18-11ea-3d03-95f57f9bf09e
gauss_camera_image = process_raw_camera_data(gauss_raw_camera_data);

# ╔═╡ a75701c4-ee18-11ea-2863-d3042e71a68b
with_gaussian_blur(gauss_camera_image)

# ╔═╡ 1ff6b5cc-ee19-11ea-2ca8-7f00c204f587
sobel_camera_image = Gray.(process_raw_camera_data(sobel_raw_camera_data));

# ╔═╡ 1bf94c00-ee19-11ea-0e3c-e12bc68d8e28
with_sobel_edge_detect(sobel_camera_image)

# ╔═╡ Cell order:
# ╠═83eb9ca0-ed68-11ea-0bc5-99a09c68f867
# ╟─8ef13896-ed68-11ea-160b-3550eeabbd7d
# ╟─ac8ff080-ed61-11ea-3650-d9df06123e1f
# ╠═911ccbce-ed68-11ea-3606-0384e7580d7c
# ╟─5f95e01a-ee0a-11ea-030c-9dba276aba92
# ╠═65780f00-ed6b-11ea-1ecf-8b35523a7ac0
# ╟─67461396-ee0a-11ea-3679-f31d46baa9b4
# ╠═74b008f6-ed6b-11ea-291f-b3791d6d1b35
# ╟─54056a02-ee0a-11ea-101f-47feb6623bec
# ╟─540ccfcc-ee0a-11ea-15dc-4f8120063397
# ╟─467856dc-eded-11ea-0f83-13d939021ef3
# ╠═56ced344-eded-11ea-3e81-3936e9ad5777
# ╟─ad6a33b0-eded-11ea-324c-cfabfd658b56
# ╠═f51333a6-eded-11ea-34e6-bfbb3a69bcb0
# ╟─b18e2c54-edf1-11ea-0cbf-85946d64b6a2
# ╟─397941fc-edee-11ea-33f2-5d46c759fbf7
# ╟─b1d5ca28-edf6-11ea-269e-75a9fb549f1d
# ╟─cf738088-eded-11ea-2915-61735c2aa990
# ╠═0ffa8354-edee-11ea-2883-9d5bfea4a236
# ╠═1f104ce4-ee0e-11ea-2029-1d9c817175af
# ╟─38dc80a0-edef-11ea-10e9-615255a4588c
# ╟─1f229ca4-edee-11ea-2c56-bb00cc6ea53c
# ╠═2a391708-edee-11ea-124e-d14698171b68
# ╟─2b1ccaca-edee-11ea-34b0-c51659f844d0
# ╟─e2863d4c-edef-11ea-1d67-332ddca03cc4
# ╠═ec5efe8c-edef-11ea-2c6f-afaaeb5bc50c
# ╟─29e10640-edf0-11ea-0398-17dbf4242de3
# ╟─6f67657e-ee1a-11ea-0c2f-3d567bcfa6ea
# ╠═38155b5a-edf0-11ea-3e3f-7163da7433fb
# ╠═73ef1d50-edf0-11ea-343c-d71706874c82
# ╟─a5f8bafe-edf0-11ea-0da3-3330861ae43a
# ╠═b6b65b94-edf0-11ea-3686-fbff0ff53d08
# ╟─d862fb16-edf1-11ea-36ec-615d521e6bc0
# ╠═e3394c8a-edf0-11ea-1bb8-619f7abb6881
# ╟─22f28dae-edf2-11ea-25b5-11c369ae1253
# ╠═8c19fb72-ed6c-11ea-2728-3fa9219eddc4
# ╠═c4761a7e-edf2-11ea-1e75-118e73dadbed
# ╠═40881910-4865-11eb-1ff7-3798bccfbc7b
# ╟─adfbe9b2-ed6c-11ea-09ac-675262f420df
# ╟─393667ca-edf2-11ea-09c5-c5d292d5e896
# ╠═9f1c6d04-ed6c-11ea-007b-75e7e780703d
# ╠═70955aca-ed6e-11ea-2330-89b4d20b1795
# ╠═ab4fc2c0-4865-11eb-0a6c-db550806559a
# ╠═e06b7fbc-edf2-11ea-1708-fb32599dded3
# ╟─5da8cbe8-eded-11ea-2e43-c5b7cc71e133
# ╟─45815734-ee0a-11ea-2982-595e1fc0e7b1
# ╟─e083b3e8-ed61-11ea-2ec9-217820b0a1b4
# ╠═c5484572-ee05-11ea-0424-f37295c3072d
# ╠═c8ecfe5c-ee05-11ea-322b-4b2714898831
# ╟─e86ed944-ee05-11ea-3e0f-d70fc73b789c
# ╟─c54ccdea-ee05-11ea-0365-23aaf053b7d7
# ╠═f6898df6-ee07-11ea-2838-fde9bc739c11
# ╠═5be9b144-ee0d-11ea-2a8d-8775de265a1d
# ╟─4d0158d0-ee0d-11ea-17c3-c169d4284acb
# ╟─f68d4a36-ee07-11ea-0832-0360530f102e
# ╠═f6991a50-ee07-11ea-0bc4-1d68eb028e6a
# ╠═f6a655f8-ee07-11ea-13b6-43ca404ddfc7
# ╠═1fc7dd0e-4869-11eb-1761-552ef1069515
# ╠═ce65f880-4868-11eb-288c-ad505b2bacf8
# ╟─c905b73e-ee1a-11ea-2e36-23b8e73bfdb6
# ╟─f6b218c0-ee07-11ea-2adb-1968c4fd473a
# ╟─f6bf64da-ee07-11ea-3efb-05af01b14f67
# ╟─25dad7ce-ee0b-11ea-3e20-5f3019dd7fa3
# ╠═9751586e-ee0c-11ea-0cbb-b7eda92977c9
# ╟─f6cc03a0-ee07-11ea-17d8-013991514d42
# ╠═63e8d636-ee0b-11ea-173d-bd3327347d55
# ╟─2cc2f84e-ee0d-11ea-373b-e7ad3204bb00
# ╟─b8f26960-ee0a-11ea-05b9-3f4bc1099050
# ╠═5de3a22e-ee0b-11ea-230f-35df4ca3c96d
# ╠═4e21e0c4-ee0b-11ea-3d65-b311ae3f98e9
# ╠═6dbf67ce-ee0b-11ea-3b71-abc05a64dc43
# ╟─846b1330-ee0b-11ea-3579-7d90fafd7290
# ╠═943103e2-ee0b-11ea-33aa-75a8a1529931
# ╟─f6d6c71a-ee07-11ea-2b63-d759af80707b
# ╠═70ac1b40-486b-11eb-2b38-1920b4dcf32b
# ╠═f6e2cb2a-ee07-11ea-06ee-1b77e34c1e91
# ╟─f6ef2c2e-ee07-11ea-13a8-2512e7d94426
# ╟─f6fc1312-ee07-11ea-39a0-299b67aee3d8
# ╟─774b4ce6-ee1b-11ea-2b48-e38ee25fc89b
# ╠═7e4aeb70-ee1b-11ea-100f-1952ba66f80f
# ╟─6a05f568-ee1b-11ea-3b6c-83b6ada3680f
# ╟─f70823d2-ee07-11ea-2bb3-01425212aaf9
# ╠═e70a84d4-ee0c-11ea-0640-bf78653ba102
# ╠═ac15e0d0-ee0c-11ea-1eaf-d7f88b5df1d7
# ╟─9604bc44-ee1b-11ea-28f8-7f7af8d0cbb2
# ╟─f714699e-ee07-11ea-08b6-5f5169861b57
# ╠═bdc2df7c-ee0c-11ea-2e9f-7d2c085617c1
# ╟─81510a30-ee0e-11ea-0062-8b3327428f9d
# ╠═6b30dc38-ed6b-11ea-10f3-ab3f121bf4b8
# ╟─e3b03628-ee05-11ea-23b6-27c7b0210532
# ╟─4139ee66-ee0a-11ea-2282-15d63bcca8b8
# ╟─e08781fa-ed61-11ea-13ae-91a49b5eb74a
# ╟─7fc8ee1c-ee09-11ea-1382-ad21d5373308
# ╠═7fcd6230-ee09-11ea-314f-a542d00d582e
# ╠═7fdb34dc-ee09-11ea-366b-ffe10d1aa845
# ╟─7fe9153e-ee09-11ea-15b3-6f24fcc20734
# ╟─80108d80-ee09-11ea-0368-31546eb0d3cc
# ╠═01070e28-ee0f-11ea-1928-a7919d452bdd
# ╟─7522f81e-ee1c-11ea-35af-a17eb257ff1a
# ╟─801d90c0-ee09-11ea-28d6-61b806de26dc
# ╠═802bec56-ee09-11ea-043e-51cf1db02a34
# ╟─b7f3994c-ee1b-11ea-211a-d144db8eafc2
# ╠═803905b2-ee09-11ea-2d52-e77ff79693b0
# ╠═80479d98-ee09-11ea-169e-d166eef65874
# ╠═805691ce-ee09-11ea-053d-6d2e299ee123
# ╟─806e5766-ee0f-11ea-1efc-d753cd83d086
# ╟─38da843a-ee0f-11ea-01df-bfa8b1317d36
# ╟─9bde9f92-ee0f-11ea-27f8-ffef5fce2b3c
# ╠═45c4da9a-ee0f-11ea-2c5b-1f6704559137
# ╟─bcf98dfc-ee1b-11ea-21d0-c14439500971
# ╟─80664e8c-ee09-11ea-0702-711bce271315
# ╠═807e5662-ee09-11ea-3005-21fdcc36b023
# ╟─808deca8-ee09-11ea-0ee3-1586fa1ce282
# ╟─809f5330-ee09-11ea-0e5b-415044b6ac1f
# ╠═4c5cdff0-48ce-11eb-10b8-691432d0012a
# ╠═ca1ac5f4-ee1c-11ea-3d00-ff5268866f87
# ╠═5cde0890-48ce-11eb-03a2-fda7bcffcfec
# ╟─ea435e58-ee11-11ea-3785-01af8dd72360
# ╟─80ab64f4-ee09-11ea-29b4-498112ed0799
# ╠═28e20950-ee0c-11ea-0e0a-b5f2e570b56e
# ╟─e9aadeee-ee1d-11ea-3525-95f6ba5fda31
# ╟─5eea882c-ee13-11ea-0d56-af81ecd30a4a
# ╠═93284f92-ee12-11ea-0342-833b1a30625c
# ╟─cf73f9f8-ee12-11ea-39ae-0107e9107ef5
# ╟─7ffd14f8-ee1d-11ea-0343-b54fb0333aea
# ╟─80b7566a-ee09-11ea-3939-6fab470f9ec8
# ╠═1c8b4658-ee0c-11ea-2ede-9b9ed7d3125e
# ╠═664b6270-48e1-11eb-1c33-3b4d2fa7d90d
# ╠═bae5f940-48e0-11eb-2762-a38b8a33bf0d
# ╟─f8bd22b8-ee14-11ea-04aa-ab16fd01826e
# ╠═576f53c0-48e0-11eb-0d9e-d9e5f1e794eb
# ╠═2a9dd06a-ee13-11ea-3f84-67bb309c77a8
# ╠═b424e2aa-ee14-11ea-33fa-35491e0b9c9d
# ╠═38eb92f6-ee13-11ea-14d7-a503ac04302e
# ╠═bc1c20a4-ee14-11ea-3525-63c9fa78f089
# ╠═24c21c7c-ee14-11ea-1512-677980db1288
# ╟─27847dc4-ee0a-11ea-0651-ebbbb3cfd58c
# ╟─b01858b6-edf3-11ea-0826-938d33c19a43
# ╟─7c1bc062-ee15-11ea-30b1-1b1e76520f13
# ╠═7c2ec6c6-ee15-11ea-2d7d-0d9401a5e5d1
# ╟─649df270-ee24-11ea-397e-79c4355e38db
# ╟─9afc4dca-ee16-11ea-354f-1d827aaa61d2
# ╠═cf6b05e2-ee16-11ea-3317-8919565cb56e
# ╟─e3616062-ee27-11ea-04a9-b9ec60842a64
# ╟─e5b6cd34-ee27-11ea-0d60-bd4796540b18
# ╟─d06ea762-ee27-11ea-2e9c-1bcff86a3fe0
# ╠═e1dc0622-ee16-11ea-274a-3b6ec9e15ab5
# ╟─efd1ceb4-ee1c-11ea-350e-f7e3ea059024
# ╠═3cd535e4-ee26-11ea-2482-fb4ad43dda19
# ╟─7c41f0ca-ee15-11ea-05fb-d97a836659af
# ╠═8b96e0bc-ee15-11ea-11cd-cfecea7075a0
# ╟─0cabed84-ee1e-11ea-11c1-7d8a4b4ad1af
# ╟─5a5135c6-ee1e-11ea-05dc-eb0c683c2ce5
# ╟─577c6daa-ee1e-11ea-1275-b7abc7a27d73
# ╠═09e8fde0-48ef-11eb-2dfa-11cad1f2ea84
# ╠═b3b85330-48ee-11eb-349a-cd5bc8f200ad
# ╠═275a99c8-ee1e-11ea-0a76-93e3618c9588
# ╠═42dfa206-ee1e-11ea-1fcd-21671042064c
# ╟─6e53c2e6-ee1e-11ea-21bd-c9c05381be07
# ╠═e7f8b41a-ee25-11ea-287a-e75d33fbd98b
# ╟─7c50ea80-ee15-11ea-328f-6b4e4ff20b7e
# ╟─8a335044-ee19-11ea-0255-b9391246d231
# ╠═e918f550-48f0-11eb-0240-9bdc1b3ad4a6
# ╠═aad67fd0-ee15-11ea-00d4-274ec3cda3a3
# ╟─8ae59674-ee18-11ea-3815-f50713d0fa08
# ╠═94c0798e-ee18-11ea-3212-1533753eabb6
# ╟─a75701c4-ee18-11ea-2863-d3042e71a68b
# ╠═d4a85820-48f2-11eb-1040-99b046126b5e
# ╠═b41383d0-48f4-11eb-2c7f-cfa803500f77
# ╟─f461f5f2-ee18-11ea-3d03-95f57f9bf09e
# ╟─7c6642a6-ee15-11ea-0526-a1aac4286cdd
# ╠═9eeb876c-ee15-11ea-1794-d3ea79f47b75
# ╟─1a0324de-ee19-11ea-1d4d-db37f4136ad3
# ╠═1bf94c00-ee19-11ea-0e3c-e12bc68d8e28
# ╟─1ff6b5cc-ee19-11ea-2ca8-7f00c204f587
# ╟─0001f782-ee0e-11ea-1fb4-2b5ef3d241e2
# ╠═1b85ee76-ee10-11ea-36d7-978340ef61e6
# ╠═477d0a3c-ee10-11ea-11cf-07b0e0ce6818
# ╟─91f4778e-ee20-11ea-1b7e-2b0892bd3c0f
# ╟─8ffe16ce-ee20-11ea-18bd-15640f94b839
# ╟─5842895a-ee10-11ea-119d-81e4c4c8c53b
# ╟─5516c800-edee-11ea-12cf-3f8c082ef0ef
# ╟─57360a7a-edee-11ea-0c28-91463ece500d
# ╟─dcb8324c-edee-11ea-17ff-375ff5078f43
# ╟─58af703c-edee-11ea-2963-f52e78fc2412
# ╟─f3d00a9a-edf3-11ea-07b3-1db5c6d0b3cf
# ╟─5aa9dfb2-edee-11ea-3754-c368fb40637c
# ╟─74d44e22-edee-11ea-09a0-69aa0aba3281
# ╟─115ded8c-ee0a-11ea-3493-89487315feb7
# ╟─dfb7c6be-ee0d-11ea-194e-9758857f7b20
# ╟─e15ad330-ee0d-11ea-25b6-1b1b3f3d7888


text:
using Dates

struct ManagementMode
    name::String
    fo_length::Dates.Minute
end

function get_fo_length(mode::ManagementMode)
    return mode.fo_length
end

PSCOPF_MODE_1 = ManagementMode("mode_1", Dates.Minute(60))
PSCOPF_MODE_2 = ManagementMode("mode_2", Dates.Minute(60))
PSCOPF_MODE_3 = ManagementMode("mode_3", Dates.Minute(60))

ext: jl

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
# Thanks to Seth Bromberger for writing almost all of this code, and giving permission to use it
REGDIR = joinpath(DEPOT_PATH[1], "registries", "General")

version_in_range(v, s) = v in Pkg.Types.VersionSpec(s)
get_pkg_dir(s) = "$REGDIR/$(uppercase(s[1]))/$s"

# "Internal" utils and build the graph.

"""
Get all the dependencies of a package by
reading Deps.toml.
"""
function get_dep_dict(s)
    deps = Dict{String, Vector{String}}()

    d = get_pkg_dir(s)
    vrange = ""
    !isfile("$d/Deps.toml") && return deps
    for l in readlines("$d/Deps.toml")
        if startswith(l, "[")
            vrange = strip(l, ['[',']','\"'])
            # println("l = $l, vrange = $vrange")
            deps[vrange] = Vector{String}()
        else
            pkg = split(l, " ")[1]
            if pkg != ""
                push!(deps[vrange], pkg)
            end
        end
    end
    return deps
end

"""
Get the latest version of a package by
reading Versions.toml.
"""
function get_latest_version(s)
    d = get_pkg_dir(s)
    maxver = v"0.0.0"
    !isfile("$d/Versions.toml") && return maxver
    for l in readlines("$d/Versions.toml")
        if startswith(l, "[")
            v = VersionNumber(strip(l, ['[',']','\"']))
            if v > maxver
                maxver = v
            end
        end
    end
    return maxver
end

"""
Given a package name, return a vector
of all _direct_ depedencies.
"""
function get_deps(s)
    depdict = get_dep_dict(s)
    depset = Set{String}()
    maxver = get_latest_version(s)
    for (vrange, deps) in depdict
        # println("vrange = $vrange, deps = $deps")
        if version_in_range(maxver, vrange)
            union!(depset, deps)
        end
    end
    return collect(depset)
end

function make_pkg_list(r, omit_packages)
    f = "$r/Registry.toml"
    t = Pkg.TOML.parsefile(f)["packages"]

    omitset = Set(omit_packages)
    pkglist = Set(v["name"] for (_, v) in t)
    for p in pkglist
        # println("making $p")
        depset = Set(get_deps(p))
        union!(pkglist, depset)
    end
    setdiff!(pkglist, omitset)
    return sort(collect(pkglist))
end

function build_graph(r=REGDIR; omit_packages=[])
    pkglist = make_pkg_list(r, omit_packages)
    pkgrev = Dict((k, v) for (v, k) in enumerate(pkglist))

    g = MetaDiGraph(length(pkglist))
    for p in pkglist
        v = pkgrev[p]
        set_prop!(g, v, :name, p)
        deps = get_deps(p)
        setdiff!(deps, omit_packages)
        for d in deps
            w = pkgrev[d]
            add_edge!(g, v, w)
        end
    end
    set_indexing_prop!(g, :name)
    return (g, reverse(g), pkglist, pkgrev)
end

function __init__()
    global g, rev_g, pkglist, pkgrev = build_graph()
end

#===
"Customer-facing" functions.
===#

"""
See all n-th order dependents of a package.
"""
function get_dependents(s, n = 1)
    s_index = pkgrev[s]
    new_g = egonet(rev_g, s_index, n)
    dependents = [props(new_g, dependent)[:name] for dependent in unique(vertices(new_g)[2:end])]
end

"""
See all n-th order dependencies of a package.
"""
function get_dependencies(s, n = 1)
    s_index = pkgrev[s]
    new_g = egonet(g, s_index, n)
    dependents = [props(new_g, dependent)[:name] for dependent in unique(vertices(new_g)[2:end])]
end


text:
function test_gencode()
    index = gencode_load("GRCh37")
    genes = gencode_locate(index, "chr1", 12300)
    if length(genes)==0
        return false
    end
    g = genes[1]
    if g["gene"]!="DDX11L1" || g["number"]!=1 || g["type"]!="intron"
        return false
    end
    genes = gencode_genes(index, "TP53")
    g = genes[1]
    if g.name!="TP53" || g.chr!="chr17" || g.start_pos!=7565097 || g.end_pos!=7590856
        return false
    end
    return true
end

@test test_gencode()


text:
#  MIT License
#
#  Copyright (c) 2020 Michael J Simms. All rights reserved.
#
#  Permission is hereby granted, free of charge, to any person obtaining a copy
#  of this software and associated documentation files (the "Software"), to deal
#  in the Software without restriction, including without limitation the rights
#  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
#  copies of the Software, and to permit persons to whom the Software is
#  furnished to do so, subject to the following conditions:
#
#  The above copyright notice and this permission notice shall be included in all
#  copies or substantial portions of the Software.
#
#  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
#  SOFTWARE.

include("Peaks.jl")

using Pkg
Pkg.add("ArgParse")
using ArgParse
Pkg.add("CSV")
using CSV

function read_accelerometer_csv(filename::String)
    data = []
    data = CSV.read(filename)
    ts = data[1]
    x = data[2]
    y = data[3]
    z = data[4]
    ts, x, y, z
end

# Parses the command line arguments
function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--csv"
            help = "A CSV file with accelerometer data"
            arg_type = String
            default = "../data/pullups.csv"
        "--threshold"
            help = "Peaks must be above this threshold value to be considered"
            arg_type = Number
            default = 0.0
    end

    return parse_args(s)
end

parsed_args = parse_commandline()
ts, x, y, z = read_accelerometer_csv(parsed_args["csv"])
threshold = parsed_args["threshold"]

x_peaks = Peaks.find_peaks_over_threshold(x, threshold)
y_peaks = Peaks.find_peaks_over_threshold(y, threshold)
z_peaks = Peaks.find_peaks_over_threshold(z, threshold)

println("X-Axis Peaks")
for peak in x_peaks
    println("{ (", peak.left_trough.x, ", ", peak.left_trough.y, "), (", peak.right_trough.x, ", ", peak.right_trough.y, ")}")
end
println("Y-Axis Peaks")
for peak in y_peaks
    println("{ (", peak.left_trough.x, ", ", peak.left_trough.y, "), (", peak.right_trough.x, ", ", peak.right_trough.y, ")}")
end
println("Z-Axis Peaks")
for peak in z_peaks
    println("{ (", peak.left_trough.x, ", ", peak.left_trough.y, "), (", peak.right_trough.x, ", ", peak.right_trough.y, ")}")
end

max_stars_repo_name: msimms/PeakFinder

max_issues_repo_name: msimms/PeakFinder

max_forks_repo_name: msimms/PeakFinder


text:
@testset "724.find-pivot-index.jl" begin
    @test pivot_index([1, 7, 3, 6, 5, 6]) == 4
    @test pivot_index([1, 2, 3]) == 0
    @test pivot_index(Int[]) == 0
    @test pivot_index([-1, 1]) == 0
end


text:
#
# Function to read atomic coordinates from a PDB file.
#
function getpdbcoords!(pdbfile,atoms)
  file = open(pdbfile,"r")
  lines = readlines(file)
  local iatom = 0
  for line in lines
    data = split(line)
    if length(data) < 1
      continue
    end
    if data[1] == "ATOM" || data[1] == "HETATM"
      iatom = iatom + 1
      coor = [ parse(Float64,line[31:38]), parse(Float64,line[39:46]),  parse(Float64,line[47:54]) ]
      a = atoms[iatom]
      atoms[iatom] = Atom(a.index,a.residue,a.resid,a.name,a.resname,a.segname,a.type,a.element,a.chain,
                          a.charge,a.mass,a.eps,a.rmin,a.eps14,a.rmin14,a.backbone,
                          coor)
    end
  end
end


text:
module LowRankArithmetic

using Combinatorics, LinearAlgebra, UnPack, TSVD
import Base: +, -, *, size, Matrix, getindex, hcat, vcat, axes, round, broadcasted, BroadcastStyle
import LinearAlgebra: rank, adjoint, svd, qr

abstract type AbstractLowRankRepresentation end

# catch-all-cases
axes(LRA::AbstractLowRankRepresentation) = map(Base.oneto, size(LRA))
hadamard(A::AbstractLowRankRepresentation, B::AbstractMatrix) = Matrix(A) .* B
hadamard(A::AbstractMatrix, B::AbstractLowRankRepresentation) = Matrix(A) .* B
hadamard(A,B) = A .* B

# standard linear algebra broadcasts
broadcasted(::typeof(*), A::AbstractLowRankRepresentation, B::AbstractLowRankRepresentation) = hadamard(A,B)
broadcasted(::typeof(*), A::AbstractLowRankRepresentation, b::AbstractVector) = multiply_cols(A,b)
broadcasted(::typeof(*), A::AbstractLowRankRepresentation, b::Adjoint{<:Number, <:AbstractVector}) = multiply_rows(A, transpose(b))
broadcasted(::typeof(*), A::AbstractLowRankRepresentation, b::Transpose{<:Number, <:AbstractVector}) = multiply_rows(A, transpose(b))
broadcasted(::typeof(*), A::AbstractLowRankRepresentation, B::AbstractMatrix) = Matrix(A) .* B
broadcasted(::typeof(*), A::AbstractMatrix, B::AbstractLowRankRepresentation) = A .* Matrix(B)

broadcasted(::typeof(+), A::AbstractLowRankRepresentation, α::Number) = add_scalar(A, α)
broadcasted(::typeof(+), α::Number, A::AbstractLowRankRepresentation) = add_scalar(A, α)
broadcasted(::typeof(+), A::AbstractLowRankRepresentation, b::AbstractVector) = add_to_cols(A, b)
broadcasted(::typeof(+), A::AbstractLowRankRepresentation, b::Adjoint{<:Number, <:AbstractVector}) = add_to_rows(A, transpose(b))
broadcasted(::typeof(+), A::AbstractLowRankRepresentation, b::Transpose{<:Number, <:AbstractVector}) = add_to_rows(A, transpose(b))

broadcasted(::typeof(+), b::AbstractVector, A::AbstractLowRankRepresentation) = add_to_cols(A, b)
broadcasted(::typeof(+), b::Adjoint{<:Number, <:AbstractVector}, A::AbstractLowRankRepresentation) = add_to_rows(A, transpose(b))
broadcasted(::typeof(+), b::Transpose{<:Number, <:AbstractVector}, A::AbstractLowRankRepresentation) = add_to_rows(A, transpose(b))

broadcasted(::typeof(^), A::AbstractLowRankRepresentation, d::Int) = elpow(A, d)
broadcasted(::typeof(Base.literal_pow), ::typeof(^), A::AbstractLowRankRepresentation, ::Val{d}) where d = elpow(A, d)

include("orthonormalization.jl")
include("utils.jl")
include("two_factor_representation.jl")
include("svd_like_representation.jl")
end


text:
include("tile.jl")
include("annotations.jl")
include("statusbar.jl")
include("test4d.jl")
ImageView.view(img)


text:
# ---
# title: 34. Find First and Last Position of Element in Sorted Array
# id: problem34
# author: zhwang
# date: 2022-02-13
# difficulty: Medium
# categories: Array, Binary Search
# link: <https://leetcode.com/problems/find-first-and-last-position-of-element-in-sorted-array/description/>
# hidden: true
# ---
#
# Given an array of integers `nums` sorted in ascending order, find the starting
# and ending position of a given `target` value.
#
# If `target` is not found in the array, return `[-1, -1]`.
#
# **Follow up:**  Could you write an algorithm with `O(log n)` runtime
# complexity?
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [5,7,7,8,8,10], target = 8
#     Output: [3,4]
#
#
# **Example 2:**
#
#
#
#     Input: nums = [5,7,7,8,8,10], target = 6
#     Output: [-1,-1]
#
#
# **Example 3:**
#
#
#
#     Input: nums = [], target = 0
#     Output: [-1,-1]
#
#
#
#
# **Constraints:**
#
#   * `0 <= nums.length <= 105`
#   * `-109 <= nums[i] <= 109`
#   * `nums` is a non-decreasing array.
#   * `-109 <= target <= 109`
#
#
## @lc code=start
using LeetCode

function search_range(nums::Vector{Int}, target::Int)::Vector{Int}
    function get_left_border(target)
        left, right = 1, n
        while left <= right
            mid = (left + right) >> 1
            if nums[mid] >= target
                right = mid - 1
            else
                left = mid + 1
            end
        end
        return left
    end
    n = length(nums)
    left_border = get_left_border(target)
    (left_border > n || nums[left_border] != target) && return [-1, -1]
    return [left_border - 1, get_left_border(target + 1) - 2]
end
## @lc code=end


text:
# # Magnetic field direction
#
# In this example we'll look at the influence of the magnetic field direction on
# a nighttime and daytime ionosphere propagation path.
# Results will be compared to LWPC.
#
# ## Generate scenarios
#
# We'll use the JSON I/O functionality of LongwaveModePropagator to create
# the different magnetic field direction scenarios as a `BatchInput{ExponentialInput}`.
# We're doing this to ensure identical inputs are passed to LWPC.
# Although commented out in the `generate` function below, we locally saved the
# `BatchInput` as a JSON file and used it to feed a custom Python code
# which generates LWPC input files, runs LWPC, and saves the results in an HDF5 file.
#
# First load the packages we need.

using Dates, Printf
using HDF5
using Plots

using LongwaveModePropagator
using LongwaveModePropagator: buildrun, Progress, next!

# For convenience, we'll define global `OUTPUT_RANGES` at which the electric field
# will be computed, and the magnetic field dip and azimuth angles `B_DIPS` and `B_AZS`
# (converted below from degrees to radians).

const OUTPUT_RANGES = 0:5e3:3000e3
const B_DIPS = [90.0, 60, 60, 60, 60, 60, 60, 60, 60]
const B_AZS = [0.0, 0, 45, 90, 135, 180, 225, 270, 315]  ## vertical, N, E, S, W
nothing  #hide

# The `generate` function produces the [`BatchInput`](@ref) of the scenarios.
# If we wanted, we could comment out the JSON lines near the bottom of the function
# to also write the `BatchInput` to a file.
# That's not necessary here.

function generate(hp, β)
    batch = BatchInput{ExponentialInput}()
    batch.name = "Magnetic field tests"
    batch.description = "Varying magnetic field directions: vertical, N, E, S, W."
    batch.datetime = Dates.now()

    ## Constants
    frequency = 24e3
    output_ranges = collect(OUTPUT_RANGES)
    segment_ranges = [0.0]
    hprime = [hp]
    beta = [β]

    ## Ocean
    ground_epsr = [81]
    ground_sigma = [4.0]

    b_mag = [50e-6]
    b_dips = deg2rad.(B_DIPS)
    b_azs = deg2rad.(B_AZS)

    N = length(b_azs)
    inputs = Vector{ExponentialInput}(undef, N)
    for i in 1:N
        input = ExponentialInput()

        input.name = @sprintf("%d_%.0f_%.0f", i, b_dips[i], b_azs[i])
        input.description = "Wait ionosphere with ocean ground at 24 kHz."
        input.datetime = Dates.now()
        input.segment_ranges = segment_ranges
        input.hprimes = hprime
        input.betas = beta
        input.b_mags = b_mag
        input.b_dips = [b_dips[i]]
        input.b_azs = [b_azs[i]]
        input.ground_sigmas = ground_sigma
        input.ground_epsrs = ground_epsr
        input.frequency = frequency
        input.output_ranges = output_ranges

        inputs[i] = input
    end
    batch.inputs = inputs

    ## json_str = JSON3.write(batch)
    ##
    ## open("bfields.json","w") do f
    ##     write(f, json_str)
    ## end

    return batch
end

batch = generate(82.0, 0.6);

# ## Run the model
#
# Results from my custom LWPC script are saved in an
# [HDF5 file](@__REPO_ROOT_URL__/examples/bfields_lwpc.h5).
# For easy comparison, we will also save the results from the
# Longwave Mode Propagator to an HDF5 file.
# Although not strictly necessary for a small number inputs like we have here,
# using an HDF5 file rather than the JSON output file that is generated automatically by
# LongwaveModePropagator when a file is passed to `propagate` is more robust
# for thousands or tens of thousands of input scenarios.

function runlmp(inputs, outfile)

    h5open(outfile, "cw") do fid
        ## Create Batch attributes if they don't already exist
        fid_attrs = attributes(fid)
        haskey(fid_attrs, "name") || (fid_attrs["name"] = inputs.name)
        haskey(fid_attrs, "description") || (fid_attrs["description"] = inputs.description)
        haskey(fid_attrs, "datetime") || (fid_attrs["datetime"] = string(Dates.now()))

        if haskey(fid, "outputs")
            g = fid["outputs"]
        else
            g = create_group(fid, "outputs")
        end
    end

    PM = Progress(length(inputs.inputs), 5)
    for i in eachindex(inputs.inputs)
        ## If we've already run this, we can skip ahead
        complete = h5open(outfile, "r") do fid
            g = open_group(fid, "outputs")
            haskey(g, string(i)) ? true : false
        end
        complete && (next!(PM); continue)

        output = buildrun(inputs.inputs[i])

        h5open(outfile, "r+") do fid
            g = open_group(fid, "outputs")
            o = create_group(g, string(i))
            attributes(o)["name"] = output.name
            attributes(o)["description"] = output.description
            attributes(o)["datetime"] = string(Dates.now())
            o["output_ranges"] = output.output_ranges
            o["amplitude"] = output.amplitude
            o["phase"] = output.phase
        end
        next!(PM)
    end
end
nothing  #hide

#
# `examples_dir` below will need to be modified to the directory where
# `bfields_lmp.h5` is located and where you would like `bfields_lmp.h5` to be saved.

root_dir = dirname(dirname(pathof(LongwaveModePropagator)))
examples_dir = joinpath(root_dir, "examples")
lmpfile = joinpath(examples_dir, "bfields_lmp.h5")
runlmp(batch, lmpfile)

# ## Plots
#
# To more easily compare the amplitude and phase curves, we define a function
# to process the HDF5 files into arrays of amplitude and phase
# versus range and magnetic field azimuth.

function process(outputs)
    dist = read(outputs["outputs"]["1"]["output_ranges"])
    mask = indexin(OUTPUT_RANGES, dist)

    agrid = Array{Float64}(undef, length(OUTPUT_RANGES), length(B_AZS))
    pgrid = similar(agrid)

    for i in eachindex(B_AZS)
        o = outputs["outputs"][string(i)]

        amp = read(o["amplitude"])
        phase = read(o["phase"])

        agrid[:,i] = amp[mask]
        pgrid[:,i] = rad2deg.(phase[mask])
    end

    return agrid, pgrid
end

agrid, pgrid = h5open(lmpfile, "r") do o
    agrid, pgrid = process(o)
end

nothing  #hide

# Here is the amplitude plot itself.

labels = string.(trunc.(Int, B_AZS), "°")
labels[1] = "90°, "*labels[1]
labels[2] = "60°, "*labels[2]
labels[3:end] .= "      ".*labels[3:end]

colors = [palette(:phase, length(B_AZS))...]
pushfirst!(colors, RGB(0.0, 0, 0))

plot(OUTPUT_RANGES/1000, agrid;
     linewidth=1.5, palette=colors, colorbar=false,
     xlabel="range (km)", ylabel="amplitude (dB)",
     labels=permutedims(labels), legendtitle="  dip, az", legend=true)
#md savefig("magneticfield_amplitude.png"); nothing # hide
#md # ![](magneticfield_amplitude.png)

# The amplitude corresponding to each magnetic field azimuth
# (where 0° is along the propagation direction)
# is color-coded in the plot above, but a couple of the colors appear to be missing.
# This is because 0° and 180° are identical. There is no "north/south" dependence.
# Therefore 45°/135° and 315°/225° are also identical.
#
# On the other hand, there is a significant difference between propagation
# "east" (pinks) versus "west" (greens).
# This effect can be [observed in real life](https://doi.org/10.1016/0021-9169(58)90081-3).
#
# Now we'll plot the difference from LWPC.

lwpcfile = joinpath(examples_dir, "bfields_lwpc.h5")

lagrid, lpgrid = h5open(lwpcfile, "r") do o
    lagrid, bpgrid = process(o)
end

adifference = agrid - lagrid

plot(OUTPUT_RANGES/1000, adifference;
     linewidth=1.5, palette=colors, colorbar=false,
     xlabel="range (km)", ylabel="amplitude difference (dB)",
     labels=permutedims(labels), legendtitle="  dip, az", legend=true)
#md savefig("magneticfield_diff.png"); nothing # hide
#md # ![](magneticfield_diff.png)

# The two models are a very close match.
# Where there are large differences they occur because of slight misalignment of
# nulls with respect to range from the transmitter.
# The 315°/225° line has the worst match, although it is not clear if there is
# a particular cause for this outside of the general differences between the
# two models.

# ## Daytime ionosphere
#
# Does magnetic field direction have a different level of influence with daytime ionospheres?

lmpfile = joinpath(examples_dir, "bfields_daytime_lmp.h5")
runlmp(generate(72.0, 0.3), lmpfile)

agrid, pgrid = h5open(lmpfile, "r") do o
    agrid, pgrid = process(o)
end

plot(OUTPUT_RANGES/1000, agrid;
     linewidth=1.5, palette=colors, colorbar=false,
     xlabel="range (km)", ylabel="amplitude (dB)",
     labels=permutedims(labels), legendtitle="  dip, az", legend=true)
#md savefig("magneticfield_daytime_amplitude.png"); nothing # hide
#md # ![](magneticfield_daytime_amplitude.png)

# At these lower reflection altitudes the effect of Earth's magnetic field is decreased.


text:
using MLFS
using Test
using PrettyPrint
using Setfield
using MLFS.HM
using MLFS.TypedIO

# @testset "MLFS.jl" begin
#     it = Symbol(:i, sizeof(Int) * 8)
#     t = Surf.TArrow(Surf.TVar(:Int), Surf.TVar(:Int))
#     g = empty(GlobalTC)
#     l = empty(LocalTC)
#     l = @set l.typeEnv = l.typeEnv[
#         [
#         :Int => T(Nom(it)),
#         :Str => T(Nom(:str)),
#         :Bool => T(Nom(:bool)),
#         :f => Arrow(Nom(it), Nom(it)),
#         :Type => T(Nom(:Type))
#         ]
#     ]
#     l = @set l.symmap = l.symmap[
#             :f => :f_0xa1
#         ]
#     # println(inferType(g, l, t))
#     exp = @surf_expr f(1)
#     # pprintln(exp)
#     # pprintln(MLFS.inferExpr(g, l, exp)(NoProp))
#     # Write your tests here.

#     stmts = @surf_toplevel begin
#         g :: (a -> a) where a
#         g = x -> x

#         h :: (((a -> a) where a) -> (Int, Str))
#         h = f -> (f(1), f("1"))

#         id :: ((a -> a) where a)
#         id = x -> x

#         choose :: ((b -> b -> b) where b)
#         choose = x -> y -> x

#         choose_id = choose(id)
#         # choose_id.?(hint_choose_id)

#         J :: Fn[a, Fn[b, (a, b)] where b] where a
#         J = x -> y -> (x, y)

#         # c :: Fn[a, (Int, a)] where a
#         c = J(1)
#         # c.?(hint_c)

#         pair = (c(1), c("a")) # .?(hint_pair)

#         choose_id′ :: ((Int -> Int) -> (Int -> Int))
#         choose_id′ = (choose(id)) # .?(stronger_assumptions_1)

#         choose_id′ = (choose(id)) # .?(stronger_assumptions_2)

#         choose_id′ :: Fn[Fn[c, c], Fn[c, c]] where c
#         choose_id′ =
#             choose(id) # .?(stronger_assumptions_3)

#         choose_id′ :: ((Fn[a, a] where a) -> (Int -> Int))
#             choose_id′ = (choose(id)) # .?(stronger_assumptions_4)

#         choose_id′′ :: ((a where a) -> (Int -> Int))
#         choose_id′′ = choose_id′

#         choose_id′ :: ((Fn[a, a] where a) -> (Fn[a, a] where a))
#             choose_id′ = (choose(id)) # .?(stronger_assumptions_5)

#         choose_id′′ :: ((a where a) -> (Int -> Int))
#         choose_id′′ =
#             let choose_id_ :: ((Fn[a, a] where a) -> (Int -> Int)),
#                 choose_id_ = choose(id)

#                 choose_id_
#             end

#         # choose_id′ :: ((Int -> Int) -> (Int -> Int))
#         # choose_id′ = choose_id

#         int_functional :: Fn[Fn[Int, Int], Int]
#         int_functional = f -> f(1)

#         int_number = int_functional(id) # .?(app_functional)

#         Vec :: Type[NewType(Vec)]
#         Vec = Vec

#         mkVec :: (() -> Vec[A]) where A
#         mkVec =  (@julia begin
#             T -> T[]
#         end)(A)

#         vecLen :: (Vec[A] -> Int) where A
#         vecLen = @julia begin
#             (x::Vector) -> length(x)
#         end

#         # mkVec.?(vec)

#         imVec :: implicit[Vec[Int]]
#         imVec = mkVec(())
#         # imVec.?(imVec)

#         ctxVecLen :: Fn[implicit[Vec[Int]], Int]
#         ctxVecLen =  x -> vecLen(x)

#         (+) :: Fn[Int, Fn[Int, Int]]
#         (+) = @julia(+)

#         ctxOp :: Fn[implicit[Vec[Int]], Int]
#         ctxOp = vec -> ctxVecLen.?("implicit_instance")

#         ctxVecLen.?("implicit_meth")

#         Eq :: Type[NewType(Eq)]

#         makeEqList :: implicit[(implicit[Eq[a]] -> Eq[a]) where a]
#         makeEqList = @julia(nothing)

#         getEq :: (Eq[a] -> (a, a) -> Bool) where a
#         getEq = @julia(nothing)

#         (==) :: (implicit[Eq[a]] -> (a, a) -> Bool) where a
#         (==) = eq -> ab -> getEq(eq)(ab)

#         aVec :: Vec[Int]
#         aVec = mkVec(())

#         _ = (aVec == aVec).?(instance_not_found)

#         # List :: Type[NewType(List)]
#         # List = List

#         # Nil :: List[A] where A
#         # Nil = @julia begin
#         #     ()
#         # end

#         # Cons :: (A -> List[A] -> List[A]) where A
#         # Cons = (@julia begin
#         #     t -> x -> y -> (x :: t, y)
#         # end)(A)

#         # Cons.?(cons)
#         # Addable :: Type[NewType(Add1)]

#         # constructAddable :: ((A -> A -> A) -> Addable[A]) where A
#         # constructAddable = @julia begin
#         #     function constructAddable(a::Any)
#         #         a
#         #     end
#         # end

#         # getFieldA :: implict[Field[constructAddable[A], :plus, (A -> A -> A)] where A]
#         # getFieldA = @julia begin
#         #     function getFieldA(a::Any)
#         #         a
#         #     end
#         # end

#         # add :: Fn[implicit[Addable[A]], A -> A -> A] where A
#         # add = addable -> l -> r ->
#         #     addable.plus(l)(r)

#     end

#     for f in MLFS.inferDecls(g, l, stmts)[1]
#         f()
#     end

#     show_hints(g)

# end

# include("testim.jl")
# open("testim_gen.jl") do f
#     println(read(f, String))
# end

include("a.jl")


text:
export vec_data

abstract type AbstractDerivative end

function _isapprox(x, y; kwargs...)
    kx = propertynames(x)
    ky = propertynames(y)
    all(kx .== ky) || return false
    for f=kx
        _isapprox(getproperty(x, f), getproperty(y, f); kwargs...) || return false
    end
    return true
end

_isapprox(x::AbstractArray, y::AbstractArray; kwargs...) =
    isapprox(x, y; kwargs...)


text:
import GLAbstraction: FrameBuffer, VertexArray, Buffer, Program, BufferAttachmentInfo
import GLAbstraction: textureformat_from_type_sym, getfirst, gluniform, clear!

#----------------GLAbstraction-------------------------#
function default_framebuffer(fb_size)
    return FrameBuffer(fb_size, DepthStencil{GLAbstraction.Float24,N0f8}, RGBA{N0f8},
                       Vec{2,GLushort}, RGBA{N0f8})
end

function clear!(fbo::FrameBuffer, color::RGBA)
    return clear!(fbo, (color.r, color.g, color.b, color.alpha))
end

const fullscreen_pos = [Vec3f0(-1, 1, 0), Vec3f0(-1, -1, 0), Vec3f0(1, 1, 0),
                        Vec3f0(1, -1, 0)]

const fullscreen_uv = [Vec2f0(0, 1), Vec2f0(0, 0), Vec2f0(1, 1), Vec2f0(1, 0)]

#REVIEW: not used
function glenum2julia(x::UInt32)
    x == GL_FLOAT && return Float32
    x == GL_FLOAT_VEC3 && return Vec3f0
    x == GL_FLOAT_VEC4 && return Vec4f0
    x == GL_FLOAT_MAT2 && return Mat2f0
    x == GL_FLOAT_MAT3 && return Mat3f0
    x == GL_FLOAT_MAT4 && return Mat4f0

    x == GL_DOUBLE && return f64
    x == GL_DOUBLE_VEC3 && return Vec3{f64}
    x == GL_DOUBLE_VEC4 && return Vec4{f64}
    x == GL_DOUBLE_MAT2 && return Mat2{f64}
    x == GL_DOUBLE_MAT3 && return Mat3{f64}
    x == GL_DOUBLE_MAT4 && return Mat4{f64}

    x == GL_INT && return i32
    x == GL_INT_VEC2 && return Vec2{i32}
    x == GL_INT_VEC3 && return Vec3{i32}
    x == GL_INT_VEC4 && return Vec4{i32}

    x == GL_UNSIGNED_INT && return u32
    x == GL_UNSIGNED_INT_VEC2 && return Vec2{u32}
    x == GL_UNSIGNED_INT_VEC3 && return Vec3{u32}
    return x == GL_UNSIGNED_INT_VEC4 && return Vec4{u32}
end

function mergepop!(d1, d2)
    t = SymAnyDict()

    d = isempty(d2) ? SymAnyDict() : Dict(d2)
    for (key, val) in d1
        t[key] = pop!(d, key, val)
    end
    d2 = [d...]
    return t
end

Base.ndims(::Type{<:Colorant}) = 1
Base.size(::Type{<:Number}) = 1

#this stuff is here because of world age stuff I think
function gluniform(location::Integer, x::Mat4{Float32})
    return glUniformMatrix4fv(location, 1, GL_FALSE, reinterpret(Float32, [x]))
end

function uniformfunc(typ::DataType, dims::Tuple{Int})
    return Symbol(string("glUniform", first(dims), GLAbstraction.opengl_postfix(typ)))
end

function uniformfunc(typ::DataType, dims::Tuple{Int,Int})
    M, N = dims
    return Symbol(string("glUniformMatrix", M == N ? "$M" : "$(M)x$(N)",
                         opengl_postfix(typ)))
end

function gluniform(location::Integer, x::FSA) where {FSA}
    GLAbstraction.glasserteltype(FSA)
    xref = [x]
    return gluniform(location, xref)
end

@generated function gluniform(location::Integer, x::Vector{FSA}) where {FSA}
    GLAbstraction.glasserteltype(eltype(FSA))
    func = uniformfunc(eltype(FSA), size(FSA))
    callexpr = if ndims(FSA) == 2
        :($func(location, length(x), GL_FALSE, xref))
    else
        :($func(location, length(x), xref))
    end
    quote
        xref = reinterpret(eltype(FSA), x)
        $callexpr
    end
end

gluniform(location::Integer, x::Int) = gluniform(location, GLint(x))
gluniform(location::Integer, x::UInt) = gluniform(location, GLuint(x))
gluniform(location::Integer, x::AbstractFloat) = gluniform(location, GLfloat(x))

#----------------- Composite OpenGL Calls -------------#
const none = !any

function glEnableDepth(f::Function)
    glEnable(GL_DEPTH_TEST)
    if f == all
        glDepthFunc(GL_ALWAYS)
    elseif f == none
        glDepthFunc(GL_NEVER)
    elseif f == <
        glDepthFunc(GL_LESS)
    elseif f == <=
        glDepthFunc(GL_LEQUAL)
    elseif f == >
        glDepthFunc(GL_GREATER)
    elseif f == >=
        glDepthFunc(GL_GEQUAL)
    elseif f == ==
        glDepthFunc(GL_EQUAL)
    elseif f == !=
        glDepthFunc(GL_NOTEQUAL)
    else
        @error "Function $f has no OpenGL DepthFunc equivalent."
    end
end

glDisableDepth() = glDisable(GL_DEPTH_TEST)

function glEnableCullFace(s::Symbol)
    glEnable(GL_CULL_FACE)
    if s == :front
        glCullFace(GL_FRONT)
    elseif s == :back
        glCullFace(GL_BACK)
    elseif s == :front_back
        glCullFace(GL_FRONT_AND_BACK)
    else
        @error "Symbol $s has no OpenGL CullFace equivalent."
    end
end

glDisableCullFace() = glDisable(GL_CULL_FACE)

function generate_buffers(program::Program, divisor::GLint; name_buffers...)
    buflen  = 0
    buffers = BufferAttachmentInfo[]
    for (name, val) in pairs(name_buffers)
        loc = attribute_location(program, name)
        if loc != INVALID_ATTRIBUTE
            buflen = buflen == 0 ? length(val) : buflen
            vallen = length(val)
            if vallen == buflen
                push!(buffers,
                      BufferAttachmentInfo(name, loc, Buffer(val; usage = GL_DYNAMIC_DRAW),
                                           divisor))
            elseif !isa(val, Vector)
                push!(buffers,
                      BufferAttachmentInfo(name, loc,
                                           Buffer(fill(val, buflen);
                                                  usage = GL_DYNAMIC_DRAW), divisor))
            end
        end
    end
    return buffers
end
## StaticArrays
const Area   = Rect2D
const X_AXIS = Vec3f0(1.0f0, 0.0, 0.0)
const Y_AXIS = Vec3f0(0.0, 1.0f0, 0.0)
const Z_AXIS = Vec3f0(0.0, 0.0, 1.0f0)

GeometryBasics.Vec3f0(v::Vec2, nr = 0.0) = Vec3f0(v..., nr)
GeometryBasics.Vec3f0(v::Vec3) = Vec3f0(v...)
GeometryBasics.Vec4f0(v::Vec2, n1 = 0.0, n2 = 0.0) = Vec4f0(v..., n1, n2)
GeometryBasics.Vec4f0(v::Vec3, n1 = 0.0) = Vec4f0(v..., n1)

#----------------GeometryBasics-------------------------#
struct Cone{T} <: GeometryPrimitive{3,T}
    origin::Point3{T}
    extremity::Point3{T}
    r::T
end

GeometryBasics.origin(c::Cone{T}) where {T}    = c.origin
GeometryBasics.extremity(c::Cone{T}) where {T} = c.extremity
GeometryBasics.radius(c::Cone{T}) where {T}    = c.r
GeometryBasics.height(c::Cone{T}) where {T}    = norm(c.extremity - c.origin)
GeometryBasics.direction(c::Cone{T}) where {T} = (c.extremity .- c.origin) ./ GeometryBasics.height(c)

function GeometryBasics.rotation(c::Cone{T}) where {T}
    d3 = direction(c)
    u = GeometryBasics.@SVector [d3[1], d3[2], d3[3]]
    if abs(u[1]) > 0 || abs(u[2]) > 0
        v = GeometryBasics.@MVector [u[2], -u[1], T(0)]
    else
        v = GeometryBasics.@MVector [T(0), -u[3], u[2]]
    end
    normalize!(v)
    w = GeometryBasics.@SVector [u[2] * v[3] - u[3] * v[2], -u[1] * v[3] + u[3] * v[1],
                                 u[1] * v[2] - u[2] * v[1]]
    return hcat(v, w, u)
end

function GeometryBasics.decompose(PT::Type{Point{3,T}}, c::Cone, resolution = 30) where {T}
    isodd(resolution) && (resolution = 2 * div(resolution, 2))
    resolution = max(8, resolution)
    nbv = div(resolution, 2)
    M = GeometryBasics.rotation(c)
    h = GeometryBasics.height(c)
    position = 1
    vertices = Vector{PT}(undef, nbv + 2)
    for j in 1:nbv
        phi = T((2π * (j - 1)) / nbv)
        vertices[j] = PT(M * Point{3,T}(c.r * cos(phi), c.r * sin(phi), 0)) + PT(c.origin)
    end
    vertices[end-1] = PT(c.origin)
    vertices[end] = PT(c.extremity)
    return vertices
end

function GeometryBasics.decompose(::Type{FT}, c::Cone,
                                  resolution = 30) where {FT<:GeometryBasics.AbstractFace}
    isodd(resolution) && (resolution = 2 * div(resolution, 2))
    resolution = max(8, resolution)
    nbv = div(resolution, 2)
    indexes = Vector{FT}(undef, resolution)
    index = 1
    for j in 1:nbv-1
        indexes[index] = (index, nbv + 1, index + 1)
        indexes[index+nbv] = (index, index + 1, nbv + 2)
        index += 1
    end
    # indexes[index] = (index,  index + 1, nbv+1)
    # indexes[index+nbv] = (index,  nbv+2, index + 1)
    indexes[nbv] = (1, nbv + 2, nbv)
    indexes[end] = (1, nbv, nbv + 1)
    # indexes[end] = (1, 1, 1)
    return indexes
end

struct Arrow{T} <: GeometryPrimitive{3,T}
    origin::Point3{T}
    extremity::Point3{T}
    r::T #cylinder ratio
    length_ratio::T # height = height_cylinder * (1+length_ratio)
    radius_ratio::T # cone_radius = radius_ratio * r
end

GeometryBasics.origin(c::Arrow)    = c.origin
GeometryBasics.extremity(c::Arrow) = c.extremity
GeometryBasics.radius(c::Arrow)    = c.r
GeometryBasics.height(c::Arrow)    = norm(c.extremity - c.origin)
GeometryBasics.direction(c::Arrow) = (c.extremity .- c.origin) ./ GeometryBasics.height(c)

function GeometryBasics.Cylinder(c::Arrow)
    return Cylinder(c.origin, direction(c) * height(c) / (1 + c.length_ratio), c.r)
end

function Cone(c::Arrow)
    return Cone(direction(c) * height(c) / (1 + c.length_ratio), c.extremity,
                c.r * c.radius_ratio)
end

function GeometryBasics.decompose(PT::Type{<:Point}, c::Arrow, resolution = 30)
    return [decompose(PT, Cylinder(c), resolution); decompose(PT, Cone(c), resolution)]
end

function GeometryBasics.decompose(::Type{FT}, c::Arrow,
                                  resolution = 30) where {FT<:GeometryBasics.AbstractFace}
    cylinder_indices = decompose(FT, Cylinder(c), resolution)
    last_id = maximum(maximum.(cylinder_indices))
    cone_indices = [id .+ last_id for id in decompose(FT, Cone(c), resolution)]
    return [cylinder_indices; cone_indices]
end

# const INSTANCEpD_MESHES = Dict{Symbol, BasicMesh}()
#
# struct InstancedAttributeMesh{D, T, FD, FT, AT <: NamedTuple} <: AbstractGlimpseMesh
#     basic      ::BasicMesh{D, T, FD, FT}
#     attributes ::Vector{AT}
# end
#
# function InstancedAttributeMesh(basic_symbol::Symbol, attributes::NamedTuple)
#     if !haskey(INSTANCED_MESHES, basic_symbol)
#         error("No instanced mesh `$basic_symbol` found.")
#     end
#     return InstancedAttributeMesh(INSTANCED_MESHES[basic_symbol], basic_symbol, attributes)
# end
#
# InstancedAttributeMesh(basic_symbol::Symbol, attributes::NamedTuple, geometry, args...) =
#     InstancedAttributeMesh(basic_symbol, attributes, BasicMesh(geometry, args...))
#
# InstancedAttributeMesh(basic_symbol::Symbol, geometry, args...; attributes...) =
#     InstancedAttributeMesh(basic_symbol, NamedTuple{keys(attributes)}(values(attributes)), BasicMesh(geometry, args...))
#
# basicmesh(mesh::InstancedAttributeMesh) = mesh.basic

#----------------------GLFW----------------------------#
glfw_destroy_current_context() = GLFW.DestroyWindow(GLFW.GetCurrentContext())

"""
Standard window hints for creating a plain context without any multisampling
or extra buffers beside the color buffer
"""
const GLFW_DEFAULT_WINDOW_HINTS = [(GLFW.SAMPLES, 0), (GLFW.DEPTH_BITS, 32),
                                   (GLFW.ALPHA_BITS, 8), (GLFW.RED_BITS, 8),
                                   (GLFW.GREEN_BITS, 8), (GLFW.BLUE_BITS, 8),
                                   (GLFW.STENCIL_BITS, 0), (GLFW.AUX_BUFFERS, 0),
                                   (GLFW.SCALE_TO_MONITOR, 1)]

function glfw_standard_screen_resolution()
    return GLFW.GetPrimaryMonitor() |>
           GLFW.GetMonitorPhysicalSize |>
           values .|>
           x -> div(x, 1)
end

## Colorutils
import ColorTypes: RGBA, Colorant, RGB
ColorTypes.RGBA(x::T) where {T<:Real} = RGBA{T}(x, x, x, x)
ColorTypes.RGBA{T}(x) where {T<:Real} = RGBA{T}(T(x), T(x), T(x), T(x))

Base.length(::Type{<:RGBA}) = 4
Base.size(x::Type{<:Colorant}) = (length(x),)
Base.size(x::Type{<:RGB{Float32}}) = (3,)
Base.ndims(x::Type{Colorant}) = 1
Base.convert(::Type{RGB{Float32}}, x::Union{AbstractVector, <:NTuple{3}}) = RGB{Float32}(x...)

const RGBAf0 = RGBA{Float32}
const RGBf0 = RGB{Float32}
const BLUE = RGBf0(0.0, 0.0, 1.0)
const GREEN = RGBf0(0.0, 1.0, 0.0)
const RED = RGBf0(1.0, 0.0, 0.0)
const BLACK = RGBf0(0.0, 0.0, 0.0)
const DEFAULT_COLOR = RGBf0(0.0, 0.4, 0.8)
const GREY = RGBf0(0.4, 0.4, 0.4)


text:
# Tests the vectorized version of cel:
#include("cel_bulirsch.jl")
# include("../src/cel_bulirsch.jl")
#if VERSION >= v"0.7"
#  using Test
#else
#  using Base.Test
#end

import Limbdark: cel_bulirsch

# Randomizer seed
using Random
Random.seed!(42)

# Function which tests cel_bulirsch both numerically, one-by-one and vectorized:
function test_cel!(kc,p,a,b)
k2 = 1.0-kc^2
# Evaluate all three at once:
ell2 = cel_bulirsch(k2,kc,p,a[1],a[2],a[3],b[1],b[2],b[3])
ell1 = zeros(3)
nphi = 100000
if p < 0
#  Make same transformations as in cel_bulirsch for negative p:
  p0 = sqrt((kc*kc-p)/(1-p))
  a0 = (a[1]-b[1])/(1-p)
  b0 = -(b[1]-a[1]*p)/(1-p)^2*(1-kc^2)/p0+a0*p0
  a[1] = a0
  b[1] = b0*p0
  p = p0*p0
end
# Numerically compute the integrals:
dphi = .5*pi/nphi
phi = linearspace(0.5*dphi,pi/2-.5*dphi,nphi)
ell3 = zeros(3); cphi2 = cos.(phi).^2; sphi2=sin.(phi).^2
den = dphi./sqrt.(cphi2+kc*kc*sphi2)
ell3[1] = sum((a[1]*cphi2+b[1]*sphi2)./(cphi2+p*sphi2).*den)
ell3[2] = sum((a[2]*cphi2+b[2]*sphi2).*den)
ell3[3] = sum((a[3]*cphi2+b[3]*sphi2).*den)

# First, evaluate one by one:
ell1[1] = cel_bulirsch(k2,kc,p,a[1],b[1])
ell1[2] = cel_bulirsch(k2,kc,1.0,a[2],b[2])
ell1[3] = cel_bulirsch(k2,kc,1.0,a[3],b[3])

return ell1,ell2,ell3
end

@testset "cel_bulirsch" begin
  # Test negative p values:
  p = -rand(); a=rand(3); b=rand(3)
  kc = rand()
  ell1,ell2,ell3 = test_cel!(kc,p,a,b)
  println("ell1: ",ell1)
  println("ell2: ",ell2)
  println("ell3: ",ell3)
  for i=1:3
    @test isapprox(ell1[i],ell2[i])
#    @test isapprox(ell1[i],ell3[i],atol=1e-6)
    @test isapprox(ell1[i],ell3[i])
  end
  # Test positive p values:
  p = rand(); a=rand(3); b=rand(3)
  kc = rand()
  ell1,ell2,ell3 = test_cel!(kc,p,a,b)
  println("ell1: ",ell1)
  println("ell2: ",ell2)
  println("ell3: ",ell3)
  for i=1:3
    @test isapprox(ell1[i],ell2[i])
#    @test isapprox(ell1[i],ell3[i],atol=1e-6)
    @test isapprox(ell1[i],ell3[i])
  end
  # Test small kc values:
  p = rand(); a=rand(3); b=rand(3)
  kc = 1e-4
  ell1,ell2,ell3 = test_cel!(kc,p,a,b)
  println("ell1: ",ell1)
  println("ell2: ",ell2)
  println("ell3: ",ell3)
  for i=1:3
    @test isapprox(ell1[i],ell2[i])
#    @test isapprox(ell1[i],ell3[i],atol=1e-6)
    @test isapprox(ell1[i],ell3[i])
  end
  # Test large kc values:
  p = rand(); a=rand(3); b=rand(3)
  kc = 1.0-1e-4
  ell1,ell2,ell3 = test_cel!(kc,p,a,b)
  println("ell1: ",ell1)
  println("ell2: ",ell2)
  println("ell3: ",ell3)
  for i=1:3
    @test isapprox(ell1[i],ell2[i])
#    @test isapprox(ell1[i],ell3[i],atol=1e-6)
    @test isapprox(ell1[i],ell3[i])
  end
end


text:
include("03_one_site_state_custom.jl")

println("
#######################################################
# Tutorial 4                                          #
#                                                     #
# 1-site operator                                     #
#######################################################
")

#
# Priming basics
#

# Operator as a linear map:
#
# O: ℂᴺ → ℂᴺ
#
# Consider ℂ² → ℂ²
#
# Decorations and labels of the Hilbert space
# like primes help label the input and output spaces
# to make them unique.
#
# Convenient way to specify a copy
# of an Index that lives in contravariant/
# covariant (bra and ket) spaces.
@show prime(i) ≠ i
# `i'` is just compact notation for `prime(i)`
@show prime(i) == i'
@show i' ≠ i

#
# Making operators
#

Z = ITensor(i', dag(i))
Z[i' => 1, i => 1] = 1
Z[i' => 2, i => 2] = -1

Z = ITensor([1 0; 0 -1], i', dag(i))

X = ITensor([0 1; 1 0], i', dag(i))

# Predefined operators:
Z = op("Z", i)

X = op("X", i)

# Apply the operator as a map
XZp = X * Zp

@visualize X * Zp

# Apply the operator as a map
# of the input vector.
XZp = X * Zp

# Take an inner product?
# dag(Zm) * XZp
# No, outer product!

@show XZp ≠ Zm
@show XZp == Zm'

# Take an inner product:
# Need to prime with:
@show (dag(Zm)' * XZp)[]
@show inner(Zm', XZp)
@show inner(Zp', XZp)

# Take an inner product
# We can also use `noprime`:
@show (dag(Zm) * noprime(XZp))[]
@show inner(Zm, noprime(XZp))
@show inner(Zp, noprime(XZp))

# `apply` is shorthand for applying the operator
# and mapping it back to the original space:
XZp = apply(X, Zp)

# Same as `noprime(X * Zp)`:
@show XZp == noprime(X * Zp)

@show inner(Zm, apply(X, Zp))

# Compute expectation values/matrix elements
@show (dag(Zm)' * X * Zp)[]
@show inner(Zm', X, Zp)

@visualize dag(Zm)' * X * Zp

@visualize dag(Zm)' * X * Zp backend="Makie"

@show inner(Zp', X, Zm)

@visualize dag(Zp)' * X * Zm

# More matrix elements
@show inner(Xp', X, Xp)
@show inner(Xm', X, Xm)


text:
using Test
using CBOR
using DataStructures
import CBOR: Tag, decode, encode, SmallInteger, UndefLength

# Taken (and modified) from Appendix A of RFC 7049

two_way_test_vectors = [
    SmallInteger(0) => hex2bytes("00"),
    SmallInteger(1) => hex2bytes("01"),
    SmallInteger(10) => hex2bytes("0a"),
    SmallInteger(23) => hex2bytes("17"),
    SmallInteger(24) => hex2bytes("1818"),
    SmallInteger(25) => hex2bytes("1819"),

    UInt8(100) => hex2bytes("1864"),
    UInt16(1000) => hex2bytes("1903e8"),
    UInt32(1000000) => hex2bytes("1a000f4240"),
    UInt64(1000000000000) => hex2bytes("1b000000e8d4a51000"),
    # UInt128(18446744073709551615) => hex2bytes("1bffffffffffffffff"),
    # Int128(-18446744073709551616) => hex2bytes("3bffffffffffffffff"),
    SmallInteger(-1) => hex2bytes("20"),
    SmallInteger(-10) => hex2bytes("29"),
    Int8(-100) => hex2bytes("3863"),
    Int16(-1000) => hex2bytes("3903e7"),

    0.0f0 => hex2bytes("fa00000000"),
    -0.0f0 => hex2bytes("fa80000000"),
    1.0f0 => hex2bytes("fa3f800000"),
    1.1 => hex2bytes("fb3ff199999999999a"),
    1.5f0 => hex2bytes("fa3fc00000"),
    65504f0 => hex2bytes("fa477fe000"),
    100000f0 => hex2bytes("fa47c35000"),
    Float32(3.4028234663852886e+38) => hex2bytes("fa7f7fffff"),
    1.0e+300 => hex2bytes("fb7e37e43c8800759c"),
    Float32(5.960464477539063e-8) => hex2bytes("fa33800000"),
    Float32(0.00006103515625) => hex2bytes("fa38800000"),
    -4f0 => hex2bytes("fac0800000"),
    -4.1 => hex2bytes("fbc010666666666666"),

    false => hex2bytes("f4"),
    true => hex2bytes("f5"),
    nothing => hex2bytes("f6"),
    Undefined() => hex2bytes("f7"),

    Tag(0, "2013-03-21T20:04:00Z") =>
        hex2bytes("c074323031332d30332d32315432303a30343a30305a"),
    Tag(1, SmallInteger(1363896240)) => hex2bytes("c11a514b67b0"),
    Tag(1, 1363896240.5) => hex2bytes("c1fb41d452d9ec200000"),
    Tag(23, hex2bytes("01020304")) => hex2bytes("d74401020304"),
    Tag(24, hex2bytes("6449455446")) => hex2bytes("d818456449455446"),
    Tag(32, "http://www.example.com") =>
        hex2bytes("d82076687474703a2f2f7777772e6578616d706c652e636f6d"),

    UInt8[] => hex2bytes("40"),
    hex2bytes("01020304") => hex2bytes("4401020304"),

    "" => hex2bytes("60"),
    "a" => hex2bytes("6161"),
    "IETF" => hex2bytes("6449455446"),
    "\"\\" => hex2bytes("62225c"),
    "\u00fc" => hex2bytes("62c3bc"),
    "\u6c34" => hex2bytes("63e6b0b4"),

    [] => hex2bytes("80"),
    SmallInteger[1, 2, 3] => hex2bytes("83010203"),
    [SmallInteger(1), SmallInteger[2, 3], SmallInteger[4, 5]] => hex2bytes("8301820203820405"),
    SmallInteger[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25] =>
        hex2bytes("98190102030405060708090a0b0c0d0e0f101112131415161718181819"),

    Dict() => hex2bytes("a0"),
    OrderedDict(SmallInteger(1)=>SmallInteger(2), SmallInteger(3)=>SmallInteger(4)) => hex2bytes("a201020304"),
    OrderedDict("a"=>SmallInteger(1), "b"=>SmallInteger[2, 3]) => hex2bytes("a26161016162820203"),
    OrderedDict("a"=>"A", "b"=>"B", "c"=>"C", "d"=>"D", "e"=>"E") =>
        hex2bytes("a56161614161626142616361436164614461656145"),

    ["a", Dict("b"=>"c")] => hex2bytes("826161a161626163")
]

# Some annoying problems with Small integers and indefinite lenght arrays
# for round tripping
cbor_equal(a, b) = a == b
cbor_equal(a::Vector{String}, b::String) = join(a, "") == b
cbor_equal(a::Vector{Vector{UInt8}}, b::Vector{UInt8}) = vcat(a...) == b
function cbor_equal(a::AbstractDict, b::AbstractDict)
    for (ka, kb) in zip(keys(a), keys(b))
        ka == kb || return false
    end
    for (ka, kb) in zip(values(a), values(b))
        ka == kb || return false
    end
    return true
end

#=
The problem is, we want to preserver Julia types for non Base types that directly
map to basic CBOR protocol types. So we can't define encode(io::IO, x::AbstractDict)
since that would mean we can't preserver the type of any custom dict type.
But, since the CBOR protocol is expecting ordered dicts, Julia's default dict type
is unordered. I think it still makes more sense to use Julia's dict type instead of
taking the dependencies on DataStructures. But to pass the byte test, we need
an ordered dict type, so we overload it just here!
=#

function CBOR.encode(io::IO, x::OrderedDict)
    CBOR.encode_length(io::IO, CBOR.TYPE_5, x)
    for (key, value) in x
        encode(io, key)
        encode(io, value)
    end
end
function CBOR.decode_internal(io::IO, ::Val{CBOR.TYPE_5})
    return OrderedDict(CBOR.decode_ntimes(io) do io
        CBOR.decode_internal(io) => CBOR.decode_internal(io)
    end...)
end

@testset "two way" begin
    for (data, bytes) in two_way_test_vectors
        @test cbor_equal(data, decode(encode(data)))
        @test isequal(bytes, encode(data))
        @test cbor_equal(data, decode(bytes))
    end
end

bytes_to_data_test_vectors = [
    hex2bytes("fa7fc00000") => NaN32,
    hex2bytes("fa7f800000") => Inf32,
    hex2bytes("faff800000") => -Inf32,
    hex2bytes("fb7ff8000000000000") => NaN,
    hex2bytes("fb7ff0000000000000") => Inf,
    hex2bytes("fbfff0000000000000") => -Inf
]

@testset "bytes to data" begin
    for (bytes, data) in bytes_to_data_test_vectors
        @test isequal(data, decode(bytes))
    end
end

iana_test_vector = [
    BigInt(18446744073709551616) => hex2bytes("c249010000000000000000"),
    BigInt(-18446744073709551617) => hex2bytes("c349010000000000000000")
]

@testset "BigInt iana" begin
    for (data, bytes) in iana_test_vector
        @test isequal(bytes, encode(data))
        @test isequal(data, decode(bytes))
    end
end

indef_length_coll_test_vectors = [
    ["Hello", " ", "world"] =>
        hex2bytes("7f6548656c6c6f612065776f726c64ff"),
    Vector{UInt8}.(["Hello", " ", "world"]) =>
        hex2bytes("5f4548656c6c6f412045776f726c64ff"),

    [SmallInteger(1), 2.3, "Twiddle"] =>
        hex2bytes("9f01fb40026666666666666754776964646c65ff"),

    OrderedDict(SmallInteger(1)=>SmallInteger(2), 3.2=>"3.2") =>
        hex2bytes("bf0102fb400999999999999a63332e32ff")
]

@testset "ifndef length collections" begin
    @testset "basic types" begin
        @test UInt8[1] == decode(encode(UndefLength(UInt8[1])))
        @test "hi" == decode(encode("hi"))
        @test [1, "2", UInt8[1]] == decode(encode([1, "2", UInt8[1]]))
        @test Dict("a" => 2) == decode(encode(Dict("a" => 2)))
    end

    for (data, bytes) in indef_length_coll_test_vectors
        @test isequal(bytes, encode(UndefLength(data)))
        @test cbor_equal(data, decode(bytes))
    end
end

#=
From the docs about undef length byte strings
5F              -- Start indefinite-length byte string
     44           -- Byte string of length 4
        aabbccdd  -- Bytes content
     43           -- Byte string of length 3
        eeff99    -- Bytes content
     FF           -- "break"

  After decoding, this results in a single byte string with seven
  bytes: 0xaabbccddeeff99.
=#
@testset "undef length bytestring" begin
    @test decode(hex2bytes("5f44aabbccdd43eeff99FF")) == hex2bytes("aabbccddeeff99")
end

# tests from the readme
function producer(ch::Channel)
    for i in 1:10
        put!(ch,i*i)
    end
end
@testset "indefinite length readme" begin
    iter = Channel(producer)

    @test ((1:10) .* (1:10)) == decode(encode(UndefLength(iter)))
end

function cubes(ch::Channel)
    for i in 1:10
        put!(ch,i)       # key
        put!(ch,i*i*i)   # value
    end
end

@testset "indefinite length readme Dict" begin
    bytes = encode(UndefLength{Pair}(Channel(cubes)))
    @test Dict(zip(1:10, (1:10) .^ 3)) == decode(bytes)
end

function producer(ch::Channel)
    for c in ["F", "ire", " ", "and", " ", "Blo", "od"]
        put!(ch, c)
    end
end

@testset "indefinite length readme String" begin
    bytes = encode(UndefLength{String}(Channel(producer)))
    @test decode(bytes) == "Fire and Blood"
end


text:
include("types.jl")
include("score.jl")
include("io.jl")
@load "library.jld2"

max_stars_repo_name: JuliaTagBot/play

max_issues_repo_name: JuliaTagBot/play

max_forks_repo_name: JuliaTagBot/play


text:
# import Base.String
# function String(envelope::AbstractEnvelope)
#     io = IOBuffer()
#     serialize(io,envelope)
#     String(take!(io))
# end

### I could create a Voucher and then take out valuable info
# """
# Verifies that the signature and calculates id of the public key.
# """
# function unwrap(envelope::Envelope,notary::Notary)
#     id = notary.verify("$(envelope.data)",envelope.signature)
#     return envelope.data, (envelope.uuid,id)
# end

# struct Envelope <: AbstractEnvelope
#     uuid::UUID
#     data
#     signature
# end

#Envelope(str::AbstractString) = deserilize(IOBuffer(str))

### The Certificates

# struct ID <: AbstractID
#     name
#     date # of birth
#     id # The id of the public key. One uses the hash function from the Community.
# end

# Certificate could also be made as Envelope{ID}.
# struct Certificate <: AbstractEnvelope
#     uuid::UUID ### Now it is the deme
#     data::AbstractID
#     signature#::AbstractSignature
# end

#Certificate(id::AbstractID,signer::AbstractSigner) = Certificate(signer.uuid,id,signer.sign("$id"))
#Certificate(str::AbstractString) = deserilize(IOBuffer(str)) # We could do a little better here

# function register(uuid::UUID,certificate::Certificate)
#     com = community(uuid)
#     com.register(certificate)
# end

#export Certificate, Envelope, unwrap