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

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

# ╔═╡ 65780f00-ed6b-11ea-1ecf-8b35523a7ac0
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add([
			Pkg.PackageSpec(name="Images", version="0.22.4"),
			Pkg.PackageSpec(name="ImageMagick", version="0.7"),
			Pkg.PackageSpec(name="PlutoUI", version="0.7"),
			Pkg.PackageSpec(name="HypertextLiteral", version="0.5"),
			])
	Pkg.add("ImageIO")

	using Images
	using ImageIO
	using PlutoUI
	using HypertextLiteral
end

# ╔═╡ 29dfe3d6-c353-4081-8192-b12f374bf702
filter!(LOAD_PATH) do path
	path != "@v#.#"
end;

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

# **Homework 1** - _images and arrays_
`18.S191`, Spring 2021

`Due Date`: **Friday Feb 26, 2021 at 11:59pm EST**

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 = "Qling", kerberos_id = "rfhklwt")

# 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)@163.com)
"""

# ╔═╡ 5f95e01a-ee0a-11ea-030c-9dba276aba92
md"""
#### Intializing packages

_When running this notebook for the first time, this could take up to 15 minutes. Hang in there!_
"""

# ╔═╡ 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"""
$(html"<br>")
#### 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)

# ╔═╡ 5da8cbe8-eded-11ea-2e43-c5b7cc71e133
begin
	colored_line(x::Vector) = hcat(Gray.(Float64.(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)

# ╔═╡ 77adb065-bfd4-4680-9c2a-ad4d92689dbf
md"#### Exerise 1.2
👉 Make a function `my_sum` using a `for` loop, which computes the total of a vector of numbers."

# ╔═╡ bd907ee1-5253-4cae-b5a5-267dac24362a
function my_sum(xs)
	res = xs[1]

	for i in 2:length(xs)
		res += xs[i]
	end

	return res
end

# ╔═╡ 6640110a-d171-4b32-8d12-26979a36b718
my_sum([1,2,3])

# ╔═╡ cf738088-eded-11ea-2915-61735c2aa990
md"#### Exerise 1.3
👉 Use your `my_sum` function to write a function `mean`, which computes the mean/average of a vector of numbers."

# ╔═╡ 0ffa8354-edee-11ea-2883-9d5bfea4a236
function mean(xs)
	res = length(xs)

	return my_sum(xs) / res
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) # replace `missing` with your code!

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

# ╔═╡ ea8d92f8-159c-4161-8c54-bab7bc00f290
md"""
> ### Note about _mutation_
> There are two ways to think about this exercise, you could _modify_ the original vector, or you can _create a new vector_. We often prefer the second version, so that the original data is preserved. We generally only use code of the first variant in the most performance-sensitive parts of a program, as it requires more care to write and use correctly. _**Be careful not to get carried away in optimizing code**, especially when learning a new language!_
>
> There is a convention among Julians that functions that modify their argument have a `!` in their name. For example, `sort(x)` returns a sorted _copy_ of `x`, while `sort!(x)` _modifies_ `x` to be sorted.
>
> #### Tips for writing non-mutating code
> 1. _Rewriting_ an existing mutating function to be non-mutating can feel like a 'tedious' and 'inefficient' process. Often, instead of trying to **rewrite** a mutating function, it's best to take a step back and try to think of your problem as _constructing something new_. Instead of a `for` loop, it might make more sense to use **descriptive** primitives like [broadcasting with the dot syntax](https://docs.julialang.org/en/v1/manual/functions/#man-vectorized) (also for [math operators](https://docs.julialang.org/en/v1/manual/mathematical-operations/#man-dot-operators)), and [map and filter](https://www.youtube.com/watch?v=_O-HBDZMLrM).
>
>
> 2. If a mutating algorithm makes the most sense for your problem, then you can first use `copy` to create a copy of an array, and then modify that copy.
>
> We will cover this topic more in the later exercises!

"""

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

# ╔═╡ d6ddafdd-1a44-48c7-b49a-554073cdf331
test_vect = let

	# feel free to change your test case here!
	to_create = [-1.0, -1.5, 8.5]

	####
	# this cell is a bit funky to deal with a common pitfall from last year
	# it regenerates the vector if you accidentally wrote a mutating function

	# don't worry about how it works for this exercise!
	to_create
end

# ╔═╡ 29e10640-edf0-11ea-0398-17dbf4242de3
md"To verify our function, let's check that the mean of the `demean(test_vect)` is 0: (_Due to floating-point round-off error it may *not* be *exactly* 0._)"

# ╔═╡ 38155b5a-edf0-11ea-3e3f-7163da7433fb
demeaned_test_vect = demean(test_vect)

# ╔═╡ 1267e961-5b75-4b55-8080-d45316a03b9b
mean(demeaned_test_vect)

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

👉 Generate a vector of 100 elements. Where:
- the center 20 elements are set to `1`, and
- all other elements are `0`.
"""

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

# ╔═╡ 4a5e9d2c-dd90-4bb0-9e31-3f5c834406b4
create_bar()

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

# ╔═╡ 59414833-a108-4b1e-9a34-0f31dc907c6e
url = "https://user-images.githubusercontent.com/6933510/107239146-dcc3fd00-6a28-11eb-8c7b-41aaf6618935.png"

# ╔═╡ c5484572-ee05-11ea-0424-f37295c3072d
philip_filename = download(url) # download to a local file. The filename is returned

# ╔═╡ c8ecfe5c-ee05-11ea-322b-4b2714898831
philip = load(philip_filename)

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

# ╔═╡ 6ccd8902-0dd9-4335-a11a-ee7f9a1a6c09
philip_head = philip[470:800, 140:410]

# ╔═╡ 15088baa-c337-405d-8885-19a6e2bfd6aa
md"""
Recall from [Section 1.1](https://computationalthinking.mit.edu/Spring21/week1/) that in Julia, an _image_ is just a 2D array of color objects:
"""

# ╔═╡ 7ad4f6bc-588e-44ab-8ca9-c87cd1048442
typeof(philip)

# ╔═╡ a55bb5ca-600b-4aa0-b95f-7ece20845c9b
md"""
Every pixel (i.e. _element of the 2D array_) is of the `RGB` type:
"""

# ╔═╡ c5dc0cc8-9305-47e6-8b20-a9f8ef867799
philip_pixel = philip[100,100]

# ╔═╡ de772e21-0bea-4fd2-868a-9a7d32550bc9
typeof(philip_pixel)

# ╔═╡ 21bdc692-91ee-474d-ae98-455913a2342e
md"""
To get the values of its individual color channels, use the `r`, `g` and `b` _attributes_:
"""

# ╔═╡ 2ae3f379-96ce-435d-b863-deba4586ec71
philip_pixel.r, philip_pixel.g, philip_pixel.b

# ╔═╡ c49ba901-d798-489a-963c-4cc113c7abfd
md"""
And to create an `RGB` object yourself:
"""

# ╔═╡ 93451c37-77e1-4d4f-9788-c2a3da1401ee
RGB(0.1, 0.4, 0.7)

# ╔═╡ f52e4914-2926-4a42-9e45-9caaace9a7db
md"""
#### Exerise 2.1
👉 Write a function **`get_red`** that takes a single pixel, and returns the value of its red channel.
"""

# ╔═╡ a8b2270a-600c-4f83-939e-dc5ab35f4735
function get_red(pixel::AbstractRGB)
	return pixel.r
end

# ╔═╡ c320b39d-4cea-4fa1-b1ce-053c898a67a6
get_red(RGB(0.8, 0.1, 0.0))

# ╔═╡ d8cf9bd5-dbf7-4841-acf9-eef7e7cabab3
md"""
#### Exerise 2.2
👉 Write a function **`get_reds`** (note the extra `s`) that accepts a 2D color array called `image`, and returns a 2D array with the red channel value of each pixel. (The result should be a 2D array of _numbers_.) Use your function `get_red` from the previous exercise.
"""

# ╔═╡ ebe1d05c-f6aa-437d-83cb-df0ba30f20bf
function get_reds(image::AbstractMatrix)
	return get_red.(image)
end

# ╔═╡ c427554a-6f6a-43f1-b03b-f83239887cee
get_reds(philip_head)

# ╔═╡ 4fd07e01-2f8b-4ec9-9f4f-8a9e5ff56fb6
md"""

Great! By extracting the red channel value of each pixel, we get a 2D array of numbers. We went from an image (2D array of RGB colors) to a matrix (2D array of numbers).

#### Exerise 2.3
Let's try to visualize this matrix. Right now, it is displayed in text form, but because the image is quite large, most rows and columns don't fit on the screen. Instead, a better way to visualize it is to **view a number matrix as an image**.

This is easier than you might think! We just want to map each number to an `RGB` object, and the result will be a 2D array of `RGB` objects, which Julia will display as an image.

First, let's define a function that turns a _number_ into a _color_.
"""

# ╔═╡ 97c15896-6d99-4292-b7d7-4fcd2353656f
function value_as_color(x)

	return RGB(x, 0, 0)
end

# ╔═╡ cbb9bf41-4c21-42c7-b0e0-fc1ce29e0eb1
value_as_color(0.0), value_as_color(0.6), value_as_color(1.0)

# ╔═╡ 3f1a670b-44c2-4cab-909c-65f4ae9ed14b
md"""
👉 Use the functions `get_reds` and `value_as_color` to visualize the red channel values of `philip_head`. Tip: Like in previous exercises, use broadcasting ('dot syntax') to apply a function _element-wise_.

Use the ➕ button at the bottom left of this cell to add more cells.
"""

# ╔═╡ 21ba6e75-55a2-4614-9b5d-ea6378bf1d98
value_as_color.([0.0, 0.6, 1.0])

# ╔═╡ f7825c18-ff28-4e23-bf26-cc64f2f5049a
md"""

#### Exerise 2.4
👉 Write four more functions, `get_green`, `get_greens`, `get_blue` and `get_blues`, to be the equivalents of `get_red` and `get_reds`. Use the ➕ button at the bottom left of this cell to add new cells.
"""

# ╔═╡ 74b636a1-2505-4ce1-be45-b309954d0d59
get_green(pixel::AbstractRGB) = pixel.g

# ╔═╡ 0c773b57-55d5-4573-9c9d-214d4295a3c6
get_greens(image::AbstractMatrix) = get_green.(image)

# ╔═╡ e7768654-aeaa-40e8-afa3-57d8490082ac
get_blue(pixel::AbstractRGB) = pixel.b

# ╔═╡ f7706780-5ffc-405b-9980-dad7987f3670
get_blues(image::AbstractMatrix) = get_blue.(image)

# ╔═╡ c54ccdea-ee05-11ea-0365-23aaf053b7d7
md"""
#### Exercise 2.5
👉 Write a function **`mean_color`** that accepts an object called `image`. It should calculate the mean amounts of red, green and blue in the image and return the average color. Be sure to use functions from previous exercises!
"""

# ╔═╡ d9dfbbcc-dc02-4be3-9ca6-19f7df299f19
mean(get_blues(philip)) |> typeof

# ╔═╡ f6898df6-ee07-11ea-2838-fde9bc739c11
function mean_color(image)
	r = get_reds(image) |> mean
	g = get_greens(image) |> mean
	b = get_blues(image) |> mean
	return RGB(r, g, b)
end

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

# ╔═╡ 5f6635b4-63ed-4a62-969c-bd4084a8202f
md"""
_At the end of this homework, you can see all of your filters applied to your webcam image!_
"""

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

	return RGB(1 - color.r, 1 - color.g, 1 - color.b)
end

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

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

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

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

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

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

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

# ╔═╡ 55b138b7-19fb-4da1-9eb1-1e8304528251
md"""
_At the end of this homework, you can see all of your filters applied to your webcam image!_
"""

# ╔═╡ f68d4a36-ee07-11ea-0832-0360530f102e
md"""
#### Exercise 3.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.
"""

# ╔═╡ fbd1638d-8d7a-4d12-aff9-9c160cc3fd74
function quantize(x::Number)

	return floor(10x) / 10
end

# ╔═╡ 7720740e-2d2b-47f7-98fd-500ed3eee479
md"""
#### Intermezzo: _multiple methods_

In Julia, we often write multiple methods for the same function. When a function is called, the method is chosen based on the input arguments. Let's look at an example:

These are two _methods_ to the same function, because they have

> **the same name, but different input types**
"""

# ╔═╡ 90421bca-0804-4d6b-8bd0-3ddbd54cc5fe
function double(x::Number)

	return x * 2
end

# ╔═╡ b2329e4c-6204-453e-8998-2414b869b808
function double(x::Vector)

	return [x..., x...]
end

# ╔═╡ 23fcd65f-0182-41f3-80ec-d85b05136c47
md"""
When we call the function `double`, Julia will decide which method to call based on the given input argument!
"""

# ╔═╡ 5055b74c-b98d-41fa-a0d8-cb36200d82cc
double(24)

# ╔═╡ 8db17b2b-0cf9-40ba-8f6f-2e53be7b6355
double([1,2,37])

# ╔═╡ a8a597e0-a01c-40cd-9902-d56430afd938
md"""
We call this **multiple dispatch**, and it is one of Julia's key features. Throughout this course, you will see lots of real-world application, and learn to use multiple dispatch to create flexible and readable abstractions!
"""

# ╔═╡ f6b218c0-ee07-11ea-2adb-1968c4fd473a
md"""
#### Exercise 3.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`.

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. Use your previous method for `quantize`!
"""

# ╔═╡ 04e6b486-ceb7-45fe-a6ca-733703f16357
function quantize(color::AbstractRGB)

	return RGB(quantize(color.r), quantize(color.g), quantize(color.b))
end

# ╔═╡ f6bf64da-ee07-11ea-3efb-05af01b14f67
md"""
#### Exercise 3.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.)
"""

# ╔═╡ 13e9ec8d-f615-4833-b1cf-0153010ccb65
quantize(image::AbstractMatrix) = quantize.(image)

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

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

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

# ╔═╡ f6d6c71a-ee07-11ea-2b63-d759af80707b
md"""
#### Exercise 3.5
👉 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. (Julia has a built-in `clamp` function, or you can write your own function.)
"""

# ╔═╡ f38b198d-39cf-456f-a841-1ba08f206010
function noisify(x::Number, s)
	seed = rand() * 2s - s

	return clamp(x + seed, 0, 1)
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.

Use your previous method for `noisify`. _(Remember that Julia chooses which method to use based on the input arguments. So to call the method from the previous exercise, the first argument should be a number.)_
"""

# ╔═╡ db4bad9f-df1c-4640-bb34-dd2fe9bdce18
function noisify(color::AbstractRGB, s)

	return RGB(noisify(color.r, s), noisify(color.g, s), noisify(color.b, s))
end

# ╔═╡ 0000b7f8-4c43-4dd8-8665-0dfe59e74c0a
md"""
Noise strength:
"""

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

# ╔═╡ 48de5bc2-72d3-11eb-3fd9-eff2b686cb75
md"""
> ### Note about _array comprehension_
> At this point, you already know of a few ways to make a new array based on one that already exists.
> 1. you can use a for loop to go through a array
> 1. you can use function broadcasting over a array
> 1. you can use _**array comprehension**_!
>
> The third option you are about to see demonstrated below and following the following syntax:
>
> ```[function_to_apply(args) for args in some_iterable_of_your_choice]```
>
> This creates a new iterable that matches what you iterate through in the second part of the comprehension. Below is an example with `for` loops through two iterables that creates a 2-dimensional `Array`.
"""

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

# ╔═╡ 21a5885d-00ab-428b-96c3-c28c98c4ca6d
function noisify(image::AbstractMatrix, s)

	return noisify.(image, s)
end

# ╔═╡ 1ea53f41-b791-40e2-a0f8-04e13d856829
noisify(0.5, 0.1) # edit this test case!

# ╔═╡ 7e4aeb70-ee1b-11ea-100f-1952ba66f80f
(original=color_red, with_noise=noisify(color_red, color_noise))

# ╔═╡ 8e848279-1b3e-4f32-8c0c-45693d12de96
[
	noisify(color_red, strength)
	for
		strength in 0 : 0.05 : 1,
		row in 1:10
]'

# ╔═╡ d896b7fd-20db-4aa9-bbcf-81b1cd44ec46
md"""

#### Exerise 3.6
Move the slider below to set the amount of noise applied to the image of Philip.
"""

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

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

# ╔═╡ 9604bc44-ee1b-11ea-28f8-7f7af8d0cbb2
if philip_noise == 1
	md"""
	> #### What's this?
	>
	> The noise intensity is `1.0`, but we can still recognise Philip in the picture...
	>
	> 👉 Modify the definition of the slider to go further than `1.0`.
	"""
end

# ╔═╡ 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 1
"""

# ╔═╡ e87e0d14-43a5-490d-84d9-b14ece472061
md"""
### Results
"""

# ╔═╡ ee5f21fb-1076-42b6-8926-8bbb6ed0ad67
function custom_filter(pixel::AbstractRGB)

	# your code here!

	return pixel
end

# ╔═╡ 9e5a08dd-332a-486b-94ab-15c49e72e522
function custom_filter(image::AbstractMatrix)

	return custom_filter.(image)
end

# ╔═╡ 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

# ╔═╡ 756d150a-b7bf-4bf5-b372-5b0efa80d987
md"## Function library

Just some helper functions used in the notebook."

# ╔═╡ 4bc94bec-da39-4f8a-82ee-9953ed73b6a4
hint(text) = Markdown.MD(Markdown.Admonition("hint", "Hint", [text]))

# ╔═╡ b1d5ca28-edf6-11ea-269e-75a9fb549f1d
md"""
You can find out more about any function (like `rand`) by clicking on the Live Docs in the bottom right of this Pluto window, and typing a function name in the top.

![image](https://user-images.githubusercontent.com/6933510/107848812-c934df80-6df6-11eb-8a32-663d802f5d11.png)

![image](https://user-images.githubusercontent.com/6933510/107848846-0f8a3e80-6df7-11eb-818a-7271ecb9e127.png)

We recommend that you leave the window open while you work on Julia code. It will continually look up documentation for anything you type!

#### Help, I don't see the Live Docs!

Try the following:

🙋 **Are you viewing a static preview?** The Live Docs only work if you _run_ the notebook. If you are reading this on our course website, then click the button in the top right to run the notebook.

🙋 **Is your screen too small?** Try resizing your window or zooming out.
""" |> hint

# ╔═╡ 24090306-7395-4f2f-af31-34f7486f3945
hint(md"""Check out this page for a refresher on basic Julia syntax:

	[Basic Julia Syntax](https://computationalthinking.mit.edu/Spring21/basic_syntax/)""")

# ╔═╡ aa1ff74a-4e78-4ef1-8b8d-3a60a168cf6d
hint(md"""
In [Section 1.1](https://computationalthinking.mit.edu/Spring21/week1/), we drew a red square on top of the image Philip with a simple command...
""")

# ╔═╡ 50e2b0fb-b06d-4ac1-bdfb-eab833466736
md"""
This exercise can be quite difficult if you use a `for` loop or list comprehension.

Instead, you should use the [dot syntax](https://docs.julialang.org/en/v1/manual/functions/#man-vectorized) to apply a function _element-wise_ to an array. For example, this is how you get the square root of `3`:

```
sqrt(3)
```

and this is how you get the square roots of 1, 2 and 3:

```
sqrt.([1, 2, 3])
```

""" |> hint

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

# ╔═╡ 8ce6ad06-819c-4af5-bed7-56ecc08c97be
almost(text) = Markdown.MD(Markdown.Admonition("warning", "Almost there!", [text]))

# ╔═╡ dfa40e89-03fc-4a7a-825e-92d67ee217b2
still_missing(text=md"Replace `missing` with your answer.") = Markdown.MD(Markdown.Admonition("warning", "Here we go!", [text]))

# ╔═╡ 086ec1ff-b62d-4566-9973-5b2cc3353409
keep_working(text=md"The answer is not quite right.") = Markdown.MD(Markdown.Admonition("danger", "Keep working on it!", [text]))

# ╔═╡ 2f6fb3a6-bb5d-4c7a-9297-84dd4b16c7c3
yays = [md"Fantastic!", md"Splendid!", 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."]

# ╔═╡ c22f688b-dc04-4a94-b541-fe06266c5446
correct(text=rand(yays)) = Markdown.MD(Markdown.Admonition("correct", "Got it!", [text]))

# ╔═╡ 09102183-f9fb-4d89-b4f9-5d76af7b8e90
let
	result = get_red(RGB(0.2, 0.3, 0.4))
	if ismissing(result)
		still_missing()
	elseif isnothing(result)
		keep_working(md"Did you forget to write `return`?")
	elseif result == 0.2
		correct()
	else
		keep_working()
	end
end

# ╔═╡ 63ac142e-6d9d-4109-9286-030a02c900b4
let
	test = [RGB(0.2, 0, 0)   RGB(0.6, 0, 0)]
	result = get_reds(test)

	if ismissing(result)
		still_missing()
	elseif isnothing(result)
		keep_working(md"Did you forget to write `return`?")
	elseif result == [ 0.2  0.6 ]
		correct()
	else
		keep_working()
	end
end

# ╔═╡ 80a4cb23-49c9-4446-a3ec-b2203128dc27
let
	result = invert(RGB(1.0, 0.5, 0.25)) # I chose these values because they can be represented exactly by Float64
	shouldbe = RGB(0.0, 0.5, 0.75)

	if ismissing(result)
		still_missing()
	elseif isnothing(result)
		keep_working(md"Did you forget to write `return`?")
	elseif !(result isa AbstractRGB)
		keep_working(md"You need to return a _color_, i.e. an object of type `RGB`. Use `RGB(r, g, b)` to create a color with channel values `r`, `g` and `b`.")
	elseif !(result == shouldbe)
		keep_working()
	else
		correct()
	end
end

# ╔═╡ a6d9635b-85ed-4590-ad09-ca2903ea8f1d
let
	result = quantize(RGB(.297, .1, .0))

	if ismissing(result)
		still_missing()
	elseif isnothing(result)
		keep_working(md"Did you forget to write `return`?")
	elseif !(result isa AbstractRGB)
		keep_working(md"You need to return a _color_, i.e. an object of type `RGB`. Use `RGB(r, g, b)` to create a color with channel values `r`, `g` and `b`.")
	elseif result != RGB(0.2, .1, .0)
		keep_working()
	else
		correct()
	end
end

# ╔═╡ 31ef3710-e4c9-4aa7-bd8f-c69cc9a977ee
let
	result = noisify(0.5, 0)

	if ismissing(result)
		still_missing()
	elseif isnothing(result)
		keep_working(md"Did you forget to write `return`?")
	elseif result == 0.5

		results = [noisify(0.9, 0.1) for _ in 1:1000]

		if 0.8 ≤ minimum(results) < 0.81 && 0.99 ≤ maximum(results) ≤ 1
			result = noisify(5, 3)

			if result == 1
				correct()
			else
				keep_working(md"The result should be restricted to the range ``[0,1]``.")
			end
		else
			keep_working()
		end
	else
		keep_working(md"What should `noisify(0.5, 0)` be?")
		correct()
	end
end

# ╔═╡ ab3d1b70-88e8-4118-8d3e-601a8a68f72d
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 eltype(random_vect) != Float64
	almost(md"""
		You generated a vector of random integers. For the remaining exercises, we want a vector of `Float64` numbers.

		The (optional) first argument to `rand` specifies the **type** of elements to generate. For example: `rand(Bool, 10)` generates 10 values that are either `true` or `false`. (Try it!)
		""")
elseif length(random_vect) != 10
	keep_working(md"`random_vect` does not have the correct size.")
elseif length(Set(random_vect)) != 10
	keep_working(md"`random_vect` is not 'random enough'")
else
	correct(md"Well done! You can run your code again to generate a new vector!")
end

# ╔═╡ e0bfc973-2808-4f84-b065-fb3d05401e30
if !@isdefined(my_sum)
	not_defined(:my_sum)
else
	let
		result = my_sum([1,2,3])
		if ismissing(result)
			still_missing()
		elseif isnothing(result)
			keep_working(md"Did you forget to write `return`?")
		elseif result != 6
			keep_working()
		else
			correct()
		end
	end
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

# ╔═╡ adf476d8-a334-4b35-81e8-cc3b37de1f28
if !@isdefined(mean)
	not_defined(:mean)
else
	let
		input = Float64[1,2,3]
		result = demean(input)

		if input === result
			almost(md"""
			It looks like you **modified** `xs` inside the function.

			It is preferable to avoid mutation inside functions, because you might want to use the original data again. For example, applying `demean` to a dataset of sensor readings would **modify** the original data, and the rest of your analysis would be erroneous.

			""")
		elseif ismissing(result)
			still_missing()
		elseif isnothing(result)
			keep_working(md"Did you forget to write `return`?")
		elseif !(result isa AbstractVector) || length(result) != 3
			keep_working(md"Return a vector of the same size as `xs`.")
		elseif abs(sum(result) / 3) < 1e-10
			correct()
		else
			keep_working()
		end
	end
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

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

		result = mean_color(input)
		shouldbe = RGB(1.0, 1.0, 0.5)

		if ismissing(result)
			still_missing()
		elseif isnothing(result)
			keep_working(md"Did you forget to write `return`?")
		elseif !(result isa AbstractRGB)
			keep_working(md"You need to return a _color_, i.e. an object of type `RGB`. Use `RGB(r, g, b)` to create a color with channel values `r`, `g` and `b`.")
		elseif !(result == shouldbe)
			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

# ╔═╡ 8cb0aee8-5774-4490-9b9e-ada93416c089
todo(text) = HTML("""<div
	style="background: rgb(220, 200, 255); padding: 2em; border-radius: 1em;"
	><h1>TODO</h1>$(repr(MIME"text/html"(), text))</div>""")

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

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

# ╔═╡ e083b3e8-ed61-11ea-2ec9-217820b0a1b4
md"""
 $(bigbreak)

## **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.
"""

# ╔═╡ f6cc03a0-ee07-11ea-17d8-013991514d42
md"""
 $(bigbreak)

## Exercise 3 - _More filters_

In the previous exercises, we learned how to use Julia's _dot syntax_ to apply a function _element-wise_ to an array. In this exercise, we will use this to write more image filters, that you can then apply to your own webcam image!

#### Exercise 3.1
👉 Write a function `invert` that inverts a color, i.e. sends $(r, g, b)$ to $(1 - r, 1-g, 1-b)$.
"""

# ╔═╡ 4139ee66-ee0a-11ea-2282-15d63bcca8b8
md"""
$(bigbreak)
### Camera input
"""

# ╔═╡ 87dabfd2-461e-4769-ad0f-132cb2370b88
md"""
$(bigbreak)
### Write your own filter!
"""

# ╔═╡ 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

# ╔═╡ 20402780-426b-4caa-af8f-ff1e7787b7f9
@bind cam_data camera_input()

# ╔═╡ 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

# ╔═╡ ed9fb2ac-2680-42b7-9b00-591e45a5e105
cam_image = process_raw_camera_data(cam_data)

# ╔═╡ d38c6958-9300-4f7a-89cf-95ca9e899c13
mean_color(cam_image)

# ╔═╡ 82f1e006-60fe-4ad1-b9cb-180fafdeb4da
invert.(cam_image)

# ╔═╡ 54c83589-b8c6-422a-b5e9-d8e0ee72a224
quantize(cam_image)

# ╔═╡ 18e781f8-66f3-4216-bc84-076a08f9f3fb
noisify(cam_image, .5)

# ╔═╡ ebf3193d-8c8d-4425-b252-45067a5851d9
[
	invert.(cam_image)      quantize(cam_image)
	noisify(cam_image, .5)  custom_filter(cam_image)
]

# ╔═╡ 8917529e-fa7a-412b-8aea-54f92f6270fa
custom_filter(cam_image)

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

# ╔═╡ Cell order:
# ╟─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
# ╟─29dfe3d6-c353-4081-8192-b12f374bf702
# ╟─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
# ╠═5da8cbe8-eded-11ea-2e43-c5b7cc71e133
# ╟─77adb065-bfd4-4680-9c2a-ad4d92689dbf
# ╠═bd907ee1-5253-4cae-b5a5-267dac24362a
# ╠═6640110a-d171-4b32-8d12-26979a36b718
# ╟─e0bfc973-2808-4f84-b065-fb3d05401e30
# ╟─24090306-7395-4f2f-af31-34f7486f3945
# ╟─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
# ╟─ea8d92f8-159c-4161-8c54-bab7bc00f290
# ╠═ec5efe8c-edef-11ea-2c6f-afaaeb5bc50c
# ╟─d6ddafdd-1a44-48c7-b49a-554073cdf331
# ╟─29e10640-edf0-11ea-0398-17dbf4242de3
# ╠═1267e961-5b75-4b55-8080-d45316a03b9b
# ╠═38155b5a-edf0-11ea-3e3f-7163da7433fb
# ╟─adf476d8-a334-4b35-81e8-cc3b37de1f28
# ╟─a5f8bafe-edf0-11ea-0da3-3330861ae43a
# ╠═b6b65b94-edf0-11ea-3686-fbff0ff53d08
# ╠═4a5e9d2c-dd90-4bb0-9e31-3f5c834406b4
# ╟─d862fb16-edf1-11ea-36ec-615d521e6bc0
# ╟─aa1ff74a-4e78-4ef1-8b8d-3a60a168cf6d
# ╟─e3394c8a-edf0-11ea-1bb8-619f7abb6881
# ╟─e083b3e8-ed61-11ea-2ec9-217820b0a1b4
# ╠═59414833-a108-4b1e-9a34-0f31dc907c6e
# ╠═c5484572-ee05-11ea-0424-f37295c3072d
# ╠═c8ecfe5c-ee05-11ea-322b-4b2714898831
# ╟─e86ed944-ee05-11ea-3e0f-d70fc73b789c
# ╠═6ccd8902-0dd9-4335-a11a-ee7f9a1a6c09
# ╟─15088baa-c337-405d-8885-19a6e2bfd6aa
# ╠═7ad4f6bc-588e-44ab-8ca9-c87cd1048442
# ╟─a55bb5ca-600b-4aa0-b95f-7ece20845c9b
# ╠═c5dc0cc8-9305-47e6-8b20-a9f8ef867799
# ╠═de772e21-0bea-4fd2-868a-9a7d32550bc9
# ╟─21bdc692-91ee-474d-ae98-455913a2342e
# ╠═2ae3f379-96ce-435d-b863-deba4586ec71
# ╟─c49ba901-d798-489a-963c-4cc113c7abfd
# ╠═93451c37-77e1-4d4f-9788-c2a3da1401ee
# ╟─f52e4914-2926-4a42-9e45-9caaace9a7db
# ╠═a8b2270a-600c-4f83-939e-dc5ab35f4735
# ╠═c320b39d-4cea-4fa1-b1ce-053c898a67a6
# ╟─09102183-f9fb-4d89-b4f9-5d76af7b8e90
# ╟─d8cf9bd5-dbf7-4841-acf9-eef7e7cabab3
# ╠═ebe1d05c-f6aa-437d-83cb-df0ba30f20bf
# ╠═c427554a-6f6a-43f1-b03b-f83239887cee
# ╟─63ac142e-6d9d-4109-9286-030a02c900b4
# ╟─50e2b0fb-b06d-4ac1-bdfb-eab833466736
# ╟─4fd07e01-2f8b-4ec9-9f4f-8a9e5ff56fb6
# ╠═97c15896-6d99-4292-b7d7-4fcd2353656f
# ╠═cbb9bf41-4c21-42c7-b0e0-fc1ce29e0eb1
# ╟─3f1a670b-44c2-4cab-909c-65f4ae9ed14b
# ╠═21ba6e75-55a2-4614-9b5d-ea6378bf1d98
# ╟─f7825c18-ff28-4e23-bf26-cc64f2f5049a
# ╠═74b636a1-2505-4ce1-be45-b309954d0d59
# ╠═0c773b57-55d5-4573-9c9d-214d4295a3c6
# ╠═e7768654-aeaa-40e8-afa3-57d8490082ac
# ╠═f7706780-5ffc-405b-9980-dad7987f3670
# ╟─c54ccdea-ee05-11ea-0365-23aaf053b7d7
# ╠═d9dfbbcc-dc02-4be3-9ca6-19f7df299f19
# ╠═f6898df6-ee07-11ea-2838-fde9bc739c11
# ╠═5be9b144-ee0d-11ea-2a8d-8775de265a1d
# ╟─4d0158d0-ee0d-11ea-17c3-c169d4284acb
# ╟─5f6635b4-63ed-4a62-969c-bd4084a8202f
# ╟─f6cc03a0-ee07-11ea-17d8-013991514d42
# ╠═63e8d636-ee0b-11ea-173d-bd3327347d55
# ╟─80a4cb23-49c9-4446-a3ec-b2203128dc27
# ╟─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
# ╟─55b138b7-19fb-4da1-9eb1-1e8304528251
# ╟─f68d4a36-ee07-11ea-0832-0360530f102e
# ╠═fbd1638d-8d7a-4d12-aff9-9c160cc3fd74
# ╠═f6a655f8-ee07-11ea-13b6-43ca404ddfc7
# ╟─c905b73e-ee1a-11ea-2e36-23b8e73bfdb6
# ╟─7720740e-2d2b-47f7-98fd-500ed3eee479
# ╠═90421bca-0804-4d6b-8bd0-3ddbd54cc5fe
# ╠═b2329e4c-6204-453e-8998-2414b869b808
# ╟─23fcd65f-0182-41f3-80ec-d85b05136c47
# ╠═5055b74c-b98d-41fa-a0d8-cb36200d82cc
# ╠═8db17b2b-0cf9-40ba-8f6f-2e53be7b6355
# ╟─a8a597e0-a01c-40cd-9902-d56430afd938
# ╟─f6b218c0-ee07-11ea-2adb-1968c4fd473a
# ╠═04e6b486-ceb7-45fe-a6ca-733703f16357
# ╟─a6d9635b-85ed-4590-ad09-ca2903ea8f1d
# ╟─f6bf64da-ee07-11ea-3efb-05af01b14f67
# ╠═13e9ec8d-f615-4833-b1cf-0153010ccb65
# ╟─25dad7ce-ee0b-11ea-3e20-5f3019dd7fa3
# ╠═9751586e-ee0c-11ea-0cbb-b7eda92977c9
# ╟─f6d6c71a-ee07-11ea-2b63-d759af80707b
# ╠═f38b198d-39cf-456f-a841-1ba08f206010
# ╠═1ea53f41-b791-40e2-a0f8-04e13d856829
# ╟─31ef3710-e4c9-4aa7-bd8f-c69cc9a977ee
# ╟─f6ef2c2e-ee07-11ea-13a8-2512e7d94426
# ╟─f6fc1312-ee07-11ea-39a0-299b67aee3d8
# ╠═db4bad9f-df1c-4640-bb34-dd2fe9bdce18
# ╟─0000b7f8-4c43-4dd8-8665-0dfe59e74c0a
# ╠═774b4ce6-ee1b-11ea-2b48-e38ee25fc89b
# ╠═7e4aeb70-ee1b-11ea-100f-1952ba66f80f
# ╟─48de5bc2-72d3-11eb-3fd9-eff2b686cb75
# ╠═8e848279-1b3e-4f32-8c0c-45693d12de96
# ╟─f70823d2-ee07-11ea-2bb3-01425212aaf9
# ╠═21a5885d-00ab-428b-96c3-c28c98c4ca6d
# ╟─d896b7fd-20db-4aa9-bbcf-81b1cd44ec46
# ╠═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
# ╟─4139ee66-ee0a-11ea-2282-15d63bcca8b8
# ╠═20402780-426b-4caa-af8f-ff1e7787b7f9
# ╟─ed9fb2ac-2680-42b7-9b00-591e45a5e105
# ╟─e87e0d14-43a5-490d-84d9-b14ece472061
# ╠═d38c6958-9300-4f7a-89cf-95ca9e899c13
# ╠═82f1e006-60fe-4ad1-b9cb-180fafdeb4da
# ╠═54c83589-b8c6-422a-b5e9-d8e0ee72a224
# ╠═18e781f8-66f3-4216-bc84-076a08f9f3fb
# ╠═ebf3193d-8c8d-4425-b252-45067a5851d9
# ╟─87dabfd2-461e-4769-ad0f-132cb2370b88
# ╠═8917529e-fa7a-412b-8aea-54f92f6270fa
# ╠═ee5f21fb-1076-42b6-8926-8bbb6ed0ad67
# ╠═9e5a08dd-332a-486b-94ab-15c49e72e522
# ╟─91f4778e-ee20-11ea-1b7e-2b0892bd3c0f
# ╟─8ffe16ce-ee20-11ea-18bd-15640f94b839
# ╟─5842895a-ee10-11ea-119d-81e4c4c8c53b
# ╟─756d150a-b7bf-4bf5-b372-5b0efa80d987
# ╟─4bc94bec-da39-4f8a-82ee-9953ed73b6a4
# ╟─8ce6ad06-819c-4af5-bed7-56ecc08c97be
# ╟─dfa40e89-03fc-4a7a-825e-92d67ee217b2
# ╟─086ec1ff-b62d-4566-9973-5b2cc3353409
# ╟─2f6fb3a6-bb5d-4c7a-9297-84dd4b16c7c3
# ╟─c22f688b-dc04-4a94-b541-fe06266c5446
# ╟─ab3d1b70-88e8-4118-8d3e-601a8a68f72d
# ╟─8cb0aee8-5774-4490-9b9e-ada93416c089
# ╟─115ded8c-ee0a-11ea-3493-89487315feb7
# ╟─dfb7c6be-ee0d-11ea-194e-9758857f7b20
# ╟─e15ad330-ee0d-11ea-25b6-1b1b3f3d7888
# ╟─83eb9ca0-ed68-11ea-0bc5-99a09c68f867


text:
#
# Graphs have a collapsable header and plotting area
# The header has fields for:
#   horizontal range, vertical range
#   cursor position,
#
# A graph renders data from a data collection.

module Graphs

using CImGui
using CImGui.CSyntax
using CImGui.CSyntax.CStatic
using CImGui.GLFWBackend
using CImGui.OpenGLBackend
using CImGui.GLFWBackend.GLFW
using CImGui.OpenGLBackend.ModernGL
using CImGui: ImVec2, ImVec4, IM_COL32, ImU32, ImColor
using Printf

using ..Gui

const GRAPH_WINDOW_WIDTH = Gui.GUI_WIDTH - 10
const GRAPH_WINDOW_HEIGHT = 200

const MAX_VERTICAL_BARS = 250
const LINE_THICKNESS = 1.0

const GRAY = 64
const YELLOW = IM_COL32(255, 255, 0, 255)
const GREEN = IM_COL32(0, 255, 0, 255)
const GREEN_TRAN = IM_COL32(0, 255, 0, 128)
const WHITE_TRAN = IM_COL32(255, 255, 255, 100)
const WHITE = IM_COL32(255, 255, 255, 255)
const LIME_GREEN = IM_COL32(166, 255, 77, 255)
const BLUE = IM_COL32(0, 0, 255, 255)
const LIGHT_BLUE = IM_COL32(0, 169, 224, 255) # Pantone light blue
const ORANGE = IM_COL32(255, 128, 0, 255)
const GREY = IM_COL32(100, 100, 100, 255)
const LIGHT_GREY = IM_COL32(200, 200, 200, 255)

abstract type AbstractGraph end

export linear, lerp

# linear returns 0->1 for a "value" between min and max.
# Generally used to map from view/data-space to unit-space
function linear(min::Float64, max::Float64, value::Float64)
   	if min < 0
        return 1.0 - ((value - max) / (min - max))
    end

   	(value - min) / (max - min)
end

# lerp returns a the value between min and max given t = 0->1
# Call linear() first to get "t".
# Generally used to map from unit-space to window-space.
# A graph will then map to local-space of the graph.
function lerp(min::Float64, max::Float64, t::Float64)
    min * (1.0 - t) + max * t
end

# Map from sample-space to unit-space where unit-space is 0->1
function map_sample_to_unit(v::Float64, min::Float64, max::Float64)
    linear(min, max, v)
end

# Map from unit-space to window-space
function map_unit_to_window(v::Float64, min::Float64, max::Float64)
    lerp(min, max, v)
end

# Local = graph-space
function map_window_to_local(x::Float64, y::Float64, offsets::CImGui.LibCImGui.ImVec2)
    (offsets.x + x, offsets.y + y)
end

function scroll_velocity(scroll::Float64)
    sign(scroll) * exp(sign(scroll) * scroll)
end

include("spike_scatter_graph.jl")
include("soma_ap_fast_graph.jl")
include("soma_ap_slow_graph.jl")
include("soma_psp_graph.jl")
include("synapse_weights_graph.jl")
include("synapse_surge_graph.jl")
include("synapse_psp_graph.jl")
# include("synapse_input_graph.jl")

spikes_graph = Graphs.SpikeScatterGraph()
soma_apFast_graph = Graphs.SomaAPFastGraph()
soma_apSlow_graph = Graphs.SomaAPSlowGraph()
soma_psp_graph = Graphs.SomaPSPGraph()
synapse_weights_graph = Graphs.SynapseWeightGraph()
synapse_surge_graph = Graphs.SynapseSurgeGraph()
synapse_psp_graph = Graphs.SynapsePSPGraph()
# synapse_input_graph = Graphs.SynapticInputGraph()

end # Module ---------------------------------------

ext: jl


text:
# everything related to a DG semidiscretization in 2D using MPI,
# currently limited to Lobatto-Legendre nodes

# TODO: MPI dimension agnostic
# TODO: MPI, adapt to different real types (and AD!)
mutable struct MPICache
  mpi_neighbor_ranks::Vector{Int}
  mpi_neighbor_interfaces::Vector{Vector{Int}}
  mpi_send_buffers::Vector{Vector{Float64}}
  mpi_recv_buffers::Vector{Vector{Float64}}
  mpi_send_requests::Vector{MPI.Request}
  mpi_recv_requests::Vector{MPI.Request}
  n_elements_by_rank::OffsetArray{Int, 1, Array{Int, 1}}
  n_elements_global::Int
  first_element_global_id::Int
end

function MPICache()
  mpi_neighbor_ranks = Vector{Int}(undef, 0)
  mpi_neighbor_interfaces = Vector{Vector{Int}}(undef, 0)
  mpi_send_buffers = Vector{Vector{Float64}}(undef, 0)
  mpi_recv_buffers = Vector{Vector{Float64}}(undef, 0)
  mpi_send_requests = Vector{MPI.Request}(undef, 0)
  mpi_recv_requests = Vector{MPI.Request}(undef, 0)
  n_elements_by_rank = OffsetArray(Vector{Int}(undef, 0), 0:-1)
  n_elements_global = 0
  first_element_global_id = 0

  MPICache(mpi_neighbor_ranks, mpi_neighbor_interfaces,
           mpi_send_buffers, mpi_recv_buffers,
           mpi_send_requests, mpi_recv_requests,
           n_elements_by_rank, n_elements_global,
           first_element_global_id)
end

# TODO: MPI dimension agnostic
function start_mpi_receive!(mpi_cache::MPICache)

  for (index, d) in enumerate(mpi_cache.mpi_neighbor_ranks)
    mpi_cache.mpi_recv_requests[index] = MPI.Irecv!(
      mpi_cache.mpi_recv_buffers[index], d, d, mpi_comm())
  end

  return nothing
end

# TODO: MPI dimension agnostic
function start_mpi_send!(mpi_cache::MPICache, mesh, equations, dg, cache)
  data_size = nvariables(equations) * nnodes(dg)^(ndims(mesh) - 1)

  for d in 1:length(mpi_cache.mpi_neighbor_ranks)
    send_buffer = mpi_cache.mpi_send_buffers[d]

    for (index, interface) in enumerate(mpi_cache.mpi_neighbor_interfaces[d])
      first = (index - 1) * data_size + 1
      last =  (index - 1) * data_size + data_size

      if cache.mpi_interfaces.remote_sides[interface] == 1 # local element in positive direction
        @views send_buffer[first:last] .= vec(cache.mpi_interfaces.u[2, :, :, interface])
      else # local element in negative direction
        @views send_buffer[first:last] .= vec(cache.mpi_interfaces.u[1, :, :, interface])
      end
    end
  end

  # Start sending
  for (index, d) in enumerate(mpi_cache.mpi_neighbor_ranks)
    mpi_cache.mpi_send_requests[index] = MPI.Isend(
      mpi_cache.mpi_send_buffers[index], d, mpi_rank(), mpi_comm())
  end

  return nothing
end

# TODO: MPI dimension agnostic
function finish_mpi_send!(mpi_cache::MPICache)
  MPI.Waitall!(mpi_cache.mpi_send_requests)
end

# TODO: MPI dimension agnostic
function finish_mpi_receive!(mpi_cache::MPICache, mesh, equations, dg, cache)
  data_size = nvariables(equations) * nnodes(dg)^(ndims(mesh) - 1)

  # Start receiving and unpack received data until all communication is finished
  d, _ = MPI.Waitany!(mpi_cache.mpi_recv_requests)
  while d != 0
    recv_buffer = mpi_cache.mpi_recv_buffers[d]

    for (index, interface) in enumerate(mpi_cache.mpi_neighbor_interfaces[d])
      first = (index - 1) * data_size + 1
      last =  (index - 1) * data_size + data_size

      if cache.mpi_interfaces.remote_sides[interface] == 1 # local element in positive direction
        @views vec(cache.mpi_interfaces.u[1, :, :, interface]) .= recv_buffer[first:last]
      else # local element in negative direction
        @views vec(cache.mpi_interfaces.u[2, :, :, interface]) .= recv_buffer[first:last]
      end
    end

    d, _ = MPI.Waitany!(mpi_cache.mpi_recv_requests)
  end

  return nothing
end

# This method is called when a SemidiscretizationHyperbolic is constructed.
# It constructs the basic `cache` used throughout the simulation to compute
# the RHS etc.
function create_cache(mesh::ParallelTreeMesh{2}, equations,
                      dg::DG, RealT, uEltype)
  # Get cells for which an element needs to be created (i.e. all leaf cells)
  leaf_cell_ids = local_leaf_cells(mesh.tree)

  elements = init_elements(leaf_cell_ids, mesh, equations, dg.basis, RealT, uEltype)

  interfaces = init_interfaces(leaf_cell_ids, mesh, elements)

  mpi_interfaces = init_mpi_interfaces(leaf_cell_ids, mesh, elements)

  boundaries = init_boundaries(leaf_cell_ids, mesh, elements)

  mortars = init_mortars(leaf_cell_ids, mesh, elements, dg.mortar)

  mpi_cache = init_mpi_cache(mesh, elements, mpi_interfaces, nvariables(equations), nnodes(dg))

  cache = (; elements, interfaces, mpi_interfaces, boundaries, mortars,
             mpi_cache)

  # Add specialized parts of the cache required to compute the volume integral etc.
  cache = (;cache..., create_cache(mesh, equations, dg.volume_integral, dg, uEltype)...)
  cache = (;cache..., create_cache(mesh, equations, dg.mortar, uEltype)...)

  return cache
end

function init_mpi_cache(mesh, elements, mpi_interfaces, nvars, nnodes)
  mpi_cache = MPICache()

  init_mpi_cache!(mpi_cache, mesh, elements, mpi_interfaces, nvars, nnodes)
  return mpi_cache
end

function init_mpi_cache!(mpi_cache, mesh, elements, mpi_interfaces, nvars, nnodes)
  mpi_neighbor_ranks, mpi_neighbor_interfaces =
    init_mpi_neighbor_connectivity(elements, mpi_interfaces, mesh)

  mpi_send_buffers, mpi_recv_buffers, mpi_send_requests, mpi_recv_requests =
    init_mpi_data_structures(mpi_neighbor_interfaces, ndims(mesh), nvars, nnodes)

  # Determine local and total number of elements
  n_elements_by_rank = Vector{Int}(undef, mpi_nranks())
  n_elements_by_rank[mpi_rank() + 1] = nelements(elements)
  MPI.Allgather!(MPI.UBuffer(n_elements_by_rank, 1), mpi_comm())
  n_elements_by_rank = OffsetArray(n_elements_by_rank, 0:(mpi_nranks() - 1))
  n_elements_global = MPI.Allreduce(nelements(elements), +, mpi_comm())
  @assert n_elements_global == sum(n_elements_by_rank) "error in total number of elements"

  # Determine the global element id of the first element
  first_element_global_id = MPI.Exscan(nelements(elements), +, mpi_comm())
  if mpi_isroot()
    # With Exscan, the result on the first rank is undefined
    first_element_global_id = 1
  else
    # On all other ranks we need to add one, since Julia has one-based indices
    first_element_global_id += 1
  end
  # TODO reuse existing structures
  @pack! mpi_cache = mpi_neighbor_ranks, mpi_neighbor_interfaces,
                     mpi_send_buffers, mpi_recv_buffers,
                     mpi_send_requests, mpi_recv_requests,
                     n_elements_by_rank, n_elements_global,
                     first_element_global_id
end

# Initialize connectivity between MPI neighbor ranks
function init_mpi_neighbor_connectivity(elements, mpi_interfaces, mesh::TreeMesh2D)
  tree = mesh.tree

  # Determine neighbor ranks and sides for MPI interfaces
  neighbor_ranks = fill(-1, nmpiinterfaces(mpi_interfaces))
  # The global interface id is the smaller of the (globally unique) neighbor cell ids, multiplied by
  # number of directions (2 * ndims) plus direction minus one
  global_interface_ids = fill(-1, nmpiinterfaces(mpi_interfaces))
  for interface_id in 1:nmpiinterfaces(mpi_interfaces)
    orientation = mpi_interfaces.orientations[interface_id]
    remote_side = mpi_interfaces.remote_sides[interface_id]
    # Direction is from local cell to remote cell
    if orientation == 1 # MPI interface in x-direction
      if remote_side == 1 # remote cell on the "left" of MPI interface
        direction = 1
      else # remote cell on the "right" of MPI interface
        direction = 2
      end
    else # MPI interface in y-direction
      if remote_side == 1 # remote cell on the "left" of MPI interface
        direction = 3
      else # remote cell on the "right" of MPI interface
        direction = 4
      end
    end
    local_element_id = mpi_interfaces.local_element_ids[interface_id]
    local_cell_id = elements.cell_ids[local_element_id]
    remote_cell_id = tree.neighbor_ids[direction, local_cell_id]
    neighbor_ranks[interface_id] = tree.mpi_ranks[remote_cell_id]
    if local_cell_id < remote_cell_id
      global_interface_ids[interface_id] = 2 * ndims(tree) * local_cell_id + direction - 1
    else
      global_interface_ids[interface_id] = (2 * ndims(tree) * remote_cell_id +
                                            opposite_direction(direction) - 1)
    end
  end

  # Get sorted, unique neighbor ranks
  mpi_neighbor_ranks = unique(sort(neighbor_ranks))

  # Sort interfaces by global interface id
  p = sortperm(global_interface_ids)
  neighbor_ranks .= neighbor_ranks[p]
  interface_ids = collect(1:nmpiinterfaces(mpi_interfaces))[p]

  # For each neighbor rank, init connectivity data structures
  mpi_neighbor_interfaces = Vector{Vector{Int}}(undef, length(mpi_neighbor_ranks))
  for (index, d) in enumerate(mpi_neighbor_ranks)
    mpi_neighbor_interfaces[index] = interface_ids[findall(x->(x == d), neighbor_ranks)]
  end

  # Sanity check that we counted all interfaces exactly once
  @assert sum(length(v) for v in mpi_neighbor_interfaces) == nmpiinterfaces(mpi_interfaces)

  return mpi_neighbor_ranks, mpi_neighbor_interfaces
end

# TODO: MPI dimension agnostic
# Initialize MPI data structures
function init_mpi_data_structures(mpi_neighbor_interfaces, ndims, nvars, n_nodes)
  data_size = nvars * n_nodes^(ndims - 1)
  mpi_send_buffers = Vector{Vector{Float64}}(undef, length(mpi_neighbor_interfaces))
  mpi_recv_buffers = Vector{Vector{Float64}}(undef, length(mpi_neighbor_interfaces))
  for index in 1:length(mpi_neighbor_interfaces)
    mpi_send_buffers[index] = Vector{Float64}(undef, length(mpi_neighbor_interfaces[index]) * data_size)
    mpi_recv_buffers[index] = Vector{Float64}(undef, length(mpi_neighbor_interfaces[index]) * data_size)
  end

  mpi_send_requests = Vector{MPI.Request}(undef, length(mpi_neighbor_interfaces))
  mpi_recv_requests = Vector{MPI.Request}(undef, length(mpi_neighbor_interfaces))

  return mpi_send_buffers, mpi_recv_buffers, mpi_send_requests, mpi_recv_requests
end

function rhs!(du, u, t,
              mesh::ParallelTreeMesh{2}, equations,
              initial_condition, boundary_conditions, source_terms,
              dg::DG, cache)
  # Start to receive MPI data
  @timed timer() "start MPI receive" start_mpi_receive!(cache.mpi_cache)

  # Prolong solution to MPI interfaces
  @timed timer() "prolong2mpiinterfaces" prolong2mpiinterfaces!(
    cache, u, mesh, equations, dg)

  # Start to send MPI data
  @timed timer() "start MPI send" start_mpi_send!(
    cache.mpi_cache, mesh, equations, dg, cache)

  # Reset du
  @timed timer() "reset ∂u/∂t" du .= zero(eltype(du))

  # Calculate volume integral
  @timed timer() "volume integral" calc_volume_integral!(
    du, u, mesh,
    have_nonconservative_terms(equations), equations,
    dg.volume_integral, dg, cache)

  # Prolong solution to interfaces
  # TODO: Taal decide order of arguments, consistent vs. modified cache first?
  @timed timer() "prolong2interfaces" prolong2interfaces!(
    cache, u, mesh, equations, dg)

  # Calculate interface fluxes
  @timed timer() "interface flux" calc_interface_flux!(
    cache.elements.surface_flux_values, mesh,
    have_nonconservative_terms(equations), equations,
    dg, cache)

  # Prolong solution to boundaries
  @timed timer() "prolong2boundaries" prolong2boundaries!(
    cache, u, mesh, equations, dg)

  # Calculate boundary fluxes
  @timed timer() "boundary flux" calc_boundary_flux!(
    cache, t, boundary_conditions, mesh, equations, dg)

  # Prolong solution to mortars
  @timed timer() "prolong2mortars" prolong2mortars!(
    cache, u, mesh, equations, dg.mortar, dg)

  # Calculate mortar fluxes
  @timed timer() "mortar flux" calc_mortar_flux!(
    cache.elements.surface_flux_values, mesh,
    have_nonconservative_terms(equations), equations,
    dg.mortar, dg, cache)

  # Finish to receive MPI data
  @timed timer() "finish MPI receive" finish_mpi_receive!(
    cache.mpi_cache, mesh, equations, dg, cache)

  # Calculate MPI interface fluxes
  @timed timer() "MPI interface flux" calc_mpi_interface_flux!(
    cache.elements.surface_flux_values, mesh,
    have_nonconservative_terms(equations), equations,
    dg, cache)

  # Calculate surface integrals
  @timed timer() "surface integral" calc_surface_integral!(
    du, mesh, equations, dg, cache)

  # Apply Jacobian from mapping to reference element
  @timed timer() "Jacobian" apply_jacobian!(
    du, mesh, equations, dg, cache)

  # Calculate source terms
  @timed timer() "source terms" calc_sources!(
    du, u, t, source_terms, equations, dg, cache)

  # Finish to send MPI data
  @timed timer() "finish MPI send" finish_mpi_send!(cache.mpi_cache)

  return nothing
end

function prolong2mpiinterfaces!(cache, u,
                                mesh::ParallelTreeMesh{2},
                                equations, dg::DG)
  @unpack mpi_interfaces = cache

  @threaded for interface in eachmpiinterface(dg, cache)
    local_element = mpi_interfaces.local_element_ids[interface]

    if mpi_interfaces.orientations[interface] == 1 # interface in x-direction
      if mpi_interfaces.remote_sides[interface] == 1 # local element in positive direction
        for j in eachnode(dg), v in eachvariable(equations)
          mpi_interfaces.u[2, v, j, interface] = u[v,          1, j, local_element]
        end
      else # local element in negative direction
        for j in eachnode(dg), v in eachvariable(equations)
          mpi_interfaces.u[1, v, j, interface] = u[v, nnodes(dg), j, local_element]
        end
      end
    else # interface in y-direction
      if mpi_interfaces.remote_sides[interface] == 1 # local element in positive direction
        for i in eachnode(dg), v in eachvariable(equations)
          mpi_interfaces.u[2, v, i, interface] = u[v, i,          1, local_element]
        end
      else # local element in negative direction
        for i in eachnode(dg), v in eachvariable(equations)
          mpi_interfaces.u[1, v, i, interface] = u[v, i, nnodes(dg), local_element]
        end
      end
    end
  end

  return nothing
end

function calc_mpi_interface_flux!(surface_flux_values,
                                  mesh::ParallelTreeMesh{2},
                                  nonconservative_terms::Val{false}, equations,
                                  dg::DG, cache)
  @unpack surface_flux = dg
  @unpack u, local_element_ids, orientations, remote_sides = cache.mpi_interfaces

  @threaded for interface in eachmpiinterface(dg, cache)
    # Get local neighboring element
    element = local_element_ids[interface]

    # Determine interface direction with respect to element:
    if orientations[interface] == 1 # interface in x-direction
      if remote_sides[interface] == 1 # local element in positive direction
        direction = 1
      else # local element in negative direction
        direction = 2
      end
    else # interface in y-direction
      if remote_sides[interface] == 1 # local element in positive direction
        direction = 3
      else # local element in negative direction
        direction = 4
      end
    end

    for i in eachnode(dg)
      # Call pointwise Riemann solver
      u_ll, u_rr = get_surface_node_vars(u, equations, dg, i, interface)
      flux = surface_flux(u_ll, u_rr, orientations[interface], equations)

      # Copy flux to local element storage
      for v in eachvariable(equations)
        surface_flux_values[v, i, direction, element] = flux[v]
      end
    end
  end

  return nothing
end


text:
const eventMouse = Union{Gtk.GdkEventButton, Gtk.GdkEventMotion, Gtk.GdkEventCrossing}
const eventKey = Gtk.GdkEventKey

macro progress(e::Expr, msg, dialog)
    return quote
        if !disableLoadingScreen && $dialog !== nothing
            set_gtk_property!($dialog, :secondary_text, $msg)
            LoadingSpinner.animate($dialog)
            waitForIdle()
        end

        $e
    end |> esc
end

# Stupid hack to make @gtktype macro work from here
function GtkType(name::Symbol)
    expr = Meta.parse("@gtktype $name")
    Base.eval(Gtk, expr)

    return getfield(Gtk, name)
end

# Radio Menu Item

GtkType(:GtkRadioMenuItem)

function RadioMenuItem(name::String, group::Any=C_NULL)
    return Gtk.GtkRadioMenuItem(
        ccall((:gtk_radio_menu_item_new_with_label, Gtk.libgtk), Ptr{GObject}, (Ptr{Nothing}, Ptr{UInt8}), group, name)
    )
end

function RadioMenuItem(name::String, widget::Gtk.GtkRadioMenuItem)
    return Gtk.GtkRadioMenuItem(
        ccall((:gtk_radio_menu_item_new_with_label_from_widget, Gtk.libgtk), Ptr{GObject}, (Ptr{GObject}, Ptr{UInt8}), widget, name)
    )
end

function connectRadioGroupSignal(func::Function, radios::Array{Gtk.GtkRadioMenuItem, 1}; onlyCallbackActive::Bool=true)
    for radio in radios
        @guarded signal_connect(radio, "toggled") do widget
            active = get_gtk_property(widget, :active, Bool)
            index = findfirst(isequal(widget), radios)

            if active && onlyCallbackActive || !onlyCallbackActive
                func(index)
            end
        end
    end
end

function connectRadioGroupSignal(funcs::Array{Function, 1}, radios::Array{Gtk.GtkRadioMenuItem, 1}; onlyCallbackActive::Bool=true)
    for radio in radios
        @guarded signal_connect(radio, "toggled") do widget
            active = get_gtk_property(widget, :active, Bool)
            index = findfirst(isequal(widget), radios)

            if active && onlyCallbackActive || !onlyCallbackActive
                funcs[index]()
            end
        end
    end
end

# Image Menu Item

GtkType(:GtkImageMenuItem)

function setImageMenuItemIcon!(item::Gtk.GtkImageMenuItem, image::Gtk.GtkImage)
    ccall((:gtk_image_menu_item_set_image, Gtk.libgtk), Ptr{Nothing}, (Ptr{GObject}, Ptr{GObject}), item, image)
end

function ImageMenuItem(name::String)
    return Gtk.GtkImageMenuItem(
        ccall((:gtk_image_menu_item_new_with_label , Gtk.libgtk), Ptr{GObject}, (Ptr{UInt8},), name)
    )
end

function ImageMenuItem(name::String, image::Gtk.GtkImage)
    item = ImageMenuItem(name)
    setImageMenuItemIcon!(item, image)

    return item
end

# Check Menu Item

GtkType(:GtkCheckMenuItem)

function CheckMenuItem(name::String, group::Any=C_NULL)
    return Gtk.GtkCheckMenuItem(
        ccall((:gtk_check_menu_item_new_with_label, Gtk.libgtk), Ptr{GObject}, (Ptr{UInt8},), name)
    )
end

# Window Cursor

# Since our objects are kinda badly thrown together we need to cache them to prevent duplication errors
const cursorCache = Dict{Tuple{Gtk.GtkWindow, Ptr{GObject}}, Gtk.GObject}()

# https://developer.gnome.org/gdk3/stable/gdk3-Cursors.html#gdk-cursor-new-from-name Name list
function newCursorFromName(window::Gtk.GtkWindow, name::String)
    ptr = ccall((:gdk_cursor_new_from_name, Gtk.libgdk), Ptr{GObject}, (Ptr{GObject}, Ptr{UInt8}), Gtk.GAccessor.display(window), name)

    if ptr != C_NULL && !haskey(cursorCache, (window, ptr))
        cursorCache[(window, ptr)] = Gtk.GLib.GObjectLeaf(ptr)
    end

    return get(cursorCache, (window, ptr), nothing)
end

function setWindowCursor!(window::Gtk.GtkWindow, cursor::Gtk.GObject)
    return ccall((:gdk_window_set_cursor, Gtk.libgdk), Ptr{GObject}, (Ptr{GObject}, Ptr{GObject}), get_gtk_property(window, :window, Gtk.GObject), cursor)
end

setWindowCursor!(window::Gtk.GtkWindow, name::String) = setWindowCursor!(window, newCursorFromName(window, name))

function getWindowCursor(window::Gtk.GtkWindow)
    ptr = ccall((:gdk_window_get_cursor, Gtk.libgdk), Ptr{GObject}, (Ptr{GObject},), get_gtk_property(window, :window, Gtk.GObject))

    if ptr != C_NULL && !haskey(cursorCache, (window, ptr))
        cursorCache[(window, ptr)] = Gtk.GLib.GObjectLeaf(ptr)
    end

    return get(cursorCache, (window, ptr), newCursorFromName(window, "default"))
end

function changeCursor(func::Function, window::Gtk.GtkWindow, startName::String, stopName::String="")
    startCursor = newCursorFromName(window, startName)
    stopCursor = isempty(stopName) ? getWindowCursor(window) : newCursorFromName(window, stopName)

    setWindowCursor!(window, startCursor)
    func()
    setWindowCursor!(window, stopCursor)
end

# General

function setDefaultDirection(direction::Integer)
    return ccall((:gtk_widget_set_default_direction, Gtk.libgtk), Ptr{}, (Cint,), direction)
end

function initComboBox!(widget::Gtk.GtkComboBoxText, choices::Array{String, 1})
    append!(widget, choices)
    set_gtk_property!(widget, :active, 0)
end

function dropdownString(value::Any)
    return value === nothing ? "nothing" : string(value)
end

function setComboIndex!(widget::Gtk.GtkComboBoxText, choices::Array{T, 1}, item::H; allowCustom::Bool=true) where {T, H}
    if !allowCustom && !(item in choices)
        throw(ArgumentError("The selected value does not exist for the given choices"))
    end

    item = dropdownString(item)
    choices = dropdownString.(choices)

    if !(item in choices)
        push!(choices, item)
        push!(widget, item)
    end

    index = findfirst(isequal(item), choices)
    Gtk.GLib.@sigatom set_gtk_property!(widget, :active, something(index, 1) - 1)

    return index
end

function setEntryText!(entry::Gtk.GtkEntryLeaf, value::Any; updatePlaceholder::Bool=true)
    text = string(value)

    Gtk.GLib.@sigatom GAccessor.text(entry, text)

    if updatePlaceholder
        Gtk.GLib.@sigatom GAccessor.placeholder_text(entry, text)
    end
end

function convertString(str::String, as::Type=String)
    if as == String
        return str

    elseif as == Nothing
        return nothing

    elseif as == Char
        if length(str) == 1
            return str[1]

        else
            error("Expected string with length 1, got $(length(str))")
        end

    elseif as == Bool
        return parse(as, str)

    elseif as <: Integer
        return parse(Int, str)

    elseif as <: AbstractFloat
        return parse(Float64, str)

    else
        error("Unsupported return type")
    end
end

function getEntryText(entry::Gtk.GtkEntryLeaf, as::Type=String)
    str = Gtk.GLib.@sigatom Gtk.bytestring(GAccessor.text(entry))

    return convertString(str, as)
end

function setTextViewText!(textView::Gtk.GtkTextViewLeaf, value::String="")
    GAccessor.buffer(textView, TextBuffer(text=value))
end

function getTextViewText(textView::Gtk.GtkTextViewLeaf)
    return join(collect(GAccessor.buffer(textView)))
end

function addProviderForScreen(screen::Any, provider::Gtk.GtkCssProviderLeaf)
    return ccall((:gtk_style_context_add_provider_for_screen, Gtk.libgtk), Nothing, (Ptr{Nothing}, Ptr{GObject}, Cuint), screen, provider, 1)
end

function setDialogFolder(dialog::GObject, path::String)
    return ccall((:gtk_file_chooser_set_current_folder, Gtk.libgtk), Cint, (Ptr{GObject}, Ptr{UInt8}), dialog, path)
end

function setDialogFilename(dialog::GObject, path::String)
    return ccall((:gtk_file_chooser_set_current_name, Gtk.libgtk), Cint, (Ptr{GObject}, Ptr{UInt8}), dialog, path)
end

function getProgressDialog(; title::String="", description::String="", filename::String="", stylesheet::String="", parent::Any=GtkNullContainer(), buttons::Tuple{}=())
    pixbuf = isempty(filename) ? Pixbuf(width=1, height=1, has_alpha=true) : Pixbuf(filename=filename, preserve_aspect_ratio=true)
    image = Image(pixbuf, name="dialog_image")

    dlg = GtkMessageDialog(title, (),
        Gtk.GtkDialogFlags.DESTROY_WITH_PARENT, Gtk.GtkMessageType.INFO, parent,
        secondary_text=description, image=image
    )

    if !isempty(stylesheet)
        screen = Gtk.GAccessor.screen(dlg)
        provider = CssProviderLeaf(filename=stylesheet)

        addProviderForScreen(screen, provider)
    end

    return dlg
end

function setProgressDialogPixbuf!(dlg::Gtk.GObject, pixbuf::Union{Gtk.GdkPixbuf, Ptr{GObject}})
    image = get_gtk_property(dlg, :image, GtkImage)

    set_gtk_property!(image, :visible, true)
    set_gtk_property!(image, :pixbuf, pixbuf)
end

useNativeFileDialogs = isdefined(Gtk, :save_dialog_native) && Gtk.libgtk_version >= v"3.20.0"

# Copy of Gtk.jl definition, but calls a function before running the widget
function infoDialog(func::Function, message::String, parent::Gtk.GtkWindow=Ahorn.window)
    w = GtkMessageDialogLeaf(ccall((:gtk_message_dialog_new, Gtk.libgtk), Ptr{Gtk.GObject},
        (Ptr{Gtk.GObject}, Cint, Cint, Cint, Ptr{UInt8}),
        parent, Gtk.GtkDialogFlags.DESTROY_WITH_PARENT,
        Gtk.GtkMessageType.INFO, Gtk.GtkButtonsType.CLOSE, C_NULL))

    set_gtk_property!(w, :text, message)

    func(w)

    run(w)
    Gtk.destroy(w)
end

infoDialog(message::String, parent::Gtk.GtkWindow=Ahorn.window) = infoDialog(w -> (), message, parent)
topMostInfoDialog(message::String, parent::Gtk.GtkWindow=Ahorn.window) = infoDialog(w -> GAccessor.keep_above(w, true), message, parent)

# Copy of Gtk.jl definition, but calls a function before running the widget
askDialog(message::String, parent::Gtk.GtkWindow=Ahorn.window) = askDialog(message, "No", "Yes", parent)
askDialog(func::Function, message::String, parent::Gtk.GtkWindow=Ahorn.window) = askDialog(func, message, "No", "Yes", parent)

function askDialog(func::Function, message::String, noText::String, yesText::String, parent::Gtk.GtkWindow=Ahorn.window)
    dlg = GtkMessageDialog(message, ((noText, 0), (yesText, 1)),
            Gtk.GtkDialogFlags.DESTROY_WITH_PARENT, Gtk.GtkMessageType.QUESTION, parent)

    func(dlg)

    response = run(dlg)
    Gtk.destroy(dlg)

    return response == 1
end

askDialog(message::String, noText::String, yesText::String, parent::Gtk.GtkWindow=Ahorn.window) = askDialog(dlg -> (), message, noText, yesText, parent)
topMostAskDialog(message::String, parent::Gtk.GtkWindow=Ahorn.window) = askDialog(dlg -> GAccessor.keep_above(dlg, true), message, parent)

# Modified version of standard open dialog
function openDialog(title::AbstractString, parent=GtkNullContainer(), filters::Union{AbstractVector, Tuple}=String[]; folder::String="", native::Bool=useNativeFileDialogs, kwargs...)
    local dlg

    if native
        dlg = GtkFileChooserNative(title, parent, GConstants.GtkFileChooserAction.OPEN, "_Open", "_Cancel"; kwargs...)

    else
        dlg = GtkFileChooserDialog(title, parent, GConstants.GtkFileChooserAction.OPEN, (
                ("_Cancel", GConstants.GtkResponseType.CANCEL),
                ("_Open",   GConstants.GtkResponseType.ACCEPT)
            ); kwargs...)
    end

    dlgp = GtkFileChooser(dlg)

    if !isempty(filters)
        Gtk.makefilters!(dlgp, filters)
    end

    if !isempty(folder)
        setDialogFolder(dlgp, folder)
    end

    response = run(dlg)
    multiple = get_gtk_property(dlg, :select_multiple, Bool)
    local selection

    if response == GConstants.GtkResponseType.ACCEPT
        if multiple
            filename_list = ccall((:gtk_file_chooser_get_filenames, Gtk.libgtk), Ptr{Gtk._GSList{String}}, (Ptr{GObject},), dlgp)
            selection = String[f for f in Gtk.GList(filename_list, true)]

        else
            selection = Gtk.bytestring(GAccessor.filename(dlgp))
        end

    else
        if multiple
            selection = String[]

        else
            selection = ""
        end
    end

    native ? Gtk.GLib.gc_unref(dlg) : Gtk.destroy(dlg)

    return selection
end

# Modified version of standard save dialog
function saveDialog(title::AbstractString, parent=GtkNullContainer(), filters::Union{AbstractVector, Tuple}=String[]; folder::String="", filename::String="", warnOverwrite::Bool=true, native::Bool=useNativeFileDialogs, kwargs...)
    local dlg

    if native
        dlg = GtkFileChooserNative(title, parent, GConstants.GtkFileChooserAction.SAVE, "_Save", "_Cancel"; kwargs...)

    else
        dlg = GtkFileChooserDialog(title, parent, GConstants.GtkFileChooserAction.SAVE, (
                ("_Cancel", GConstants.GtkResponseType.CANCEL),
                ("_Save",   GConstants.GtkResponseType.ACCEPT)
            ); kwargs...)
    end

    dlgp = GtkFileChooser(dlg)

    if !isempty(folder)
        setDialogFolder(dlgp, folder)
    end

    if !isempty(filename)
        setDialogFilename(dlgp, filename)
    end

    if !isempty(filters)
        Gtk.makefilters!(dlgp, filters)
    end

    ccall((:gtk_file_chooser_set_do_overwrite_confirmation, Gtk.libgtk), Nothing, (Ptr{GObject}, Cint), dlg, warnOverwrite)
    response = run(dlg)

    if response == GConstants.GtkResponseType.ACCEPT
        selection = Gtk.bytestring(GAccessor.filename(dlgp))

    else
        selection = ""
    end

    native ? Gtk.GLib.gc_unref(dlg) : Gtk.destroy(dlg)

    return selection
end

function scalePixbufSimple(pixbuf::Gtk.GdkPixbuf, width::Integer=width(pixbuf), height::Integer=height(pixbuf), method::Int32=GdkInterpType.NEAREST)
    return ccall((:gdk_pixbuf_scale_simple, Gtk.libgdkpixbuf), Ptr{GObject}, (Ptr{GObject}, Cint, Cint, Cint), pixbuf, width, height, method)
end

function pixbufFromSurface(surface::Cairo.CairoSurface, alpha::Bool=true)
    return Pixbuf(data=[Ahorn.ARGB(Ahorn.argb32ToRGBATuple(c)...) for c in getSurfaceData(surface)], has_alpha=alpha)
end

getSurfaceData(surface::Cairo.CairoSurfaceImage) = surface.data

function getSurfaceData(surface::Cairo.CairoSurface)
    imageSurface = CairoImageSurface(zeros(UInt32, floor(Int, height(surface)), floor(Int, width(surface))), Cairo.FORMAT_ARGB32)
    imageSurfaceCtx = CairoContext(imageSurface)

    drawImage(imageSurfaceCtx, surface, 0, 0)

    return imageSurface.data
end

const cairoGCCache = Dict{Int, Cairo.CairoContext}()

function getSurfaceContext(surface::Cairo.CairoSurface)::Cairo.CairoContext
    key = Int(surface.ptr)

    if get(cairoGCCache, key, C_NULL) == C_NULL
        delete!(cairoGCCache, surface)
    end

    return get!(cairoGCCache, key) do
        creategc(surface)
    end
end

function deleteSurface(surface::Cairo.CairoSurface)
    key = Int(surface.ptr)
    if haskey(cairoGCCache, key)
        ctx = cairoGCCache[key]

        delete!(cairoGCCache, key)

        Cairo.destroy(ctx)
        Cairo.destroy(surface)
    end
end

# Borrowed from Gtk.jl repo, build is not installable yet
function g_timeout_add(interval::Integer, cb::Function, user_data::CT) where CT
    callback = @cfunction($cb, Cint, (Ref{CT},) )

    ref, deref = Gtk.GLib.gc_ref_closure(user_data)

    return ccall((:g_timeout_add_full, Gtk.libglib), Cint, (Cint, UInt32, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}), 0, UInt32(interval), callback, ref, deref)
end

# Borrowed from Gtk.jl repo, build is not installable yet
function g_idle_add(cb::Function, user_data::CT) where CT
    callback = @cfunction($cb, Cint, (Ref{CT},) )
    ref, deref = Gtk.GLib.gc_ref_closure(user_data)

    return ccall((:g_idle_add_full, Gtk.libglib), Cint, (Cint, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}), 0, callback, ref, deref)
end

function addIdleCallback(cb::Function, data::CT, timeout::Number=-1) where CT
    if timeout == -1
        g_idle_add(data -> Cint(cb(data)), data)

    else
        g_timeout_add(round(Int, timeout * 1000), data -> Cint(cb(data)), data)
    end
end

waitForIdle() = addIdleCallback(data -> false, nothing)

MenuItemsTypes = Union{Gtk.GtkMenuItem, Gtk.GtkImageMenuItem}

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
immutable WaveletLayer <: Mocha.Layer
    bank::AbstractBank
    bottoms::Vector{Symbol}
    name::AbstractString
    tops::Vector{Symbol}
end

function WaveletLayer( ;
        bank::AbstractBank = NullBank(),
        bottoms::Vector{Symbol} = Symbol[],
        name::AbstractString = "wavelets",
        tops::Vector{Symbol} = Symbol[])
    @assert !isa(bank, NullBank)
    @assert length(bottoms) == 1
    @assert length(tops) == 1
    return WaveletLayer(bank, bottoms, name, tops)
end

Mocha.can_do_bp(::WaveletLayer) = false
Mocha.has_neuron(::WaveletLayer) = false
Mocha.has_param(::WaveletLayer) = false


text:
using MakieLayout
using Makie

function axhline!(la::LAxis, y; kwargs...)
    points = lift(la.limits) do lims
        x1 = lims.origin[1]
        x2 = lims.origin[1] + lims.widths[1]
        Point2f0.([(x1, y), (x2, y)])
    end

    lines!(la, points; xautolimit=false, kwargs...)
end

function axvline!(la::LAxis, x; kwargs...)
    points = lift(la.limits) do lims
        y1 = lims.origin[2]
        y2 = lims.origin[2] + lims.widths[2]
        Point2f0.([(x, y1), (x, y2)])
    end

    lines!(la, points; yautolimit=false, kwargs...)
end

function abline!(la::LAxis, a, b; kwargs...)
    points = lift(la.limits) do lims
        x1 = lims.origin[1]
        x2 = lims.origin[1] + lims.widths[1]
        Point2f0.([(x1, a * x1 + b), (x2, a * x2 + b)])
    end

    lines!(la, points; xautolimit=false, yautolimit=false, kwargs...)
end

function dynamicfunc!(la::LAxis, func; n=nothing, kwargs...)

    points = lift(la.limits, la.layoutnodes.computedbbox) do lims, bbox

        if isnothing(n)
            xwidth = width(bbox)
            n = Int(round(xwidth))
        end

        x1 = lims.origin[1]
        x2 = lims.origin[1] + lims.widths[1]
        xs = LinRange(x1, x2, n)
        ys = func.(xs)
        [Point2f0(x, y) for (x, y) in zip(xs, ys)]
    end

    lines!(la, points; xautolimit=false, yautolimit=false, kwargs...)
end

function dynamicimage!(la::LAxis, func; n=300, kwargs...)
    pixels = lift(la.limits) do lims
        x1, y1 = lims.origin
        x2, y2 = lims.origin .+ lims.widths
        xs = LinRange(x1, x2, n)
        ys = LinRange(y1, y2, n)

        pixels = Float32[func(x, y) for x in xs, y in ys]
    end

    xrange = lift(lims -> lims.origin[1] .+ [0, lims.widths[1]], la.limits)
    yrange = lift(lims -> lims.origin[2] .+ [0, lims.widths[2]], la.limits)

    image!(la, xrange, yrange, pixels; xautolimit=false, yautolimit=false, kwargs...)
end

begin
    scene = Scene(resolution = (1000, 1000));
    screen = display(scene)
    campixel!(scene);

    nrows = 1
    ncols = 1

    maingl = GridLayout(scene, 1, 1, alignmode=Outside(30))

    gridgl = maingl[1, 1] = GridLayout(
        nrows, ncols)

    la = gridgl[1, 1] = LAxis(scene)

    # scatter!(la, rand(100, 2) .* 50, markersize=20)
    axhline!(la, 0)
    axvline!(la, 0)
    abline!(la, 1, 0)

    dynamicfunc!(la, sin, color=:red)
    dynamicfunc!(la, cos, color=:blue)

    nothing
end

# save("layout.png", scene)


text:
"""
Runs Impute 5 on chromosome 20 over 3 chunks.
Each chunk imputes 20 million SNPs using 4 or 3 threads.
Each script runs a different chunk and are initialized simultaneously.
"""
function impute5_1000genomes_chr20(
    reffile::AbstractString,
    tgtfile::AbstractString,
    )
    chr = 20                # chromosome 20
    allale_start = 1        # first SNP is 60479
    allele_end = 64000000   # last SNP is 62965167
    chunk_size = 22000000   # each chunk imputes 22cM
    nchunks = 3             # chunk 1 imputes first 22cM region, chunk 2 imputes next 22cM ...etc
    cthreads = [4, 3, 3]    # chunk 1 uses 4 threads, chunk 2 uses 3 threads...etc

    # start 3 processes, each imputing different regions, using 3-4 threads
    Threads.@threads for i in 1:nchunks
        cur_start = (i - 1) * chunk_size + 1
        cur_end = i == nchunks ? allele_end : chunk_size * i
        cur_threads = cthreads[i]
        run(`/home/biona001/impute5_v1.1.4/impute5_1.1.4_static
            --h $reffile
            --g $tgtfile
            --r $chr:$cur_start-$cur_end
            --o impute5.chr$chr.result.chunk$i.vcf.gz
            --threads $cur_threads`)
    end
end

"""
Runs Impute 5 on chromosome 10 over 7 chunks.
Each chunk imputes 20 million SNPs using 1 or 2 thread.
"""
function impute5_1000genomes_chr10(
    reffile::AbstractString,
    tgtfile::AbstractString,
    )
    chr = 10                # chromosome 10
    allale_start = 1        # first SNP is 60494
    allele_end = 140000000  # last SNP is 135524373
    chunk_size = 20000000   # each chunk imputes 20cM
    nchunks = 7            # chunk 1 imputes first 20cM region, chunk 2 imputes next 20cM ...etc
    cthreads = [2, 2, 2, 1, 1, 1, 1]  # each chunk runs 1 threads

    # run each chunk with 1 thread
    Threads.@threads for i in 1:nchunks
        cur_start = (i - 1) * chunk_size + 1
        cur_end = i == nchunks ? allele_end : chunk_size * i
        cur_threads = cthreads[i]
        run(`/home/biona001/impute5_v1.1.4/impute5_1.1.4_static
            --h $reffile
            --g $tgtfile
            --r $chr:$cur_start-$cur_end
            --o impute5.chr$chr.result.chunk$i.vcf.gz
            --threads $cur_threads`)
    end
end

chr = parse(Int, ARGS[1])
reffile = ARGS[2]
tgtfile = ARGS[3]

if Threads.nthreads() == 1
    error("start Julia with >1 thread for this script to work!")
end

if chr == 20
    _, t, _, _, _ = @timed impute5_1000genomes_chr20(reffile, tgtfile)
elseif chr == 10
    _, t, _, _, _ = @timed impute5_1000genomes_chr10(reffile, tgtfile)
end

println("Total elapsed time = $t seconds")

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
const registernames = ["w", "x", "y", "z"]
struct AsmLine
    op::String
    idx1::Int
    idx2::Int
    n::Int
    function AsmLine(str::String)
        s = split(strip(str), r"\s+")
        b = length(s) > 2 ? (s[3] in registernames ? s[3] : tryparse(Int, s[3])) : nothing
        idx1 = findfirst(registernames .== s[2])
        idx2 = (b isa AbstractString) ? findfirst(registernames .== b) : -1
        n = b isa Integer ? b : typemax(Int)
        return new(s[1], idx1, idx2, n)
    end
end

""" Speedup: if we fail to pop a 26-multiplier due to bad input, z becomes large, so drop those"""
function day24(maxz = 26^4)
    rstates = [[[0, 0, 0, 0], [0, 0]]]
    lines = readlines("AoCdata/AoC_2021_day24.txt")
    for (i, line) in enumerate(lines)
        asm = AsmLine(line)
        if asm.op == "inp"
            d = Dict{Vector{Int},Vector{Int}}()
            for (reg, (mn, mx)) in rstates
                (w, x, y, z) = reg
                k = [0, x, y, z]
                if haskey(d, k)
                    if d[k][1] > mn || d[k][1] == 0
                        d[k][1] = mn
                    end
                    if d[k][2] < mx
                        d[k][2] = mx
                    end
                elseif z <= maxz
                    d[k] = [mn, mx]
                end
            end
            rstates = [[k, v] for (k, v) in d]
            newstates = empty(rstates)
            for ((w, x, y, z), (mn, mx)) in rstates, i in 1:9
                push!(newstates, [[i, x, y, z], [mn * 10 + i, mx * 10 + i]])
            end
            rstates = newstates
        elseif asm.op == "add"
            for i in eachindex(rstates)
                rstates[i][1][asm.idx1] += asm.idx2 > 0 ? rstates[i][1][asm.idx2] : asm.n
            end
        elseif asm.op == "mul"
            for i in eachindex(rstates)
                rstates[i][1][asm.idx1] *= asm.idx2 > 0 ? rstates[i][1][asm.idx2] : asm.n
            end
        elseif asm.op == "div"
            for i in eachindex(rstates)
                rstates[i][1][asm.idx1] ÷= asm.idx2 > 0 ? rstates[i][1][asm.idx2] : asm.n
            end
        elseif asm.op == "mod"
            for i in eachindex(rstates)
                rstates[i][1][asm.idx1] %= asm.idx2 > 0 ? rstates[i][1][asm.idx2] : asm.n
            end
        elseif asm.op == "eql"
            for i in eachindex(rstates)
                p = rstates[i][1][asm.idx1]
                q = asm.idx2 > 0 ? rstates[i][1][asm.idx2] : asm.n
                rstates[i][1][asm.idx1] = p == q
            end
        else
            error("unknown op $op in line $line")
        end
    end
    filter!(x -> x[1][4] == 0, rstates)
    return maximum(x[2][2] for x in rstates), minimum(x[2][1] for x in rstates)
end

part = day24()
println("Part 1: ", part[1])
println("Part 2: ", part[2])


text:
module IOLoggingTests

using IOLogging
using Logging
using Logging: Debug, Info, Warn, Error, BelowMinLevel, with_logger, min_enabled_level
using Test

function singleLineLogging()
    @debug "DEBUG"
    @info "INFO"
    @warn "WARN"
    @error "ERROR"
end

function multiLineLogging()
    @debug "debug_line_1\ndebug_line_2"
    @info "info_line_1\ninfo_line_2"
    @warn "warn_line_1\nwarn_line_2"
    @error "error_line_1\nerror_line_2"
end

function keywordLogging()
    deb = "debug"
    inf = "info"
    war = "warn"
    err = "error"

    @debug "DEBUG" deb
    @info "INFO" inf
    @warn "WARN" war
    @error "ERROR" err
end

function setupTestlog()
    if isfile("appFlush.log")
        rm("appFlush.log")
    end
    open("appFlush.log", "w") do file
        println(file, "testlog")
    end
end

@testset "IOLogging All Tests" begin
@testset "IOLogger" begin
    @testset "Assertions" begin
        @test_nowarn IOLogger()
        logger = IOLogger(Dict(BelowMinLevel => stdout))
        @test min_enabled_level(logger) === BelowMinLevel
    end
    @testset "Logging" begin
        buf = IOBuffer()
        log = IOLogger(Dict(Info => buf))
        @testset "Singleline No Debug" begin
            with_logger(singleLineLogging, log)
            res = String(take!(buf))
            @test !occursin(r"(DEBUG)|(Debug)", res)
            @test occursin(r"\[Info::[\d\-:T\.]*\][ a-zA-Z@\[\]\d\.\\\/:]* - INFO\n", res)
            @test occursin(r"\[Warn::[\d\-:T\.]*\][ a-zA-Z@\[\]\d\.\\\/:]* - WARN\n", res)
            @test occursin(r"\[Error::[\d\-:T\.]*\][ a-zA-Z@\[\]\d\.\\\/:]* - ERROR\n", res)
        end
        @testset "Multiline No Debug" begin
            with_logger(multiLineLogging, log)
            res = String(take!(buf))
            @test (!occursin(r"(DEBUG)|(Debug)", res) &&
                    !occursin(r"debug_line_1\n", res) &&
                    !occursin(r"debug_line_2\n", res))
            @test (occursin(r"\[Info::[\d\-:T\.]*\][ a-zA-Z@\[\]\d\.\\\/:]*\n", res) &&
                    occursin(r"info_line_1\n", res) &&
                    occursin(r"info_line_2\n", res))
            @test (occursin(r"\[Warn::[\d\-:T\.]*\][ a-zA-Z@\[\]\d\.\\\/:]*\n", res) &&
                    occursin(r"warn_line_1\n", res) &&
                    occursin(r"warn_line_2\n", res))
            @test (occursin(r"\[Error::[\d\-:T\.]*\][ a-zA-Z@\[\]\d\.\\\/:]*\n", res) &&
                    occursin(r"error_line_1\n", res) &&
                    occursin(r"error_line_2\n", res))
        end
        @testset "Keyword Logging No Debug" begin
            with_logger(keywordLogging, log)
            res = String(take!(buf))

            @test !occursin(r"deb = debug", res)
            @test occursin(r"inf = info", res)
            @test occursin(r"war = warn", res)
            @test occursin(r"err = err", res)
        end
        @testset "Logging between Levels No Debug" begin
            infLog = IOBuffer()
            errLog = IOBuffer()
            logger = IOLogger(Dict(Info => infLog, Error => errLog))
            with_logger(singleLineLogging, logger)

            infres = String(take!(infLog)) # should only contain messages from INFO and WARN
            @test  (!occursin(r"DEBUG", infres) &&
                     occursin(r"INFO", infres) &&
                     occursin(r"WARN", infres) &&
                    !occursin(r"ERROR", infres))

            errres = String(take!(errLog)) # should only contain messages from ERROR
            @test  (!occursin(r"DEBUG", errres) &&
                    !occursin(r"INFO", errres) &&
                    !occursin(r"WARN", errres) &&
                     occursin(r"ERROR", errres))
        end
        @testset "Message Limits" begin
            with_logger(log) do
                for _ in 1:20
                    @info "message_Limit_test" maxlog=10
                end
            end
            res = String(take!(buf))
            @test count(c -> c == '\n', res) == 10

            # This is a different log position with the same message
            with_logger(log) do
                for _ in 1:20
                    @info "message_Limit_test" maxlog=10
                end
            end
            res = String(take!(buf))
            @test count(c -> c == '\n', res) == 10
        end
    end
end

@testset "FileLogger" begin
    mktempdir(@__DIR__) do dir
        defaultLog = "default.log"
        cd(dir) do
            @testset "Assertions" begin
                @test_nowarn FileLogger()
                file = "assertions.log"
                logger = FileLogger(Dict(BelowMinLevel => file))
                @test !isfile(file) # make sure the files don't exist until needed
                # FIXME: @test if the file exists after logging
                @test min_enabled_level(logger) === BelowMinLevel
            end
            @testset "Appending && Flushing" begin
                @testset "Regular Use" begin
                    log = FileLogger(Dict(Info => "appFlush.log"))
                    setupTestlog()
                    with_logger(log) do
                        @info "infolog"
                    end

                    lines = readlines("appFlush.log", keep = true)
                    @test lines[1] == "testlog\n"
                    @test occursin(r"\[Info::[\d\-:T\.]*\][ a-zA-Z@\[\]\d\.\\\/:]* - infolog\n", lines[2])
                end

                @testset "No Flush" begin
                    log = FileLogger(Dict(Info => "appFlush.log"), flush = false)
                    setupTestlog()
                    with_logger(log) do
                        @info "infolog"
                    end

                    lines = readlines("appFlush.log", keep = true)
                    @test length(lines) == 1
                    @test lines[1] == "testlog\n"
                end

                @testset "No Append" begin
                    log = FileLogger(Dict(Info => "appFlush.log"), append = false)
                    setupTestlog()
                    with_logger(log) do
                        @info "infolog"
                    end

                    lines = readlines("appFlush.log", keep = true)
                    @test length(lines) == 1
                    @test occursin(r"\[Info::[\d\-:T\.]*\][ a-zA-Z@\[\]\d\.\\\/:]* - infolog\n", lines[1])
                end
            end
            @testset "Logging" begin
                log = FileLogger()
                @testset "Singleline No Debug" begin
                    with_logger(singleLineLogging, log)
                    res = string(readlines(defaultLog, keep = true)...)

                    @test !occursin(r"(DEBUG)|(Debug)", res)
                    @test occursin(r"\[Info::[\d\-:T\.]*\][ a-zA-Z@\[\]\d\.\\\/:]* - INFO\n", res)
                    @test occursin(r"\[Warn::[\d\-:T\.]*\][ a-zA-Z@\[\]\d\.\\\/:]* - WARN\n", res)
                    @test occursin(r"\[Error::[\d\-:T\.]*\][ a-zA-Z@\[\]\d\.\\\/:]* - ERROR\n", res)
                end
                @testset "Multiline No Debug" begin
                    with_logger(multiLineLogging, log)
                    res = string(readlines(defaultLog, keep = true)...)

                    @test (!occursin(r"(DEBUG)|(Debug)", res) &&
                           !occursin(r"debug_line_1\n", res) &&
                           !occursin(r"debug_line_2\n", res))
                    @test (occursin(r"\[Info::[\d\-:T\.]*\][ a-zA-Z@\[\]\d\.\\\/:]*\n", res) &&
                           occursin(r"info_line_1\n", res) &&
                           occursin(r"info_line_2\n", res))
                    @test (occursin(r"\[Warn::[\d\-:T\.]*\][ a-zA-Z@\[\]\d\.\\\/:]*\n", res) &&
                           occursin(r"warn_line_1\n", res) &&
                           occursin(r"warn_line_2\n", res))
                    @test (occursin(r"\[Error::[\d\-:T\.]*\][ a-zA-Z@\[\]\d\.\\\/:]*\n", res) &&
                           occursin(r"error_line_1\n", res) &&
                           occursin(r"error_line_2\n", res))
                end
                @testset "Keyword Logging No Debug" begin
                    with_logger(keywordLogging, log)
                    res = string(readlines(defaultLog, keep = true)...)

                    @test !occursin(r"deb = debug", res)
                    @test occursin(r"inf = info", res)
                    @test occursin(r"war = warn", res)
                    @test occursin(r"err = err", res)
                end
                @testset "Logging between Levels No Debug" begin
                    infLog = "infLog"
                    errLog = "errLog"
                    logger = FileLogger(Dict(Info => infLog, Error => errLog))
                    with_logger(singleLineLogging, logger)
                    infres = string(readlines(infLog, keep = true)...)
                    errres = string(readlines(errLog, keep = true)...)

                    @test  (!occursin(r"DEBUG", infres) &&
                            occursin(r"INFO", infres) &&
                            occursin(r"WARN", infres) &&
                            !occursin(r"ERROR", infres))

                    @test  (!occursin(r"DEBUG", errres) &&
                            !occursin(r"INFO", errres) &&
                            !occursin(r"WARN", errres) &&
                            occursin(r"ERROR", errres))
                end

                @testset "Message Limits" begin
                    with_logger(log) do
                        for _ in 1:20
                            @info "message_Limit_test" maxlog=10
                        end
                    end
                    res = string(readlines(defaultLog, keep = true)...)
                    matches = collect(eachmatch(r"message_Limit_test", res))
                    @test length(matches) == 10

                    # This is a different log position with the same message
                    with_logger(log) do
                        for _ in 1:20
                            @info "message_Limit_test" maxlog=10
                        end
                    end
                    res = string(readlines(defaultLog, keep = true)...)
                    matches = collect(eachmatch(r"message_Limit_test", res))
                    @test length(matches) == 20
                end
            end
        end
    end
end
end

end


text:
# These examples are documented by the LCP notebooks from
# https://github.com/QuantEcon/SimpleDifferentialOperators.jl
# In general, these are optimal stopping problems solved
# as linear complementarity problems.

using Parameters, DiffEqOperators, Random, Suppressor, Test, PATHSolver

# Seed
Random.seed!(1793)

function LCP_split(S_)
  # setup
    u1(x) = x^2
    ρ = 0.75
    μ1(x) = -x
    σ = 0.1
    M = 300
    x = range(-5., 5., length = M)

    # construct operator and non-operator LCP objects
    BC = RobinBC((0., 1., 0.), (0., 1., 0.), x[2] - x[1], 1)
    L1 = UpwindDifference(1, 1, x[2] - x[1], M, map(t -> μ1(t), x)) * BC
    L2 = σ^2 / 2 * CenteredDifference(2, 2, x[2] - x[1], M) * BC
    u_vec = u1.(x)

    return (L1 = Array(L1)[1], L2 = Array(L2)[1], S = S_ .* ones(M), u = u_vec, ρ = ρ)
end

@testset "Split Case, Nonzero S" begin
  # setup
  @unpack L1, L2, S, u, ρ = LCP_split(0.125)
  g = z -> ρ .* z - L1 * z - L2 * z -
      u + ρ .* S - L1 * S - L2 * S

  n = 300
  lb = zeros(n)
  ub = 300. .* ones(n) # a reasonable guess?
  options(convergence_tolerance = 1e-14, output = :no,
          time_limit = 600) # 10 minute budget
  exit_code, sol_z, sol_g = @suppress solveLCP(g, lb, ub);
  # tests
  @test exit_code == :StationaryPointFound
  @test sol_z[3] ≈ 8.774937148286833 atol = 1e-5
  @test sol_z[120] ≈ 0.3021168981772892 atol = 1e-5
  @test sol_z[270] ≈ 5.729541001488482 atol = 1e-5
  @test sol_g[end-1] ≈ 3.197442310920451e-14 atol = 1e-5
end

@testset "Split Case, Zero S" begin
  # setup
  @unpack L1, L2, S, u, ρ = LCP_split(0.)
  f = z -> ρ .* z - L1 * z - L2 * z -
      u + ρ .* S - L1 * S - L2 * S

  # @unpack L, B, S, q = LCP_split(0.)
  # f = z -> B*z + q
  n = 300
  lb = zeros(n)
  ub = 300*ones(n) # a reasonable guess?
  options(convergence_tolerance = 1e-12, output = :no,
          time_limit = 600) # 10 minute budget
  exit_code, sol_z, sol_f = @suppress solveLCP(f, lb, ub);

  # tests
  @test exit_code == :StationaryPointFound
  @test sol_z[3] ≈ 8.888461348456772 atol = 1e-5
  @test sol_z[76] ≈ 2.279767804635279 atol = 1e-5
  @test sol_z[150] ≈ 0.005770703117189083 atol = 1e-5
  @test sol_z[269] ≈ 5.744079221450249 atol = 1e-5
  @test sol_f[123] ≈ 1.7763568394002505e-15 atol = 1e-5
end

nothing


text:
using ParPoly
# using Base.Test

# using SIMD
# x = 2.0
# c = randn(Float64, 16)

# g(x,c) = @estrin   x c[1] c[2] c[3] c[4] c[5] c[6] c[7] c[8]
# f(x,c) = @evalpoly x c[1] c[2] c[3] c[4] c[5] c[6] c[7] c[8]

# @test g(x,c) ≈ f(x,c)

# @code_native g(x,c)

# macroexpand(:(@estrin x c[1] c[2] c[3] c[4] c[5] c[6] c[7] c[8]))

# @time f(x,c)
# @time g(x,c)