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

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 iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 20bb3256-f1bb-4b3b-b410-89969443fd4d
begin
    using Pkg
    Pkg.activate("../../")
    using Revise
	# maybe needed since we added a few dependencies
	Pkg.instantiate()
end

# ╔═╡ db61ecb1-6ed9-44d9-b8e4-b50622efbac1
using TestImages, ImgProcMic, ImageShow, Random, PlutoTest, Plots, FourierTools, Statistics, FFTW, IndexFunArrays, SpecialFunctions, PlutoUI, Noise, Optim, ForwardDiff

# ╔═╡ 82726a23-060c-4c7d-91c1-e78a22fb98d6
md"# 1 Fourier Transforms

In the last homework we learnt that FFTW calculates the Fourier transform such that the DC frequency (zero frequency) is at the first entry of the array.

To shift the frequency to the center, we normally use `fftshift`.
However, another issue is, that the FFT also interpretes the center of the array at the first index position.
So for example, we would expect that the Fourier transform of the delta peak is a constant array!

As you can see below, the Fourier transform is only a constant array if the delta peak is at the fist entry
"

# ╔═╡ ebe90686-c72f-4a42-a644-f06babe00f41
# FFT of delta peak results in constant array
fft([1.0, 0.0, 0.0])

# ╔═╡ 5053c100-5f27-4f2f-acbe-e12c9070feac
# no constant array since delta is not at the first entry
fft([0.0, 1.0, 0.0])

# ╔═╡ 361a7a34-d70b-4597-948c-80ed913a562b
# no constant array since delta is not at the first entry
fft([0.0, 0.0, 1.0])

# ╔═╡ 653a94c9-cb89-47df-8ba5-0ce683440b78
md"
The reason for this behaviour is that the second example is a shifted array of the first one. Hence, in Fourier space this corresponds to a phase ramp!

In the last homework we programmed `ffts(x) = fftshift(fft(x))` and `iffts(x) = ifft(ifftshift(x))`.

We now also introduce `my_ft` and `my_ift`.

The qualitative meaning of the different conventions for the Fourier transform of a signal with length `N` is:
* `fft`: center in real space is at `(1,)` and in Fourier space at `(1,)`
* `ffts`: center in real space is at `(1,)` and in Fourier space at `(N ÷ 2 + 1,)`
* `ft`: center in real space is at `(N ÷ 2 + 1,)` and in Fourier space at `(N ÷ 2 + 1,)`

"

# ╔═╡ b65e0745-82c6-48e1-a956-fa71583dad15
begin
	my_ft(x) = fftshift(fft(ifftshift(x)))
	my_ift(x) = fftshift(ifft(ifftshift(x)))
end

# ╔═╡ c53325bd-eeb6-4663-867d-44fd3ace273b
md"## Task 1.1
Try to change the following test below.
You are only allowed to insert `fftshift` and `ifftshift` statements.
Don't change the order of the Fourier transform calls.
"

# ╔═╡ b57917e8-3d2e-433b-84df-6129f760f954
begin
	arr_even = [1.0,2,3,4]
	arr_odd = [1.0,2,3,4,5]
end

# ╔═╡ 6e0e223d-9107-4121-880a-8ed5d7e5e9c2
md"The Test is broken but we fix it by changing the right hand side. Do that for the red tests below as well"

# ╔═╡ 5213f82b-5e06-4cee-a0d1-21f1c5cc2998
PlutoTest.@test ffts(arr_even) ≈ fft(arr_even)

# ╔═╡ 64103e19-ce79-48b7-bac9-b429f6d423c2
# this one is fixed
PlutoTest.@test ffts(arr_even) ≈ fftshift(fft(arr_even))

# ╔═╡ f380db27-8773-413b-978f-4496fc585ae3
md"### Task 1.1
Now try to fix always the right hand side accordingly!
"

# ╔═╡ 0464d972-fcaf-4f83-be7d-073d217e8f4c
# TODO
PlutoTest.@test ift(arr_odd) ≈ fftshift(iffts(arr_odd))

# ╔═╡ 87fabbfe-e25b-467a-9fbd-ceec98ba4ed6
# TODO
PlutoTest.@test iffts(arr_odd) ≈ ifft(ifftshift(arr_odd))

# ╔═╡ c537e21b-5ed2-4cac-947a-e052474a2442
# TODO
PlutoTest.@test ffts(arr_odd) ≈ ft(fftshift(arr_odd))

# ╔═╡ 3093ae8d-27d4-4d32-8021-80fc6b0d1472
# TODO
PlutoTest.@test ifft(arr_odd) ≈ ifftshift(ift(fftshift(arr_odd)))

# ╔═╡ a43b5bb5-d56b-43c7-aa96-ed2d9f217373
md"# 2 Convolution

From calculus class we know that a convolution can be expressed via Fourier transforms like

$U * h = \mathcal{F}^{-1}\left[\mathcal{F}[U] \cdot \mathcal{F}[h] \right] = \mathcal{F}^{-1}\left[\mathcal{F}[U] \cdot H \right]$

where $*$ is a convolution and $\cdot$ an elementwise multiplication. $H$ is the OTF of the PSF $h$.

Now implement it yourself!
"

# ╔═╡ a0362131-41ad-4151-b851-d75af22791e1
"""
	my_conv(U, h)

Calculates a FFT based convolution between U and h.
The output is a real array!
So either use real valued transforms `rfft` or use `fft` with a final `real` conversion.

# Example
```julia-repl
julia> my_conv([1.0,0,0,0], [1.0, 0.5, 0.0, 0.5])
4-element Vector{ComplexF64}:
 1.0 + 0.0im
 0.5 + 0.0im
 0.0 + 0.0im
 0.5 + 0.0im
```
"""
function my_conv(U::AbstractArray{T, N}, h; center=ntuple(i -> 1, N)) where {T, N}
	# TODO
	return real(ifft(fft(U) .* fft(circshift(h, .- center .+ 1))))
end

# ╔═╡ b568e681-7eaf-4c1a-9f36-f46163c35041
begin
	img = Float32.(testimage("mandril_gray"))
	Gray.(img)
end;

# ╔═╡ cdf36a7b-54ff-4ab7-9497-81d1ed684be3
# some simple kernel
kernel = IndexFunArrays.normal(img, sigma=7);

# ╔═╡ 0d9e2e69-e51f-483c-aadd-6bde1ab0c28a
md"You should see the blurry monkey here. If not, it might be wrong."

# ╔═╡ 5436a802-9d78-48ff-af79-2cf8065fd514
Gray.(my_conv(img, kernel, center=(257, 257)))

# ╔═╡ 2eb0bea1-b6dc-417f-81b3-435fede2aa66
md"## 2 Test"

# ╔═╡ c7a4ce47-9bdf-46f3-88a4-888df19b9cb1
PlutoTest.@test my_conv([1.0,2,3,4], [1.0,0,0,0]) ≈ [1,2,3,4]

# ╔═╡ 3b479c44-4f6d-41ad-af70-04616a64154c
PlutoTest.@test my_conv([1.0,2,3,4], [0.0,1.0,0,0], center=2) ≈ [1,2,3,4]

# ╔═╡ 39807680-27cf-444d-9d93-d845baab61a4
PlutoTest.@test my_conv([1.0,2,3,4,5], [0.0,0.0,1.0,0,0], center=3) ≈ [1,2,3,4, 5]

# ╔═╡ 51bcf5d5-0405-467e-b528-275965fa1287
PlutoTest.@test my_conv(img, IndexFunArrays.delta(img, offset=(1,1))) ≈ img

# ╔═╡ da0d4ab1-b1d6-4817-ad38-6d3213c4e075
PlutoTest.@test my_conv([1, 2, 3, 4, 5, 6], [-1, 2, 2, 1, 3, 1]) ≈ [36.0, 32.0, 28.0, 30.0, 20.0, 22.0]

# ╔═╡ 09321d84-f9a7-412d-8fe4-3ece1bd90b21
md"# 3 Incoherent Image Formation

In this part we want to simulate an incoherent imaging situation (fluorescent microscope).
We simplify it by only considering an unitless `radius`.

As reminder, the qualitative procedure is as following:

We take a point source, Fourier transform it, take only the frequency inside the radius, go back to real space, take the absolute squared.

Pay attention where the frequencies should be located when you apply the `circ`.
"

# ╔═╡ ebfe7f3a-098a-41d8-a5bd-2c1e2b374fe0
begin
	"""
		circ(size, radius)

	`size` the size of the resulting array
	and `radius` the radius of the circle
	"""
	circ(size, radius) = rr2(size) .<= radius^2
end

# ╔═╡ 772fe3f9-6db8-4df4-8813-cbb334035351
"""
	calc_psf(arr_size, radius)

Calculate a simple PSF. `arr_size` is the output size and `radius` is the radius
in Fourier space.
The output array is normalized by its sum to 1.
"""
function calc_psf(arr_size, radius)
	# TODO
	rr2 = circ(arr_size, radius)
	h = abs2.(ift(rr2))
	h ./= sum(h)
	return h
end

# ╔═╡ 33b3a409-2690-4a4a-a091-cdfe6a831c59
md"r = $(@bind r PlutoUI.Slider(0.01:0.1:10.0, show_value=true))"

# ╔═╡ f209c5de-c2bd-4a0b-bbfd-6831e0254023
gray_show(calc_psf((64, 64), r))

# ╔═╡ 8b922c48-7d56-4e5b-b039-789e281c5fe1
md"r2 = $(@bind r2 PlutoUI.Slider(1:1:256, show_value=true))"

# ╔═╡ 119da3f0-0f1d-4a65-93cd-868f3c1d5f3e
h = calc_psf(size(img), r2);

# ╔═╡ 5830d782-67a3-4cae-847c-bdbdf0217aa7
# change this line such that the monkey is correctly centerd
# TODO
Gray.(my_conv(img, h, center=size(h) .÷ 2 .+ 1))

# ╔═╡ fbdcb172-b079-46a5-b8f3-f5ece30fe25a
md"## 3 Test"

# ╔═╡ 66db7400-9791-4952-9812-39b22829b29a
# large radius is a perfect optical system -> delta peak
PlutoTest.@test calc_psf((2, 2), 1000) ≈  [0 0; 0 1]

# ╔═╡ baff7f36-d18a-4434-b979-662b6d44eb46
PlutoTest.@test sum(calc_psf((13, 12), 3)) ≈ 1

# ╔═╡ 7d92a5c5-d48d-429b-85fa-63904f21fc62
PlutoTest.@test minimum(calc_psf((13, 12), 3)) ≥ 0

# ╔═╡ fd85fd4a-3e43-431c-aaa6-e92848c9e304
# compare to (approx) analytical solution
begin
	h2 = jinc.(7.25 / 9.219π * IndexFunArrays.rr((64, 64))).^2
	h2 ./= sum(h2)
	PlutoTest.@test ≈(1 .+ h2, 1 .+ calc_psf((64, 64), 7.25), rtol=0.001)
end

# ╔═╡ 08936039-7557-4ea3-8e26-5fbcdf12aec2
md"# 4 Generalized Wiener Filtering

A simpled deconvolution approach suited for Gaussian noise is the Wiener filter.
You can find the details in the slides.
Try to implement it here!
"

# ╔═╡ 7256f510-10ac-4fb2-99fc-0ffbcd45aae3
function wiener_filter(img, h, ϵ)
	# TODO
	otf = fft(h)
	wiener = conj.(otf) ./ (abs2.(otf) .+ ϵ)
	return abs.(real(ifft(fft(img) .* wiener)))
end

# ╔═╡ 9f4d5c69-c5a8-4a34-9100-e8209ede71b4
begin
	# PSF
	h3 = ifftshift(calc_psf(size(img), 20))
	gray_show(h3)
end

# ╔═╡ 2d5653cf-fa36-4a0e-99cd-4c9769b84705
begin
	img_b = my_conv(img, h3)
	Gray.(img_b)
end

# ╔═╡ 0f8555c7-a1c2-4643-912b-6816019a848a
img_gauss = add_gauss(img_b, 0.1);

# ╔═╡ 34ced1ed-82f4-4c32-9cec-086ff2a62bce
img_poisson = poisson(img_b, 20);

# ╔═╡ cbd3d6c4-bf40-4431-a7b9-5b26b411a9b9
Gray.([img_gauss img_poisson])

# ╔═╡ 29cc8e86-c109-411e-ba86-a8155f7c3a94
md"
pow1 = $(@bind pow1 Slider(-6:0.1:-0, show_value=true))

pow2 = $(@bind pow2 Slider(-6:0.1:-0, show_value=true))
"

# ╔═╡ 82542705-e24e-409e-a7bd-491ce750007e
img_gauss_wiener = wiener_filter(img_gauss, h3, 10^pow1);

# ╔═╡ e50febee-5c4b-44d7-9d78-29395a8c3ab6
img_poisson_wiener = wiener_filter(img_poisson, h3, 10^pow2);

# ╔═╡ 0b46da45-b743-40a1-a1f8-0581b6fd741a
Gray.([img_gauss_wiener  img_poisson_wiener])

# ╔═╡ f85f964d-0443-4e19-b655-036f82a0ba69
md"## 4 Test"

# ╔═╡ 94f64888-b732-4367-a905-b57de684fcf7
PlutoTest.@test wiener_filter([1.0, 2.0], [1.0, 0.0], 0) ≈ [1.0, 2.0]

# ╔═╡ 3d89c0be-a888-4e2f-8820-85e07bd6be30
PlutoTest.@test  wiener_filter([1.0, 2.0], [1.0, 0.1], 0)  ≈ [0.808081, 1.91919] rtol=1e-5

# ╔═╡ c1e8419c-c748-490b-93d3-9bbcde4f4da9
md"# 5 Gradient Descent Optimization
In this part we want to implement an optimization routine with a _strange_ sensor.

The sensor has an additive Gaussian noise part but also an cubic gain behaviour.
See the function below
"

# ╔═╡ c337e8cf-bab9-46cd-aae8-22f6265f9cb7
begin
	function strange_sensor_f(value::T) where T
		value .^2
	end

	function strange_sensor(value::T) where T
		return abs.(randn(T) * 0.10f0 + strange_sensor_f(value))
	end
end

# ╔═╡ ff2288f0-2edf-4f22-90ef-af5777676ae7
# strange output
strange_sensor.([1.0 2; 3 4])

# ╔═╡ d23fe57a-7255-4d8f-a00a-eaec953213b2
md"
First we simulate the full `img` of the mandril with that sensor.
Clearly the appearance has changed due to the quadratic behaviour but also noise is visible.
"

# ╔═╡ a0f77be2-c40c-4372-a245-45b76ddc5861
begin
	img_strange = strange_sensor.(img) # todo
	gray_show(img_strange)
end

# ╔═╡ cdc25e03-c6fd-4ee0-801e-24654a83d65d
md"## 5.1 Task
Since our sensor is quite noisy, we want to measure the image multiple times.
Try to complete `measure`.
"

# ╔═╡ 4af883c2-fdc9-458c-998b-ff0b8b5df146
"""
	measure(img, N)

Measure the `img` `N` times. Return a vector with the `N` measured images.
"""
function measure(img, N)
	# TODO
	return [strange_sensor.(img) for i = 1:N]
end

# ╔═╡ c02969bc-3b5e-4422-a3c9-3a7d00cb3ddf
imgs = measure(img, 10); # TODO: simulate 10 times

# ╔═╡ e3a8ef81-0d19-4e75-b137-6454ce262991
Gray.([reduce(hcat, imgs[1:end÷2]); reduce(hcat, imgs[end÷2+1:end])])

# ╔═╡ ebd74bcb-4187-469f-a206-ecf9041918f1
md"## 5.1 Test"

# ╔═╡ 7579c331-fa53-47f4-a479-312f2f7a3931
PlutoTest.@test typeof(measure([1.0], 12)) <: Vector

# ╔═╡ 3cca079c-1245-493b-ae6c-9db18e11835c
PlutoTest.@test length(measure([1.0], 12)) == 12

# ╔═╡ 0ae38ccc-9129-4c61-9810-dbc448e8a07f
begin
	Random.seed!(42)
	a = measure([1.0 2.0; 3 4], 3)
end;

# ╔═╡ 4f5a9a27-1893-49e4-a2e9-59ba3a2340a2
begin
		Random.seed!(42)
		b = [abs.([1.0 2.0; 3 4].^2 .+ randn((2,2)) * 0.1f0) for i = 1:3]
end;

# ╔═╡ 37ce2f96-377d-48d7-bc56-f83b0cce349c
PlutoTest.@test a ≈ b

# ╔═╡ cd0829a1-cdf1-45c9-a4dd-91bb3ec5bb03
md"## 5.2 Loss Function
Having our 10 images we would like to retrieve a best guess for the underlying image.
Taking the mean does not work in this case since the input image is modified by `strange_sensor_f`.

Therefore, we interprete the reconstruction as an optimization problem.
Try to find the corresponding pages in the lecture.

Generally, we want to minimize a loss function which looks like

$\underset{\mu}{\mathrm{argmin}}\, \mathcal L = \sum_{\text{img} \, \in \, \text{imgs}} \sum_{p_i \, \in \, \text{img}} (\text{img}[i] - f(\mu[i]))^2$

So we sum over all measured images. For each image, we additionally sum each pixel and  compare it with $f(\mu[i])$.
$f$ is the same as `strange_sensor_f` and by calling $f(\mu[i])$ we apply $f$ to the reconstruction $\mu$. Via that function call, we hope that we find a $\mu$ which fits to the measurments. So $f$ is the forward model of the sensor (without the noise part).
"

# ╔═╡ 328b5994-b1c4-4850-b8d3-2f3781bed99f
"""
	loss(imgs, μ)

Calculate the loss between the `ìmgs` and `μ`.
Basically implement the sum of the square difference value.
Don't forget to apply `strange_sensor_f` to `μ` before!

Using two for loops is perfectly fine!
"""
function loss(imgs::Vector{<:AbstractArray{T, N}}, μ::AbstractArray{T, N}) where {T, N}
	# TODO
	value = zero(T)
	for img in imgs
		value += sum( (img .- strange_sensor_f.(μ)) .^2)
	end
	return value
end

# ╔═╡ bc225836-b68f-40c4-8e12-fa9e84ffac44
loss

# ╔═╡ 78da1c47-1de7-4d95-a1e1-cd60e912a0fc
# comparison with ground truth image
loss(imgs, img)

# ╔═╡ 3d182fa9-84bb-4504-8945-bf653ce1f99d
# the ground truth is not the minimum...
# why not?
loss(imgs, img .+0.0001f0)

# ╔═╡ b56a311e-912a-4e7c-b821-4124b847194a
md"## 5.2 Test"

# ╔═╡ 840ddb7e-56ba-4d7a-8f08-cee67128315f
PlutoTest.@test loss([[1]], [2]) isa Number

# ╔═╡ 528194f2-9f74-4dda-94f8-3e4f0e515bc9
PlutoTest.@test loss([[1]], [2])  ≈ 9

# ╔═╡ a203fb4f-c3d6-4d46-9c58-2454b71b0e1b
PlutoTest.@test loss([[1], [2], [4.0]], sqrt.([2]))  ≈ 5

# ╔═╡ 063ef53f-d152-4282-af5f-7f2addce3ab0
md"## 5.3 Gradient of `strange_sensor_f`
For the optimization we later want to apply a gradient descent optimization scheme.
Hence, we need a few gradients!
"

# ╔═╡ 95392d70-07cf-47a7-8707-f1dab177e7a5
"""
	 grad_strange_sensor_f(value::T)

Calculate the gradient of `strange_sensor_f`.
Use the rules you know already from school!
"""
function grad_strange_sensor_f(value::T) where T
	2 .* value # TODO
end

# ╔═╡ 01de98fe-5537-4745-acf8-99daa0b8aa6f
grad_strange_sensor_f(1)

# ╔═╡ 75549871-14a3-4b9f-9b98-cccc34ed315d
md"## 5.3 Test"

# ╔═╡ 867b3bf0-1208-48a1-b1ee-904d11b28e1f
PlutoTest.@test grad_strange_sensor_f(0) ≈ 0

# ╔═╡ 7c2d2d6d-f933-452e-9870-7dce7ad4bf1d
# comparison with automatic differentation package!
PlutoTest.@test ForwardDiff.derivative(strange_sensor_f, 42) ≈ grad_strange_sensor_f(42)

# ╔═╡ cd363a92-a473-47f8-b41b-c5d5249fec90
md"## 5.4 Gradient of Loss
Now we come to the last part of the gradient.
We need the full gradient of the loss function.

The gradient of the loss function with respect to `μ` will be an array again.
That is plausible since the loss function accounts for all pixels. By changing a single pixel we also change the value of the loss function. Therefore, for each pixel a gradient exists describing the influence to the loss value.

You need to derive the loss with respect to $\mu$. But don't forget about the chain rule for `strange_sensor_f`!

Again, two for loops are perfectly fine!
"

# ╔═╡ 1ad4cf95-a3fb-4e45-b7bd-80ab724c158e
function gradient(imgs::Vector{<:AbstractArray{T, N}}, μ::AbstractArray{T, N}) where {T, N}
	# TODO
	∇ = zero.(μ)
	for img in imgs
		∇ .+= 2 .* (img .- strange_sensor_f.(μ)) .* (.- grad_strange_sensor_f.(μ))
	end
	return ∇
end

# ╔═╡ cb41a2f9-a0c5-4eeb-af51-9bf40a5d7fd6
gradient(imgs, img)

# ╔═╡ f72450cb-a54e-4827-9c33-7a44f561bd43
md"## 5.4 Test"

# ╔═╡ 7e6eb33d-0721-4bb6-a9ee-1cfb07a17e46
PlutoTest.@test gradient([[1.0], [2.0]], [2.5]) ≈ [95]

# ╔═╡ cf550002-13ce-4788-8206-790fb355a91b
PlutoTest.@test gradient([[3.5, 2.1], [2.0, 0.0]], [3.23, 23.2]) isa Vector

# ╔═╡ bfd35c9f-1449-4d86-bb50-9bbbd6ea2c30
PlutoTest.@test size(gradient([[3.5, 2.1], [2.0, 0.0]], [3.23, 23.2])) == (2,)

# ╔═╡ c1795dd9-f8a5-472d-bd82-7871e8603534
PlutoTest.@test  gradient([[3.5, 2.1], [2.0, 0.0]], [3.23, 23.2]) ≈ [198.526136, 99702.46399999999]

# ╔═╡ c5436e99-513b-4c51-b3ad-f82318304a3e
md"## 5.5 Gradient Descent"

# ╔═╡ b393a0f5-af03-4ded-8e22-886022f5da30
"""
	gradient_descent(imgs, μ, N_iter, step_size)

Runs a gradient descent using `gradient` and the experimental images `imgs`.
`N_iter` is the number of iterations, `step_size` is the step size of the gradient step.

The optimized output is `μ_optimized`.
"""
function gradient_descent(imgs, μ, N_iter, step_size)
	# TODO
	μ_optimized = copy(μ)
	for i = 1:N_iter
		μ_optimized .= μ_optimized .- step_size .* gradient(imgs, μ_optimized)
	end
	return μ_optimized
end

# ╔═╡ c81c5eaa-0e0f-4730-9148-ec0e7b63cdd4
μ_init = one.(imgs[1]); # TODO, what could be a good initilization?

# ╔═╡ 32f6d83a-d5db-4a21-b2c4-d8876de38c46
md"
Try to change the gradient step and the number of iterations.

Number of iterations $(@bind N_iter Slider(1:100, show_value=true))

step size $(@bind step_size Slider(1f-5:1f-5:1f-1, show_value=true))
"

# ╔═╡ f7d7cf86-ece9-406b-af78-f86177a767c4
μ_optimized = gradient_descent(imgs, μ_init, N_iter, step_size);

# ╔═╡ 2110469f-16f9-4454-996a-a949a94fffa3
# that value should get smaller with more iterations
# smaller is better
loss(imgs, μ_optimized)

# ╔═╡ 3643e55e-3205-4561-b271-058d59f46685
# this value should be larger than ↑
loss(imgs, μ_init)

# ╔═╡ 85c81f80-c959-4682-932b-26fe7f415f4d
Gray.(μ_optimized)

# ╔═╡ 09313c1d-102f-4e0b-b9c8-7d2e904a1dd6
md"## 5.5 Test"

# ╔═╡ b2e7aae8-2777-4657-ba3d-f0dfcb581ede
PlutoTest.@test gradient_descent([[1.0], [2.0]], [3.0], 3, 0.01) ≈ [1.21021] rtol=1f-4

# ╔═╡ ae027113-998b-467e-8967-cc416d480377
PlutoTest.@test loss(imgs, μ_init) > loss(imgs, μ_optimized)

# ╔═╡ 4ff0afae-af4c-45d0-83b2-4fb5eead2fa1
PlutoTest.@test loss(imgs, gradient_descent(imgs, (imgs[1] .+ imgs[2]) / 2, 15, 0.000001);) < loss(imgs, imgs[1])

# ╔═╡ 2374bfc1-cd14-4bee-b5e5-bb15e0e6e253
begin
	μ_mean = sqrt.(reduce((a,b) -> a .+ b, imgs) ./ 10)
	Gray.(μ_mean)
end;

# ╔═╡ 9360c9a0-6a0f-4ff6-a215-23ce5c3c8099
md"#### Can you beat the _mean_?"

# ╔═╡ a837b22c-61f1-4a98-a362-c08b405c4dca
PlutoTest.@test loss(imgs, μ_optimized) < loss(imgs, μ_mean)

# ╔═╡ e2aec734-ce2b-4761-b3d0-b559df8b17da
md"Note that in this particular example the mean can be proven to be the best estimator, so it is not surprising to not be able to beat the mean. Yet, if the model cannot be inverted, we still may want to use a gradient-based optimization to solve the problem."

# ╔═╡ 959b4745-7742-459a-bef0-e374ec4aec17
[Gray.(img) Gray.(μ_mean) Gray.(μ_optimized)]

# ╔═╡ b3d54170-466c-4b0d-a337-280ef4ea87f3
md"## 5.6 Test
To check wether your gradient and your loss is correct, see also the following output of Optim (a sophisticated package for optimization).
If the output is not a nice Mandrill, there is most likely something wrong with your loss/gradients.
"

# ╔═╡ eafe5f6e-18ab-4341-8435-a3c0b9423b35
begin
	g_optim!(G, x) = isnothing(G) ? nothing : G .= gradient(imgs, x)
	f_optim(x) = loss(imgs, x)
end

# ╔═╡ 13304ea6-0639-4863-827e-d66a8630bc63
begin
	μ_optim_init = copy(μ_init)
	res = optimize(f_optim, g_optim!, μ_optim_init, ConjugateGradient(), Optim.Options(iterations=50))
	μ_optim = Optim.minimizer(res)
	res
end

# ╔═╡ c685ac28-79e8-4eda-8c0f-0944a1276691
Gray.(Optim.minimizer(res))

# ╔═╡ 7a691c80-e08d-423a-a69a-572f02fdddda
loss(imgs, μ_optim)

# ╔═╡ 19a0e512-46d6-429e-93f0-05e2faf95218
md"## Can you beat Optim?"

# ╔═╡ abe5aca2-61ce-4f99-ad75-0c84293cd7b3
PlutoTest.@test loss(imgs, μ_optimized) < loss(imgs, μ_optim)

# ╔═╡ 1171004b-96a1-4ee4-bbd5-2b7f2b1d582a
# yes, but hard with an hand optimized iterations/value pair

# ╔═╡ Cell order:
# ╠═20bb3256-f1bb-4b3b-b410-89969443fd4d
# ╠═db61ecb1-6ed9-44d9-b8e4-b50622efbac1
# ╟─82726a23-060c-4c7d-91c1-e78a22fb98d6
# ╠═ebe90686-c72f-4a42-a644-f06babe00f41
# ╠═5053c100-5f27-4f2f-acbe-e12c9070feac
# ╠═361a7a34-d70b-4597-948c-80ed913a562b
# ╟─653a94c9-cb89-47df-8ba5-0ce683440b78
# ╠═b65e0745-82c6-48e1-a956-fa71583dad15
# ╟─c53325bd-eeb6-4663-867d-44fd3ace273b
# ╠═b57917e8-3d2e-433b-84df-6129f760f954
# ╟─6e0e223d-9107-4121-880a-8ed5d7e5e9c2
# ╠═5213f82b-5e06-4cee-a0d1-21f1c5cc2998
# ╠═64103e19-ce79-48b7-bac9-b429f6d423c2
# ╟─f380db27-8773-413b-978f-4496fc585ae3
# ╠═0464d972-fcaf-4f83-be7d-073d217e8f4c
# ╠═87fabbfe-e25b-467a-9fbd-ceec98ba4ed6
# ╠═c537e21b-5ed2-4cac-947a-e052474a2442
# ╠═3093ae8d-27d4-4d32-8021-80fc6b0d1472
# ╟─a43b5bb5-d56b-43c7-aa96-ed2d9f217373
# ╠═a0362131-41ad-4151-b851-d75af22791e1
# ╠═b568e681-7eaf-4c1a-9f36-f46163c35041
# ╠═cdf36a7b-54ff-4ab7-9497-81d1ed684be3
# ╟─0d9e2e69-e51f-483c-aadd-6bde1ab0c28a
# ╠═5436a802-9d78-48ff-af79-2cf8065fd514
# ╟─2eb0bea1-b6dc-417f-81b3-435fede2aa66
# ╠═c7a4ce47-9bdf-46f3-88a4-888df19b9cb1
# ╠═3b479c44-4f6d-41ad-af70-04616a64154c
# ╠═39807680-27cf-444d-9d93-d845baab61a4
# ╠═51bcf5d5-0405-467e-b528-275965fa1287
# ╠═da0d4ab1-b1d6-4817-ad38-6d3213c4e075
# ╟─09321d84-f9a7-412d-8fe4-3ece1bd90b21
# ╠═ebfe7f3a-098a-41d8-a5bd-2c1e2b374fe0
# ╠═772fe3f9-6db8-4df4-8813-cbb334035351
# ╟─33b3a409-2690-4a4a-a091-cdfe6a831c59
# ╠═f209c5de-c2bd-4a0b-bbfd-6831e0254023
# ╟─8b922c48-7d56-4e5b-b039-789e281c5fe1
# ╠═119da3f0-0f1d-4a65-93cd-868f3c1d5f3e
# ╠═5830d782-67a3-4cae-847c-bdbdf0217aa7
# ╟─fbdcb172-b079-46a5-b8f3-f5ece30fe25a
# ╠═66db7400-9791-4952-9812-39b22829b29a
# ╠═baff7f36-d18a-4434-b979-662b6d44eb46
# ╠═7d92a5c5-d48d-429b-85fa-63904f21fc62
# ╠═fd85fd4a-3e43-431c-aaa6-e92848c9e304
# ╟─08936039-7557-4ea3-8e26-5fbcdf12aec2
# ╠═7256f510-10ac-4fb2-99fc-0ffbcd45aae3
# ╠═9f4d5c69-c5a8-4a34-9100-e8209ede71b4
# ╠═2d5653cf-fa36-4a0e-99cd-4c9769b84705
# ╠═0f8555c7-a1c2-4643-912b-6816019a848a
# ╠═34ced1ed-82f4-4c32-9cec-086ff2a62bce
# ╠═cbd3d6c4-bf40-4431-a7b9-5b26b411a9b9
# ╠═29cc8e86-c109-411e-ba86-a8155f7c3a94
# ╠═82542705-e24e-409e-a7bd-491ce750007e
# ╠═e50febee-5c4b-44d7-9d78-29395a8c3ab6
# ╠═0b46da45-b743-40a1-a1f8-0581b6fd741a
# ╟─f85f964d-0443-4e19-b655-036f82a0ba69
# ╠═94f64888-b732-4367-a905-b57de684fcf7
# ╠═3d89c0be-a888-4e2f-8820-85e07bd6be30
# ╟─c1e8419c-c748-490b-93d3-9bbcde4f4da9
# ╠═c337e8cf-bab9-46cd-aae8-22f6265f9cb7
# ╠═ff2288f0-2edf-4f22-90ef-af5777676ae7
# ╟─d23fe57a-7255-4d8f-a00a-eaec953213b2
# ╠═a0f77be2-c40c-4372-a245-45b76ddc5861
# ╟─cdc25e03-c6fd-4ee0-801e-24654a83d65d
# ╠═4af883c2-fdc9-458c-998b-ff0b8b5df146
# ╠═c02969bc-3b5e-4422-a3c9-3a7d00cb3ddf
# ╠═e3a8ef81-0d19-4e75-b137-6454ce262991
# ╟─ebd74bcb-4187-469f-a206-ecf9041918f1
# ╠═7579c331-fa53-47f4-a479-312f2f7a3931
# ╠═3cca079c-1245-493b-ae6c-9db18e11835c
# ╟─0ae38ccc-9129-4c61-9810-dbc448e8a07f
# ╟─4f5a9a27-1893-49e4-a2e9-59ba3a2340a2
# ╠═37ce2f96-377d-48d7-bc56-f83b0cce349c
# ╟─cd0829a1-cdf1-45c9-a4dd-91bb3ec5bb03
# ╠═328b5994-b1c4-4850-b8d3-2f3781bed99f
# ╠═bc225836-b68f-40c4-8e12-fa9e84ffac44
# ╠═78da1c47-1de7-4d95-a1e1-cd60e912a0fc
# ╠═3d182fa9-84bb-4504-8945-bf653ce1f99d
# ╟─b56a311e-912a-4e7c-b821-4124b847194a
# ╠═840ddb7e-56ba-4d7a-8f08-cee67128315f
# ╠═528194f2-9f74-4dda-94f8-3e4f0e515bc9
# ╠═a203fb4f-c3d6-4d46-9c58-2454b71b0e1b
# ╟─063ef53f-d152-4282-af5f-7f2addce3ab0
# ╠═95392d70-07cf-47a7-8707-f1dab177e7a5
# ╠═01de98fe-5537-4745-acf8-99daa0b8aa6f
# ╟─75549871-14a3-4b9f-9b98-cccc34ed315d
# ╠═867b3bf0-1208-48a1-b1ee-904d11b28e1f
# ╠═7c2d2d6d-f933-452e-9870-7dce7ad4bf1d
# ╟─cd363a92-a473-47f8-b41b-c5d5249fec90
# ╠═1ad4cf95-a3fb-4e45-b7bd-80ab724c158e
# ╠═cb41a2f9-a0c5-4eeb-af51-9bf40a5d7fd6
# ╟─f72450cb-a54e-4827-9c33-7a44f561bd43
# ╠═7e6eb33d-0721-4bb6-a9ee-1cfb07a17e46
# ╠═cf550002-13ce-4788-8206-790fb355a91b
# ╠═bfd35c9f-1449-4d86-bb50-9bbbd6ea2c30
# ╠═c1795dd9-f8a5-472d-bd82-7871e8603534
# ╟─c5436e99-513b-4c51-b3ad-f82318304a3e
# ╠═b393a0f5-af03-4ded-8e22-886022f5da30
# ╠═c81c5eaa-0e0f-4730-9148-ec0e7b63cdd4
# ╟─32f6d83a-d5db-4a21-b2c4-d8876de38c46
# ╠═f7d7cf86-ece9-406b-af78-f86177a767c4
# ╠═2110469f-16f9-4454-996a-a949a94fffa3
# ╠═3643e55e-3205-4561-b271-058d59f46685
# ╠═85c81f80-c959-4682-932b-26fe7f415f4d
# ╟─09313c1d-102f-4e0b-b9c8-7d2e904a1dd6
# ╠═b2e7aae8-2777-4657-ba3d-f0dfcb581ede
# ╠═ae027113-998b-467e-8967-cc416d480377
# ╠═4ff0afae-af4c-45d0-83b2-4fb5eead2fa1
# ╠═2374bfc1-cd14-4bee-b5e5-bb15e0e6e253
# ╟─9360c9a0-6a0f-4ff6-a215-23ce5c3c8099
# ╠═a837b22c-61f1-4a98-a362-c08b405c4dca
# ╟─e2aec734-ce2b-4761-b3d0-b559df8b17da
# ╠═959b4745-7742-459a-bef0-e374ec4aec17
# ╟─b3d54170-466c-4b0d-a337-280ef4ea87f3
# ╠═eafe5f6e-18ab-4341-8435-a3c0b9423b35
# ╠═13304ea6-0639-4863-827e-d66a8630bc63
# ╠═c685ac28-79e8-4eda-8c0f-0944a1276691
# ╠═7a691c80-e08d-423a-a69a-572f02fdddda
# ╟─19a0e512-46d6-429e-93f0-05e2faf95218
# ╠═abe5aca2-61ce-4f99-ad75-0c84293cd7b3
# ╟─1171004b-96a1-4ee4-bbd5-2b7f2b1d582a


text:
using DataFrames
using Dates

@testset "simulate" begin
    @testset "basic" begin
        @system SSimulate(Controller) begin
            a => 1 ~ preserve(parameter)
            b(a) ~ accumulate
        end
        n = 10
        r = simulate(SSimulate, stop=n)
        @test r isa DataFrame
        @test size(r, 1) == (n+1)
        @test propertynames(r) == [:time, :a, :b]
        @test r[end, :time] == n*u"hr"
        @test r[end, :a] == 1
        @test r[end, :b] == n
        r = simulate(SSimulate, stop=n, config=(:SSimulate => :a => 2))
        @test r[end, :a] == 2
        @test r[end, :b] == 2n
        r = simulate(SSimulate, stop=n, target=[:b])
        @test size(r, 2) == 2
        @test propertynames(r) == [:time, :b]
        r = simulate(SSimulate, stop=n, index=:b, target=[:b])
        @test size(r, 2) == 1
        @test propertynames(r) == [:b]
    end

    @testset "stop number" begin
        @system SSimulateStopNumber(Controller) begin
            a => 1 ~ preserve(parameter)
            b(a) ~ accumulate
        end
        r = simulate(SSimulateStopNumber, stop=5)
        @test r[end-1, :b] == 4
        @test r[end, :b] == 5
    end

    @testset "stop number unit" begin
        @system SSimulateStopNumberUnit(Controller) begin
            a => 1 ~ preserve(parameter)
            b(a) ~ accumulate
        end
        r1 = simulate(SSimulateStopNumberUnit, stop=2u"hr")
        @test r1[end, :b] == 2 * 1
        r2 = simulate(SSimulateStopNumberUnit, stop=2u"d")
        @test r2[end, :b] == 2 * 24
        r3 = simulate(SSimulateStopNumberUnit, stop=2u"yr")
        @test r3[end, :b] == 2 * 24 * 365.25
    end

    @testset "stop count" begin
        @system SSimulateStopCount(Controller) begin
            a => 1 ~ preserve(parameter)
            b(a) ~ accumulate
            c => 5 ~ preserve::int
        end
        r = simulate(SSimulateStopCount, stop=:c)
        @test r[end-1, :b] == 4
        @test r[end, :b] == 5
    end

    @testset "stop count unit" begin
        @system SSimulateStopCountUnit(Controller) begin
            a => 1 ~ preserve(parameter)
            b(a) ~ accumulate
            c ~ preserve::Cropbox.Quantity(parameter)
        end
        r1 = simulate(SSimulateStopCountUnit, stop=:c, config=:0 => :c => 2u"hr")
        @test r1[end, :b] == 2 * 1
        r2 = simulate(SSimulateStopCountUnit, stop=:c, config=:0 => :c => 2u"d")
        @test r2[end, :b] == 2 * 24
        r3 = simulate(SSimulateStopCountUnit, stop=:c, config=:0 => :c => 2u"yr")
        @test r3[end, :b] == 2 * 24 * 365.25
    end

    @testset "stop boolean" begin
        @system SSimulateStopBoolean(Controller) begin
            a => 1 ~ preserve(parameter)
            b(a) ~ accumulate
            z(b) => b >= 5 ~ flag
        end
        r1 = simulate(SSimulateStopBoolean, stop=:z)
        @test r1[end-1, :b] == 4
        @test r1[end, :b] == 5
        r2 = simulate(SSimulateStopBoolean, stop="z")
        @test r1 == r2
        r3 = simulate(SSimulateStopBoolean, stop=s -> s.b' >= 5)
        @test r1 == r3
    end

    @testset "snap" begin
        @system SSimulateSnap(Controller) begin
            a => 1 ~ preserve(parameter)
            b(a) ~ accumulate
            f(b) => (b % 2 == 0) ~ flag
        end
        n = 10
        f(s) = s.b' % 2 == 0
        r0 = simulate(SSimulateSnap, stop=n, snap=nothing)
        @test !all(r0.b .% 2 .== 0)
        r1 = simulate(SSimulateSnap, stop=n, snap=f)
        @test all(r1.b .% 2 .== 0)
        r2 = simulate(SSimulateSnap, stop=n, snap=:f)
        @test all(r2.b .% 2 .== 0)
        r3 = simulate(SSimulateSnap, stop=n, snap="f")
        @test all(r3.b .% 2 .== 0)
    end

    @testset "snatch" begin
        @system SSimulateSnatch(Controller) begin
            a => 1 ~ preserve(parameter)
            b(a) ~ accumulate
        end
        n = 3
        i = 0
        f(D, s) = begin
            d = D[1]
            @test s.a' == d[:a] && s.b' == d[:b]
            i += 1
            d[:c] = i
        end
        r = simulate(SSimulateSnatch, stop=n, snatch=f)
        @test i == 1 + n
        @test r[!, :c] == [1, 2, 3, 4]
    end

    @testset "callback" begin
        @system SSimulateCallback(Controller) begin
            a => 1 ~ preserve(parameter)
            b(a) ~ accumulate
        end
        n = 3
        i = 0
        f(s, m) = begin
            r = m.result[end, :]
            @test s.a' == r.a && s.b' == r.b
            i += 1
        end
        simulate(SSimulateCallback, stop=n, callback=f)
        @test i == n
    end

    @testset "layout" begin
        @system SSimulateLayout(Controller) begin
            i => 1 ~ accumulate
            a(i) => i-1 ~ track
            b(i) => 2i ~ track
        end
        L = [
            (target=:a,),
            (index=[:t => "context.clock.time", "i"], target=["a", :B => :b]),
            (base="context.clock", index="time", target="step"),
            (target=:b, meta=(:c => 0, :d => :D)),
        ]
        n = 1
        r = simulate(SSimulateLayout, L, stop=n)
        @test propertynames(r[1]) == [:time, :a]
        @test propertynames(r[2]) == [:t, :i, :a, :B]
        @test propertynames(r[3]) == [:time, :step]
        @test propertynames(r[4]) == [:time, :b, :c, :d]
        @test r[1][end, :time] == r[2][end, :t] == r[3][end, :time] == r[4][end, :time] == n*u"hr"
        @test r[1][end, :a] == 0
        @test r[2][end, :B] == 2
        @test r[3][end, :step] == 1u"hr"
        @test all(r[4][!, :c] .== 0) && all(r[4][!, :d] .== :D)
    end

    @testset "layout and configs" begin
        @system SSimulateLayoutConfigs(Controller) begin
            p ~ preserve(parameter)
            i => 1 ~ accumulate
            a(i, p) => p*(i-1) ~ track
            b(i, p) => 2p*i ~ track
        end
        L = [
            (index=:i, target=:a),
            (index=:t => "context.clock.time", target=:b),
            (target=[:i, :a, :b],),
        ]
        p1, p2 = 1, 2
        C = [
            :SSimulateLayoutConfigs => :p => p1,
            :SSimulateLayoutConfigs => :p => p2,
        ]
        n = 10
        r = simulate(SSimulateLayoutConfigs, L, C, stop=n)
        @test length(r) == length(L)
        o = r[3]
        @test o[o.time .== n*u"hr", :i] == [n, n]
        @test o[o.time .== n*u"hr", :a] == [p1*(n-1), p2*(n-1)]
        @test o[o.time .== n*u"hr", :b] == [2p1*n, 2p2*n]
    end

    @testset "configs" begin
        @system SSimulateConfigs(Controller) begin
            a ~ preserve(parameter)
            b(a) ~ accumulate
        end
        p = [1, 2]
        C = @config !(:SSimulateConfigs => :a => p)
        n = 10
        r = simulate(SSimulateConfigs, configs=C, stop=n)
        @test r[r.time .== n*u"hr", :a] == p
        @test r[r.time .== n*u"hr", :b] == p .* n
    end

    @testset "meta" begin
        @system SSimulateMeta(Controller) begin
            a ~ preserve(parameter)
            b(a) ~ accumulate
        end
        C = (
            :SSimulateMeta => (a=1,),
            :Extra => (b=2, c=:C),
        )
        n = 10
        r1 = simulate(SSimulateMeta, config=C, index=(), meta=:Extra, stop=n)
        @test propertynames(r1) == [:a, :b, :c]
        @test all(r1.b .== 2)
        @test all(r1.c .== :C)
        r2 = simulate(SSimulateMeta, config=C, index=(), meta=(:Extra, :d => 0), stop=n)
        @test propertynames(r2) == [:a, :b, :c, :d]
        @test r1.b == r2.b
        @test r1.c == r2.c
        @test all(r2.d .== 0)
    end

    @testset "wildcard" begin
        @system SSimulateWildcardA begin
            x => 1 ~ preserve
            y => 2 ~ preserve
        end
        @system SSimulateWildcard(Controller) begin
            A(context) ~ ::SSimulateWildcardA
            a => 3 ~ preserve
            b => 4 ~ preserve
        end
        r1 = simulate(SSimulateWildcard; target="*")
        @test names(r1) == ["time", "a", "b"]
        @test collect(r1[1,:]) == [0u"hr", 3, 4]
        r2 = simulate(SSimulateWildcard; target="A.*")
        @test names(r2) == ["time", "A.x", "A.y"]
        @test collect(r2[1,:]) == [0u"hr", 1, 2]
        r3 = simulate(SSimulateWildcard; target=["*", "A.*"])
        @test names(r3) == ["time", "a", "b", "A.x", "A.y"]
        @test collect(r3[1,:]) == [0u"hr", 3, 4, 1, 2]
        r4 = simulate(SSimulateWildcard; target=["A.*", "*"])
        @test names(r4) == ["time", "A.x", "A.y", "a", "b"]
        @test collect(r4[1,:]) == [0u"hr", 1, 2, 3, 4]
    end

    @testset "options" begin
        @system SSimulateOptions(Controller) begin
            a ~ preserve(extern)
            b ~ ::int(override)
        end
        n = 10
        a, b = 1, 2
        o = (; a, b)
        r = simulate(SSimulateOptions, options=o, stop=n)
        @test r[end, :a] == a
        @test r[end, :b] == b
    end

    @testset "seed" begin
        @system SSimulateSeed(Controller) begin
            a => rand() ~ track
            b(a) ~ accumulate
        end
        n = 10
        r1 = simulate(SSimulateSeed, seed=0, stop=n)
        @test r1[end, :a] == 0.5392892841426182
        @test r1[end, :b] == 3.766035118243237
        r2 = simulate(SSimulateSeed, seed=0, stop=n)
        @test r1 == r2
    end

    @testset "no seed" begin
        @system SSimulateNoSeed(Controller) begin
            a => rand() ~ track
            b(a) ~ accumulate
        end
        n = 10
        r1 = simulate(SSimulateNoSeed, seed=nothing, stop=n)
        r2 = simulate(SSimulateNoSeed, seed=nothing, stop=n)
        @test r1 != r2
    end

    @testset "extractable" begin
        @system SSimulateExtractable(Controller) begin
            a => 1 ~ track
            b => :hello ~ track::sym
            c => "world" ~ track::str
            d => DateTime(2020, 3, 1) ~ track::DateTime
            e => Dict(:k => 0) ~ track::Dict
            f => [1, 2, 3] ~ track::Vector
            g => (1, 2) ~ track::Tuple
        end
        r = simulate(SSimulateExtractable)
        N = propertynames(r)
        # default = Union{Number,Symbol,AbstractString,AbstractDateTime}
        @test :a ∈ N
        @test :b ∈ N
        @test :c ∈ N
        @test :d ∈ N
        @test :e ∉ N
        @test :f ∉ N
        @test :g ∉ N
    end
end

ext: jl


text:
struct Block
    composition::Compose.Context
    connections::Union{Dict{Symbol, XYTuple}, Vector{XYTuple}}
end

function Block(x::Real, y::Real, tag::Symbol, inputs=[:in1, :in2]; kwargs...)
    logic_tags = [:NOT, :AND, :OR, :XOR, :NAND, :NOR, :XNOR]
    if tag in logic_tags
        return Gate(x, y, tag, inputs, kwargs...)
    end
end

function and_gate(x, y, inputs, output=:out, sx=1.0, sy=1.0)
    p = (x, y)
    in_pos = partition(
        p + 0.5 * (-1sx, -1sy),
        p + 0.5 * (-1sx, 1sy),
        length(inputs)
    )
    connections = Dict{Symbol, XYTuple}(
        map((label, pos) -> Pair(label, pos), inputs, in_pos)...,
        :out => p + 0.5 * (1sx, 0sy)
    )
    composition = compose(
        context(
            x-0.55sx, y-0.55sy, 1.1sx, 1.1sy,
            units = UnitBox(-1.1, -1.1, 2.2, 2.2)
        ),
        and_base()
    )
    return Block(composition, connections)
end

function IC(
        x::Real, y::Real, name::String, pins::Matrix{Symbol};
        scale_x::Real = 1.0, scale_y::Real = 1.0, text_offset::Float64 = 0.05,
        _fontsize::Measures.Length = 20pt, _font::String = "Helvetica-Bold"
    )
    M, N = size(pins)
    p = (x, y)
    sx, sy = scale_x, scale_y

    rot = Rotation(pi/2)
    if M > N
        # vertical layout
        @assert (N == 2) "pins should be ($M, 2)"
        width, height = 1.5sx, M * 0.5sy
        left = partition(p + (-0.5width, -0.5height), p + (-0.5width, 0.5height), M)
        right = partition(p + (0.5width, -0.5height), p + (0.5width, 0.5height), M)
        connections = Dict{Symbol, XYTuple}(map(Pair, pins[:], [left..., right...][:]))

        text_comp = compose(
            context(),
            map(pins[:, 1], left) do s, tp
                text(tp[1] + text_offset, tp[2], string(s), hleft, vcenter)
            end...,
            map(pins[:, 2], right) do s, tp
                text(tp[1] - text_offset, tp[2], string(s), hright, vcenter)
            end...,
            linewidth(0mm), fontsize(0.75_fontsize), font(_font), fill("black")
        )
    else
        # horizontal layout
        @assert (M == 2) "pins should be (2, $N)"
        width, height = N * 0.5sx, 1.5sy
        top = partition(p + (-0.5width, -0.5height), p + (0.5width,- 0.5height), N)
        bot = partition(p + (-0.5width, 0.5height), p + (0.5width, 0.5height), N)
        connections = Dict{Symbol, XYTuple}(
            map(Pair, permutedims(pins, (2, 1))[:], [top..., bot...][:])...
        )

        text_comp = compose(
            context(),
            map(pins[1, :], top) do s, tp
                rot = Rotation(pi/2, tp[1], tp[2] + text_offset*sy)
                text(tp[1], tp[2] + text_offset*sy, string(s), hleft, vcenter, rot)
            end...,
            map(pins[2, :], bot) do s, tp
                rot = Rotation(pi/2, tp[1], tp[2] - text_offset*sy)
                text(tp[1], tp[2] - text_offset*sy, string(s), hright, vcenter, rot)
            end...,
            linewidth(0.0mm), fontsize(0.75_fontsize), font(_font), fill("black")
        )
    end
    composition = compose(
        context(),
        text_comp,
        compose(
            context(),
            text(x, y, name, hcenter, vcenter),
            linewidth(0.0mm), fontsize(_fontsize), font(_font), fill("black")
        ),
        compose(
            context(),
            rectangle(x - 0.5width, y - 0.5height, width, height),
            fill(nothing)
        )
    )

    return Block(composition, connections)
end

# function IC16(
#         x::Real, y::Real, pins::Vector{Symbol}=Symbol[];
#         name::String = "", #rot::Symbol=:h,
#         scale_x::Real=1.0, scale_y::Real=1.0
#     )
#     p = (x, y)
#     sx, sy = scale_x, scale_y
#
#     bottom = partition(p + (-1sx, -0.5sy), p + (1sx, -0.5sy), 8)
#     top = partition(p + (1sx, 0.5sy), p + (-1sx, 0.5sy), 8)
#     connections = if isempty(pins)
#         [bottom..., top...]
#     else
#         @assert (length(pins) == 16) """
#         16 pins should be given, but only $(length(pins)) were.
#         """
#         Dict{Symbol, XYTuple}(
#         map((label, pos) -> Pair(label, pos), pins, [bottom..., top...])...
#         )
#     end
#     composition = compose(
#         context(
#             x-1.05sx, y-0.55sy, 2.1sx, 1.1sy,
#             units = UnitBox(-1.05, -0.55, 2.1, 1.1)
#         ),
#         rectangle(-1.0, -0.5, 2.0, 1.0),
#         render(0.0, 0.0, name)
#     )
#     return Block(composition, connections)
# end

function getindex(B::Block, tag::Symbol)
    if typeof(B.connections) == Dict{Symbol, XYTuple}
        return B.connections[tag]
    else
        throw(ErrorException("Block uses integer indices!"))
    end
end

function getindex(B::Block, i::Integer)
    if typeof(B.connections) == Dict{Symbol, XYTuple}
        throw(ErrorException("Block uses symbol indices!"))
    else
        return B.connections[i]
    end
end


text:
using Flight
using OrdinaryDiffEq
using SciMLBase
using LinearAlgebra
using BenchmarkTools

using StructArrays

using Flight.Kinematics
using Flight.Plotting

atm_sys = System(AtmosphereDescriptor())
p = Geographic(alt = AltOrth(13000))
x_kin = init_x0(KinInit(v_eOb_b = [300, 10, 30], Ob = p, q_nb = Ry(0.01)))
dx_pos = copy(x_kin.pos)
kin_data = f_kin!(dx_pos, x_kin)
air_data = AirData(kin_data, atm_sys)

thr = EThruster();
thr_sys = System(thr);
f_cont!(thr_sys, kin_data, air_data);
y_thr = thr_sys.y
thr_mdl = Model(thr_sys, (kin_data, air_data,))
thr_mdl.u.throttle = 1
step!(thr_mdl, 10, true)

g = AirframeGroup(left = EThruster(), right = EThruster());
g_sys = System(g);
f_cont!(g_sys, kin_data, air_data);
y_g = g_sys.y
g_mdl = Model(g_sys, (kin_data, air,))
g_mdl.u.left.throttle = 1
g_mdl.u.right.throttle = 1
step!(g_mdl, 10, true)


text:
#!/usr/bin/env julia

# Create structured grid VTK file.

using WriteVTK
import Compat: UTF8String, view
typealias FloatType Float32
const vtk_filename_noext = "structured"

function main()
    outfiles = UTF8String[]
    for dim in 2:3
        # Define grid.
        if dim == 2
            Ni, Nj, Nk = 20, 30, 1
        elseif dim == 3
            Ni, Nj, Nk = 20, 30, 40
        end

        local xyz
        if dim == 2
            xyz = zeros(FloatType, 2, Ni, Nj)
            for j = 1:Nj, i = 1:Ni
                xyz[1, i, j] = i/Ni * cos(3*pi/2 * (j-1) / (Nj-1))
                xyz[2, i, j] = i/Ni * sin(3*pi/2 * (j-1) / (Nj-1))
            end
        elseif dim == 3
            xyz = zeros(FloatType, 3, Ni, Nj, Nk)
            for k = 1:Nk, j = 1:Nj, i = 1:Ni
                xyz[1, i, j, k] = i/Ni * cos(3*pi/2 * (j-1) / (Nj-1))
                xyz[2, i, j, k] = i/Ni * sin(3*pi/2 * (j-1) / (Nj-1))
                xyz[3, i, j, k] = (k-1) / Nk
            end
        end

        # Create some scalar and vectorial data.
        q = zeros(FloatType, Ni, Nj, Nk)
        vec = zeros(FloatType, 3, Ni, Nj, Nk)

        # Just to test support for subarrays:
        p = zeros(FloatType, 2Ni, 2Nj, 2Nk)
        psub = view(p, 1:Ni, 1:Nj, 1:Nk)

        for k = 1:Nk, j = 1:Nj, i = 1:Ni
            p[i, j, k] = i*i + k
            q[i, j, k] = k*sqrt(j)
            vec[1, i, j, k] = i
            vec[2, i, j, k] = j
            vec[3, i, j, k] = k
        end

        # Create some scalar data at grid cells.
        # Note that in structured grids, the cells are the hexahedra (3D) or quads (2D)
        # formed between grid points.
        local cdata
        if dim == 2
            cdata = zeros(FloatType, Ni-1, Nj-1)
            for j = 1:Nj-1, i = 1:Ni-1
                cdata[i, j] = 2i + 20 * sin(3*pi * (j-1) / (Nj-2))
            end
        elseif dim == 3
            cdata = zeros(FloatType, Ni-1, Nj-1, Nk-1)
            for k = 1:Nk-1, j = 1:Nj-1, i = 1:Ni-1
                cdata[i, j, k] = 2i + 3k * sin(3*pi * (j-1) / (Nj-2))
            end
        end

        # Test extents (this is optional!!)
        ext = [0, Ni-1, 0, Nj-1, 0, Nk-1]

        # Initialise new vts file (structured grid).
        vtk = vtk_grid(vtk_filename_noext*"_$(dim)D", xyz; extent=ext)

        # This is also accepted:
        # vtk = vtk_grid(vtk_filename_noext, xyz[1, :], xyz[2, :]) # For 2D
        # vtk = vtk_grid(vtk_filename_noext, xyz[1, :], xyz[2, :], xyz[3, :]) # For 3D

        # Add data.
        vtk_point_data(vtk, psub, "p_values")
        vtk_point_data(vtk, q, "q_values")
        vtk_point_data(vtk, vec, "myVector")
        vtk_cell_data(vtk, cdata, "myCellData")

        # Save and close vtk file.
        append!(outfiles, vtk_save(vtk))
    end # dim loop

    println("Saved:", [" "^3 * s for s in outfiles]...)
    return outfiles::Vector{UTF8String}
end

main()


text:
module UtilWrapper

using JavaCall

const ArithmeticUtils = @jimport org.hipparchus.util.ArithmeticUtils
const BigDecimal = @jimport java.math.BigDecimal
const BigInteger = @jimport java.math.BigInteger
const BigReal = @jimport org.hipparchus.util.BigReal
const BigRealField = @jimport org.hipparchus.util.BigRealField
const CalculusFieldElement = @jimport org.hipparchus.CalculusFieldElement
const Class = @jimport java.lang.Class
const Combinations = @jimport org.hipparchus.util.Combinations
const CombinatoricsUtils = @jimport org.hipparchus.util.CombinatoricsUtils
const Comparable = @jimport java.lang.Comparable
const Comparator = @jimport java.util.Comparator
const CompositeFormat = @jimport org.hipparchus.util.CompositeFormat
const Consumer = @jimport java.util.function.Consumer
const ContinuedFraction = @jimport org.hipparchus.util.ContinuedFraction
const Decimal64 = @jimport org.hipparchus.util.Decimal64
const Decimal64Field = @jimport org.hipparchus.util.Decimal64Field
const Enum = @jimport java.lang.Enum
const EventObject = @jimport java.util.EventObject
const FastMath = @jimport org.hipparchus.util.FastMath
const Field = @jimport org.hipparchus.Field
const FieldElement = @jimport org.hipparchus.FieldElement
const FieldPosition = @jimport java.text.FieldPosition
const FieldSinCos = @jimport org.hipparchus.util.FieldSinCos
const FieldTuple = @jimport org.hipparchus.util.FieldTuple
const Incrementor = @jimport org.hipparchus.util.Incrementor
const Incrementor_MaxCountExceededCallback = @jimport org.hipparchus.util.Incrementor$MaxCountExceededCallback
const Iterable = @jimport java.lang.Iterable
const IterationEvent = @jimport org.hipparchus.util.IterationEvent
const IterationListener = @jimport org.hipparchus.util.IterationListener
const IterationManager = @jimport org.hipparchus.util.IterationManager
const JIterator = @jimport java.util.Iterator
const JString = @jimport java.lang.String
const KthSelector = @jimport org.hipparchus.util.KthSelector
const Locale = @jimport java.util.Locale
const Localizable = @jimport org.hipparchus.exception.Localizable
const MathArrays = @jimport org.hipparchus.util.MathArrays
const MathArrays_Function = @jimport org.hipparchus.util.MathArrays$Function
const MathArrays_OrderDirection = @jimport org.hipparchus.util.MathArrays$OrderDirection
const MathArrays_Position = @jimport org.hipparchus.util.MathArrays$Position
const MathContext = @jimport java.math.MathContext
const MathUtils = @jimport org.hipparchus.util.MathUtils
const MultidimensionalCounter = @jimport org.hipparchus.util.MultidimensionalCounter
const MultidimensionalCounter_JIterator = @jimport org.hipparchus.util.MultidimensionalCounter$Iterator
const Number = @jimport java.lang.Number
const NumberFormat = @jimport java.text.NumberFormat
const Object = @jimport java.lang.Object
const OpenIntToDoubleHashMap = @jimport org.hipparchus.util.OpenIntToDoubleHashMap
const OpenIntToDoubleHashMap_JIterator = @jimport org.hipparchus.util.OpenIntToDoubleHashMap$Iterator
const OpenIntToFieldHashMap = @jimport org.hipparchus.util.OpenIntToFieldHashMap
const OpenIntToFieldHashMap_JIterator = @jimport org.hipparchus.util.OpenIntToFieldHashMap$Iterator
const Optional = @jimport java.util.Optional
const Pair = @jimport org.hipparchus.util.Pair
const ParsePosition = @jimport java.text.ParsePosition
const PivotingStrategy = @jimport org.hipparchus.util.PivotingStrategy
const Precision = @jimport org.hipparchus.util.Precision
const RandomGenerator = @jimport org.hipparchus.random.RandomGenerator
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const ResizableDoubleArray = @jimport org.hipparchus.util.ResizableDoubleArray
const ResizableDoubleArray_ExpansionMode = @jimport org.hipparchus.util.ResizableDoubleArray$ExpansionMode
const RoundingMode = @jimport java.math.RoundingMode
const SinCos = @jimport org.hipparchus.util.SinCos
const Spliterator = @jimport java.util.Spliterator
const StringBuffer = @jimport java.lang.StringBuffer
const Tuple = @jimport org.hipparchus.util.Tuple

export ArithmeticUtils
export BigReal
export BigRealField
export Combinations
export CombinatoricsUtils
export CompositeFormat
export ContinuedFraction
export Decimal64
export Decimal64Field
export FastMath
export FastMathCalc
export FastMathLiteralArrays
export FieldSinCos
export FieldTuple
export Incrementor
export IterationEvent
export IterationListener
export IterationManager
export KthSelector
export MathArrays
export MathUtils
export MultidimensionalCounter
export OpenIntToDoubleHashMap
export OpenIntToFieldHashMap
export Pair
export PivotingStrategy
export Precision
export ResizableDoubleArray
export SinCos
export Tuple
export abs
export acos
export acosh
export add
export add_and_check
export add_element
export add_element_rolling
export add_elements
export add_exact
export add_iteration_listener
export asin
export asinh
export atan
export atan2
export atanh
export big_decimal_value
export binomial_coefficient
export binomial_coefficient_double
export binomial_coefficient_log
export build_array
export byte_value
export can_increment
export cbrt
export ceil
export check_binomial
export check_dimension
export check_equal_length
export check_finite
export check_non_negative
export check_not_na_n
export check_not_null
export check_order
export check_positive
export check_range_inclusive
export check_rectangular
export clear
export combinations_iterator
export comparator
export compare_to
export compute
export concatenate
export contains_key
export contract
export convolve
export copy
export copy_sign
export cos
export cos_angle
export cosh
export create
export decrement_exact
export describe_constable
export discard_front_elements
export discard_most_recent_elements
export distance
export distance1
export distance_inf
export divide
export divide_unsigned
export double_value
export ebe_add
export ebe_divide
export ebe_multiply
export ebe_subtract
export equals
export equals_including_na_n
export equals_with_relative_tolerance
export evaluate
export exp
export expm1
export factorial
export factorial_double
export factorial_log
export fire_initialization_event
export fire_iteration_performed_event
export fire_iteration_started_event
export fire_termination_event
export float_value
export floor
export floor_div
export floor_mod
export fma
export for_each
export format_double
export gcd
export get
export get_capacity
export get_class
export get_component
export get_components
export get_contraction_criterion
export get_count
export get_counts
export get_declaring_class
export get_default_number_format
export get_dimension
export get_element
export get_elements
export get_expansion_factor
export get_expansion_mode
export get_exponent
export get_field
export get_first
export get_instance
export get_iterations
export get_k
export get_key
export get_max_iterations
export get_maximal_count
export get_n
export get_num_elements
export get_one
export get_pivoting_strategy
export get_real
export get_rounding_mode
export get_runtime_class
export get_scale
export get_second
export get_size
export get_sizes
export get_source
export get_value
export get_zero
export hash
export hash_code
export hypot
export i_e_e_eremainder
export increment
export increment_exact
export increment_iteration_count
export initialization_performed
export int_value
export is_infinite
export is_mathematical_integer
export is_monotonic
export is_na_n
export is_power_of_two
export iteration_performed
export iteration_started
export iterator
export lcm
export linear_combination
export log
export log10
export log1p
export long_value
export main
export max
export min
export mul_and_check
export multiply
export multiply_exact
export multiply_full
export multiply_high
export name
export natural
export negate
export new_instance
export next_after
export next_down
export next_up
export normalize_angle
export normalize_array
export notify
export notify_all
export ordinal
export parse_and_ignore_whitespace
export parse_fixedstring
export parse_next_character
export parse_number
export pivot_index
export pow
export put
export random
export reciprocal
export reduce
export remainder
export remainder_unsigned
export remove
export remove_iteration_listener
export representable_delta
export reset
export reset_iteration_count
export rint
export root_n
export round
export safe_norm
export scalb
export scale
export scale_in_place
export select
export sequence
export set_element
export set_num_elements
export set_rounding_mode
export set_scale
export short_value
export shuffle
export signum
export sin
export sin_cos
export sinh
export size
export sort_in_place
export spliterator
export sqrt
export stirling_s2
export sub_and_check
export substitute_most_recent_element
export subtract
export subtract_exact
export tan
export tanh
export termination_performed
export to_degrees
export to_int_exact
export to_radians
export to_string
export ulp
export unique
export value_of
export values
export verify_values
export wait
export with_callback
export with_count
export with_maximal_count

include("arithmetic_utils.jl")
include("big_real.jl")
include("big_real_field.jl")
include("combinations.jl")
include("combinatorics_utils.jl")
include("composite_format.jl")
include("continued_fraction.jl")
include("decimal64.jl")
include("decimal64_field.jl")
include("fast_math.jl")
include("fast_math_calc.jl")
include("fast_math_literal_arrays.jl")
include("field_sin_cos.jl")
include("field_tuple.jl")
include("incrementor.jl")
include("iteration_event.jl")
include("iteration_listener.jl")
include("iteration_manager.jl")
include("kth_selector.jl")
include("math_arrays.jl")
include("math_utils.jl")
include("multidimensional_counter.jl")
include("open_int_to_double_hash_map.jl")
include("open_int_to_field_hash_map.jl")
include("pair.jl")
include("pivoting_strategy.jl")
include("precision.jl")
include("resizable_double_array.jl")
include("sin_cos.jl")
include("tuple.jl")

end


text:
function embiggen(arr, increase)
  result = falses(increase*2 .+ size(arr))
  result[[1:n for n=size(arr)]...] = arr
  return result
end

neighbours = [(x, y, z) for x=-1:1, y=-1:1, z=-1:1 if (x,y,z) != (0,0,0)]

function neighbour_count(arr)
  return sum(map(shift -> circshift(arr, shift), neighbours))
end

function evolve(x, n)
  x = embiggen(x, n)
  for it=1:n
    cs = neighbour_count(x)
    x = (cs .== 3) .| (x .& (cs .== 2))
  end
  return x
end

input="##..#.#.\n#####.##\n#######.\n#..#..#.\n#.#...##\n..#....#\n....#..#\n..##.#.."

init = hcat([[x == '#' for x=line] for line=split(input)]...)
init = reshape(BitArray(init), size(init)..., 1)
println(sum(evolve(init, 6)))

max_stars_repo_name: Porges/adventofcode

max_issues_repo_name: Porges/adventofcode

max_forks_repo_name: Porges/adventofcode


text:
# functions for calculating the flux for interior face integrals

@doc """
### EulerEquationMod.calcFaceFlux

  This function calculates the DG flux between a specified set of faces,
  using the solution data at the faces stored in eqn.q_face.
  Note that the flux is negated because the face integrals have a
  negative sign in the weak form.

  Conservative variables only!

  **Inputs**:

   * mesh
   * sbp
   * eqn
   * functor: the functor that calculates the flux at a node
   * interfaces: an array of type Interface that specifies which interfaces
                to calculate the flux for

  **Inputs/Outputs**:
   * face_flux: array to store the flux in, numDofPerNode x nnodesPerFace
               x length(interfaces)

  The functor must have the signature:

  `func( uL, qR, aux_vars, dxidx, nrm, flux_j, eqn.params)`

  where uL and uR are the solution values for a node on the left and right
  elements, aux_vars are the auxiliary variables for the node,
  dxidx is the scaled mapping jacobian for elementL, and nrm is the face
  normal in reference space. flux_j is the array of length numDofPerNode to be
  populated with the flux. params is eqn.params.

"""->
function calcFaceFlux( mesh::AbstractDGMesh{Tmsh},
  sbp::AbstractOperator,
  eqn::EulerData{Tsol, Tres, Tdim, :conservative},
  functor::FluxType,
  interfaces::AbstractArray{Interface,1},
  face_flux::AbstractArray{Tres, 3}) where {Tmsh, Tsol, Tres, Tdim}

  nfaces = length(interfaces)
  params = eqn.params
  for i=1:nfaces  # loop over faces
    interface_i = interfaces[i]
    for j = 1:mesh.numNodesPerFace

      eL = interface_i.elementL
      fL = interface_i.faceL
      # get components
      qL = ro_sview(eqn.q_face, :, 1, j, i)
      qR = ro_sview(eqn.q_face, :, 2, j, i)
      aux_vars = ro_sview(eqn.aux_vars_face, :, j, i)
      nrm_xy = ro_sview(mesh.nrm_face, :, j, i)
      flux_j = sview(face_flux, :, j, i)

      functor(params, qL, qR, aux_vars, nrm_xy, flux_j)
    end
  end

  return nothing
end

"""
  Like [`calcFaceFlux`](@ref), but computes the flux for a single element and
  then integrates it immediately, updating eqn.res

  **Inputs**:
   * mesh
   * sbp
   * eqn
   * opts
   * functor: a FluxType that evalutes the flux
   * interfaces: the vector of [`Interface`](@ref)s to compute the integrals for

"""
function calcFaceIntegral_nopre(
        mesh::AbstractDGMesh{Tmsh},
        sbp::AbstractOperator,
        eqn::EulerData{Tsol, Tres, Tdim, :conservative},
        opts,
        functor::FluxType,
        interfaces::AbstractArray{Interface, 1}) where {Tmsh, Tsol, Tres, Tdim}

  nfaces = length(interfaces)
  params = eqn.params

  q_faceL = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceR = zeros(q_faceL)

  flux_face = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)

  for i=1:nfaces
    iface_i = interfaces[i]
#    fill!(flux_face, 0)

    qL = ro_sview(eqn.q, :, :, iface_i.elementL)
    qR = ro_sview(eqn.q, :, :, iface_i.elementR)
    interiorFaceInterpolate!(mesh.sbpface, iface_i, qL, qR, q_faceL,
                             q_faceR)

    for j=1:mesh.numNodesPerFace
      qL_j = ro_sview(q_faceL, :, j)
      qR_j = ro_sview(q_faceR, :, j)

      eqn.aux_vars_face[1, j, i] = calcPressure(params, qL_j)
      aux_vars = ro_sview(eqn.aux_vars_face, :, j, i)

      nrm_xy = ro_sview(mesh.nrm_face, :, j, i)
      flux_j = sview(flux_face, :, j)

      functor(params, qL_j, qR_j, aux_vars, nrm_xy, flux_j)
    end  # end loop j

    resL = sview(eqn.res, :, :, iface_i.elementL)
    resR = sview(eqn.res, :, :, iface_i.elementR)
    interiorFaceIntegrate!(mesh.sbpface, iface_i, flux_face, resL, resR,
                           SummationByParts.Subtract())
  end  # end loop i

  return nothing
end

"""
  This function loops over interfaces and computes a face integral that
  uses data from all volume nodes. See [`FaceElementIntegralType`](@ref)
  for details on the integral performed.
"""
function getFaceElementIntegral(
                           mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
                           face_integral_functor::FaceElementIntegralType,
                           flux_functor::FluxType,
                           sbpface::AbstractFace,
                           interfaces::AbstractArray{Interface, 1}) where {Tmsh, Tsol, Tres, Tdim}

  params = eqn.params
  nfaces = length(interfaces)

  for i=1:nfaces
    iface = interfaces[i]
    elL = iface.elementL
    elR = iface.elementR
    qL = ro_sview(eqn.q, :, :, elL)
    qR = ro_sview(eqn.q, :, :, elR)
    aux_vars = ro_sview(eqn.aux_vars, :, :, elL)
    nrm_face = ro_sview(mesh.nrm_face, :, :, i)
    resL = sview(eqn.res, :, :, elL)
    resR = sview(eqn.res, :, :, elR)

    calcFaceElementIntegral(face_integral_functor, params, sbpface, iface, qL,
                            qR, aux_vars, nrm_face, flux_functor, resL, resR)

  end

  return nothing
end

"""
  This function loops over interfaces and computes a face integral that
  uses data from all volume nodes. See FaceElementIntegralType for details on
  the integral performed. This function works on the staggered grid.

  Inputs:
    mesh_s: mesh object on solution mesh
    mesh_f: mesh object on flux mesh
    sbp_s: SBP operator on solution mesh
    sbp_f: SBP operator on flux mesh
    face_integral_functor: a [`FaceElementIntegralType`](@doc) functor
    flux_functor: a [`FluxType`] functor that is passed to face_integral_functor
                  to use as the numerical flux function
    sbpface: an AbstractFace for the flux grid
    interfaces: vector of [`Interfaces`](@doc) that the integral will be
                computed for.

  Inputs/Outputs:
    eqn: equation object (implicitly lives on solution grid).  eqn.res is
          updated with the results.
"""
function getFaceElementIntegral(
                           mesh_s::AbstractDGMesh{Tmsh},
                           mesh_f::AbstractDGMesh{Tmsh},
                           sbp_s::AbstractOperator, sbp_f::AbstractOperator,
                           eqn::EulerData{Tsol, Tres, Tdim},
                           face_integral_functor::FaceElementIntegralType,
                           flux_functor::FluxType,
                           sbpface::AbstractFace,
                           interfaces::AbstractArray{Interface, 1}) where {Tmsh, Tsol, Tres, Tdim}
  params = eqn.params
  qf = eqn.q_flux
  res = eqn.res
  nfaces = length(interfaces)
  data = params.face_element_integral_data
  @unpack data resL_s resR_s resL_f resR_f

  aux_vars = zeros(Tres, 1, mesh_f.numNodesPerElement)

  for i=1:nfaces
    iface = interfaces[i]
    elL = iface.elementL
    elR = iface.elementR
    nrm_face = ro_sview(mesh_f.nrm_face, :, :, i)

    qf_L = ro_sview(qf, :, :, elL)
    qf_R = ro_sview(qf, :, :, elR)

    for j=1:mesh_f.numNodesPerElement
      aux_vars[1, j] = calcPressure(params, ro_sview(qf_L, :, j))
    end

    fill!(resL_f, 0.0)
    fill!(resR_f, 0.0)

    calcFaceElementIntegral(face_integral_functor, params, sbpface, iface,
                            qf_L, qf_R, aux_vars, nrm_face, flux_functor,
                            resL_f, resR_f)

    # interpolate back
    smallmatmat!(resL_f, mesh_s.I_S2F, resL_s)
    smallmatmat!(resR_f, mesh_s.I_S2F, resR_s)

    # accumulate into res
    @simd for j=1:mesh_s.numNodesPerElement
      @simd for k=1:mesh_s.numDofPerNode
        res[k, j, elL] += resL_s[k, j]
        res[k, j, elR] += resR_s[k, j]
      end
    end

  end  # end loop i

  return nothing
end

#------------------------------------------------------------------------------
# parallel face integrals (face parallel)

"""
  This function is a thin wrapper around calcSharedFaceIntegrals_inner.
  It present the interface needed by [`finishExchangeData`](@ref).
"""
function calcSharedFaceIntegrals( mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol},
                            opts, data::SharedFaceData) where {Tmsh, Tsol}

  # if there were other ways to doing the shared face integrals, there would
  # be a conditional here
  calcSharedFaceIntegrals_nopre_inner(mesh, sbp, eqn, opts, data, eqn.flux_func)
  return nothing
end

"""
  Like [`calcSharedFaceIntegrals_inner`](@ref), but performs the integration and
  updates eqn.res rather than computing the flux and storing it in
  eqn.flux_sharedface
"""
function calcSharedFaceIntegrals_nopre_inner(
                            mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres},
                            opts, data::SharedFaceData, functor::FluxType) where {Tmsh, Tsol, Tres}

  if getParallelData(data) != PARALLEL_DATA_FACE
    throw(ErrorException("cannot use calcSharedFaceIntegrals without parallel face data"))
  end

  params = eqn.params

  # calculate the flux
  idx = data.peeridx
  interfaces = data.interfaces
  qL_arr = data.q_send
  qR_arr = data.q_recv
  aux_vars_arr = eqn.aux_vars_sharedface[idx]
  nrm_arr = mesh.nrm_sharedface[idx]
  flux_j = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)

  for j=1:length(interfaces)
    interface_i = interfaces[j]
    eL = interface_i.elementL
    fL = interface_i.faceL

    # compute the flux
    for k=1:mesh.numNodesPerFace
      qL = ro_sview(qL_arr, :, k, j)
      qR = ro_sview(qR_arr, :, k, j)
      aux_vars = ro_sview(aux_vars_arr, :, k, j)
      parent(aux_vars)[1] = calcPressure(params, qL)
      nrm_xy = ro_sview(nrm_arr, :, k, j)
      flux_k = sview(flux_j, :, k)
      functor(params, qL, qR, aux_vars, nrm_xy, flux_k)
    end

    # do the integration
    res_j = sview(eqn.res, :, :, eL)
    boundaryFaceIntegrate!(mesh.sbpface, fL, flux_j, res_j, SummationByParts.Subtract())

  end

  return nothing
end

#------------------------------------------------------------------------------
# parallel face integrals (element parallel)

"""
  This function is a thin wrapper around
  [`calcSharedFaceIntegrals_element_inner`](@ref).
  It presents the interface required by [`finishExchangeData`](@ref)
"""
function calcSharedFaceIntegrals_element(mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol},
                            opts, data::SharedFaceData) where {Tmsh, Tsol}

  calcSharedFaceIntegrals_nopre_element_inner(mesh, sbp, eqn, opts, data, eqn.flux_func)

  return nothing
end

"""
  Like [`calcSharedFaceIntegrals_element_inner`](@ref), but performs the integration and
  updates eqn.res rather than computing the flux only and storing it in
  eqn.flux_sharedface
"""
function calcSharedFaceIntegrals_nopre_element_inner(
                            mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres},
                            opts, data::SharedFaceData, functor::FluxType) where {Tmsh, Tsol, Tres}

  q = eqn.q
  params = eqn.params

  # TODO: make these fields of params
  q_faceL = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceR = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)

  idx = data.peeridx
  interfaces = data.interfaces
  bndries_local = data.bndries_local
  bndries_remote = data.bndries_remote

  # eqn.q_face_send should actually contain el
  # data, but we want to use eqn.q because that's the canonical source
  # qL_arr = eqn.q_face_send[i]
  qR_arr = data.q_recv
  nrm_arr = mesh.nrm_sharedface[idx]
  aux_vars_arr = eqn.aux_vars_sharedface[idx]
  flux_face = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)

  start_elnum = mesh.shared_element_offsets[idx]

  for j=1:length(interfaces)
    iface_j = interfaces[j]
    bndryL_j = bndries_local[j]
    bndryR_j = bndries_remote[j]
    fL = bndryL_j.face

    # interpolate to face
    qL = ro_sview(q, :, :, iface_j.elementL)
    el_r = iface_j.elementR - start_elnum + 1
    qR = ro_sview(qR_arr, :, :, el_r)
    interiorFaceInterpolate!(mesh.sbpface, iface_j, qL, qR, q_faceL, q_faceR)

    # calculate flux
    for k=1:mesh.numNodesPerFace
      qL_k = ro_sview(q_faceL, :, k)
      qR_k = ro_sview(q_faceR, :, k)
      aux_vars = ro_sview(aux_vars_arr, :, k, j)
      nrm_xy = ro_sview(nrm_arr, :, k, j)
      flux_k = sview(flux_face, :, k)

      parent(aux_vars)[1] = calcPressure(params, qL_k)

      functor(params, qL_k, qR_k, aux_vars, nrm_xy, flux_k)
     end

     # do the integration
     res_j = sview(eqn.res, :, :, bndryL_j.element)
     boundaryFaceIntegrate!(mesh.sbpface, fL, flux_face, res_j, SummationByParts.Subtract())

   end  # end loop over interfaces

  return nothing
end

#------------------------------------------------------------------------------
# parallel face element integrals

"""
  This function is a thin wrapper around
  [`calcSharedFaceElementIntegrals_element_inner`](@ref),
  presenting the interface needed by [`finishExchangeData`](@ref).
  See that function for the interface details.
"""
function calcSharedFaceElementIntegrals_element(
                            mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres},
                            opts, data::SharedFaceData) where {Tmsh, Tsol, Tres}

  if opts["use_staggered_grid"]
    calcSharedFaceElementIntegralsStaggered_element_inner(mesh, mesh.mesh2,
                            sbp, mesh.sbp2, eqn, opts, data,
                            eqn.face_element_integral_func, eqn.flux_func)

  else
    calcSharedFaceElementIntegrals_element_inner(mesh, sbp, eqn, opts, data,
                            eqn.face_element_integral_func,  eqn.flux_func)
  end

  return nothing
end

"""
  This function loops over given set of shared faces and computes a face
  integral that
  uses data from all volume nodes.  See [`FaceElementIntegralType`](@ref)
  for details on the integral performed.

  **Inputs**:

   * mesh
   * sbp
   * eqn
   * opts
   * data: a SharedFaceData specifying which shared faces to compute
   * face_integral_functor
   * flux_functor

"""
function calcSharedFaceElementIntegrals_element_inner(
                            mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres},
                            opts, data::SharedFaceData,
                            face_integral_functor::FaceElementIntegralType,
                            flux_functor::FluxType) where {Tmsh, Tsol, Tres}

  if getParallelData(data) != PARALLEL_DATA_ELEMENT
    throw(ErrorException("cannot use calcSharedFaceIntegrals_element without parallel element data"))
  end

  q = eqn.q
  params = eqn.params
  # we don't care about elementR here, so use this throwaway array
  resR = Array{Tres}(mesh.numDofPerNode, mesh.numNodesPerElement)

  # get the data for the parallel interface
  idx = data.peeridx
  interfaces = data.interfaces
  bndries_local = data.bndries_local
  bndries_remote = data.bndries_remote
#    qL_arr = eqn.q_face_send[i]
  qR_arr = data.q_recv
  nrm_face_arr = mesh.nrm_sharedface[idx]

  writedlm("qR$(data.peernum)_$(mesh.myrank)_primal.dat", real(qR_arr))

  start_elnum = mesh.shared_element_offsets[idx]
  #TODO debugging
  # compute the integrals
  for j=1:length(interfaces)
    iface_j = interfaces[j]
    elL = iface_j.elementL
    elR = iface_j.elementR - start_elnum + 1  # is this always equal to j?
    qL = ro_sview(q, :, :, elL)
    qR = ro_sview(qR_arr, :, :, elR)
    aux_vars = ro_sview(eqn.aux_vars, :, :, elL)
    nrm_face = ro_sview(nrm_face_arr, :, :, j)
    resL = sview(eqn.res, :, :, elL)

    #=
    for k=1:mesh.numNodesPerElement
      for p=1:mesh.numDofPerNode
        resL[p, k] += qL[p, k] + qR[p, k]
      end
    end
    =#

    calcFaceElementIntegral(face_integral_functor, eqn.params, mesh.sbpface,
                            iface_j, qL, qR, aux_vars,
                            nrm_face, flux_functor, resL, resR)
  end  # end loop j

  return nothing
end

"""
  Like [`calcSharedFaceElementIntegrals_element_inner`](@ref), but for staggered grid.

  data.q_recv is the solution grid data.  This function interpolates it to the
  flux grid on the fly.

  **Inputs**:

   * mesh_s: solution grid mesh
   * mesh_f: flux grid mesh
   * sbp_s: SBP operator for solution grid
   * sbp_f: SBP operator for the flux grid
   * opts: options dictionary
   * data: [`SharedFaceData`](@ref)
   * face_integral_functor: [`FaceElementIntegralType`](@ref)
   * flux_functor: [`FluxType`](@ref) passed to face_integral functor

  **Inputs/Outputs**:

    * eqn: equation object (lives on solution grid)

  Aliasing restrictions: none
"""
function calcSharedFaceElementIntegralsStaggered_element_inner(
                            mesh_s::AbstractDGMesh{Tmsh},
                            mesh_f::AbstractDGMesh{Tmsh},
                            sbp_s::AbstractOperator, sbp_f::AbstractOperator,
                            eqn::EulerData{Tsol, Tres},
                            opts, data::SharedFaceData,
                            face_integral_functor::FaceElementIntegralType,
                            flux_functor::FluxType) where {Tmsh, Tsol, Tres}

  if getParallelData(data) != PARALLEL_DATA_ELEMENT
    throw(ErrorException("cannot use calcSharedFaceIntegrals_element without parallel element data"))
  end

  q = eqn.q_flux
  res = eqn.res
  params = eqn.params

  # temporary arrays needed for interpolating to the flux grid
  fdata = params.face_element_integral_data
  @unpack fdata qvars_f resL_f resR_f resL_s

  # we don't care about resR here
  aux_vars = Array{Tres}(1, mesh_f.numNodesPerElement)

  # get the data for the parallel interface
  idx = data.peeridx
  interfaces = data.interfaces
  qR_arr = data.q_recv
  nrm_face_arr = mesh_f.nrm_sharedface[idx]
  start_elnum = mesh_f.shared_element_offsets[idx]

  # compute the integrals
  for j=1:length(interfaces)
    iface_j = interfaces[j]
    elL = iface_j.elementL
    elR = iface_j.elementR - start_elnum + 1
    qL = ro_sview(q, :, :, elL)
    qR = ro_sview(qR_arr, :, :, elR)
#    aux_vars = ro_sview(eqn.aux_vars, :, :, elL)
    nrm_face = ro_sview(nrm_face_arr, :, :, j)

    # interpolate to flux grid
    interpolateElementStaggered(params, mesh_s, qR, aux_vars, qvars_f)

    fill!(resL_f, 0.0)
    # we dont carea bout resR, so dont waste time zeroing it out
    calcFaceElementIntegral(face_integral_functor, eqn.params, mesh_f.sbpface,
                            iface_j, qL, qvars_f, aux_vars, nrm_face,
                            flux_functor, resL_f, resR_f)

    # interpolate residual back to solution grid
    smallmatmat!(resL_f, mesh_s.I_S2F, resL_s)

    @simd for k=1:mesh_s.numNodesPerElement
      @simd for p=1:mesh_s.numDofPerNode
        res[p, k, elL] += resL_s[p, k]
      end
    end

  end  # end loop j

  return nothing
end

@doc """
### EulerEquationMod.interpolateFace

  This function interpolates the solution values from the internal nodes
  to the face flux points of the elements

  **Inputs**:

   * mesh: an AbstractDGMesh
   * sbp
   * eqn
   * opts
   * q: a 3D array of solution values at the nodes, numDofPerNode x
       numNodesPerElement x numEl

  **Inputs/Outputs**:

   * q_face: a 4D array of solution values at each interface,
            numDofPerNode x 2 x numfacenodes x numInterface
            q_face[:, 1, j, i] stores the q values for elementL of interface
            i node j and q_face[:, 2, j, i] stores the values for elementR

  eqn.aux_vars_face is also populated
"""->
function interpolateFace(mesh::AbstractDGMesh, sbp, eqn, opts,
                   q::Abstract3DArray, q_face::AbstractArray{Tsol, 4}) where Tsol

  # interpolate solution
  interiorfaceinterpolate!(mesh.sbpface, mesh.interfaces, q, q_face)

  # recalculte aux_vars
  for i=1:mesh.numInterfaces
    for j=1:mesh.numNodesPerFace
      q_vals = ro_sview(q_face, :, 1, j, i) # always use elementL
      eqn.aux_vars_face[1, j, i] = calcPressure(eqn.params, q_vals)
    end
  end

  if opts["parallel_data"] == PARALLEL_DATA_FACE
    for peer=1:mesh.npeers
      q_vals_p = eqn.shared_data[peer].q_send
      aux_vars = eqn.aux_vars_sharedface[peer]
      for i=1:mesh.peer_face_counts[peer]
        for j=1:mesh.numNodesPerFace
          q_vals = ro_sview(q_vals_p, :, j, i)
          aux_vars[1, j, i] = calcPressure(eqn.params, q_vals)
        end
      end
    end
  end  # else populate aux_vars_sharedface in evalSharedFaceIntegrals

  return nothing
end

#------------------------------------------------------------------------------
# Different fluxes
#------------------------------------------------------------------------------

"""
  This flux function throws an error. Useful for defaults.
"""
mutable struct ErrorFlux <: FluxType
end

function (obj::ErrorFlux)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractArray,
              F::AbstractArray{Tres}) where {Tsol, Tres}

  error("ErrorFlux called.")
  return nothing
end

mutable struct ErrorFlux_revm <: FluxType_revm
end

function (obj::ErrorFlux_revm)(params::ParamType,
                  qL::AbstractArray{Tsol,1}, qR::AbstractArray{Tsol, 1},
                  aux_vars::AbstractArray{Tres},
                  nrm::AbstractArray{Tmsh}, nrm_bar::AbstractArray{Tmsh},
                  F_bar::AbstractArray{Tres}) where {Tmsh, Tsol, Tres}

  error("ErrorFlux_revm called")

  return nothing
end

mutable struct ErrorFlux_revq <: FluxType_revq
end

function (obj::ErrorFlux_revq)(params::ParamType,
                      qL::AbstractArray{Tsol,1}, qL_bar::AbstractArray{Tsol, 1},
                      qR::AbstractArray{Tsol, 1}, qR_bar::AbstractArray{Tsol, 1},
                      aux_vars::AbstractArray{Tres}, dir::AbstractArray{Tmsh},
                      F_bar::AbstractArray{Tres}) where {Tmsh, Tsol, Tres}

  error("ErrorFlux_revq called")

  return nothing
end

mutable struct ZeroFlux <: FluxType
end

function (obj::ZeroFlux)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractArray,
              F::AbstractArray{Tres}) where {Tsol, Tres}

  fill!(F, 0)

  return nothing
end

mutable struct ZeroFlux_revm <: FluxType_revm
end

function (obj::ZeroFlux_revm)(params::ParamType,
                  qL::AbstractArray{Tsol,1}, qR::AbstractArray{Tsol, 1},
                  aux_vars::AbstractArray{Tres},
                  nrm::AbstractArray{Tmsh}, nrm_bar::AbstractArray{Tmsh},
                  F_bar::AbstractArray{Tres}) where {Tmsh, Tsol, Tres}

  return nothing
end

mutable struct ZeroFlux_revq <: FluxType_revq
end

function (obj::ZeroFlux_revq)(params::ParamType,
                      qL::AbstractArray{Tsol,1}, qL_bar::AbstractArray{Tsol, 1},
                      qR::AbstractArray{Tsol, 1}, qR_bar::AbstractArray{Tsol, 1},
                      aux_vars::AbstractArray{Tres}, dir::AbstractArray{Tmsh},
                      F_bar::AbstractArray{Tres}) where {Tmsh, Tsol, Tres}

  return nothing
end

"""
  This flux function sets F = q.  Useful for testing
"""
mutable struct IdentityFlux <: FluxType
end

function (obj::IdentityFlux)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractVector{Tmsh},
              F::AbstractVector{Tres}) where {Tsol, Tres, Tmsh}

  for i=1:length(F)
    F[i] = q[i]
  end

  return nothing
end

function (obj::IdentityFlux)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractMatrix{Tmsh},
              F::AbstractMatrix{Tres}) where {Tsol, Tres, Tmsh}

  for j=1:size(F, 2)
    for i=1:length(F)
      F[i, j] = q[i]
    end
  end

  return nothing
end

"""
  Calls the [`RoeSolver`](@ref)
"""

mutable struct RoeFlux <: FluxType
end

function (obj::RoeFlux)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractVector{Tmsh},
              F::AbstractVector{Tres}) where {Tsol, Tres, Tmsh}

  RoeSolver(params, uL, uR, aux_vars, nrm, F)
  return nothing
end

mutable struct RoeFlux_revm <: FluxType_revm
end

function (obj::RoeFlux_revm)(params::ParamType,
              uL::AbstractArray{Tsol,1}, uR::AbstractArray{Tsol, 1}, aux_vars,
              nrm::AbstractVector{Tmsh}, nrm_bar::AbstractArray{Tmsh, 1},
              flux_bar::AbstractVector{Tres}) where {Tsol, Tres, Tmsh}

  RoeSolver_revm(params, uL, uR, aux_vars, nrm, nrm_bar, flux_bar)

  return nothing
end

mutable struct RoeFlux_revq <: FluxType_revq
end

function (obj::RoeFlux_revq)(params::ParamType,
                      qL::AbstractArray{Tsol,1}, qL_bar::AbstractArray{Tsol, 1},
                      qR::AbstractArray{Tsol, 1}, qR_bar::AbstractArray{Tsol, 1},
                      aux_vars::AbstractArray{Tres}, dir::AbstractArray{Tmsh},
                      F_bar::AbstractArray{Tres}) where {Tmsh, Tsol, Tres}

  RoeSolver_revq(params, qL, qL_bar, qR, qR_bar, aux_vars, dir, F_bar)

  return nothing
end

"""
  Calls [`calcEulerFlux_standard`](@ref)
"""
mutable struct LFFlux <: FluxType
end

function (obj::LFFlux)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractVector{Tmsh},
              F::AbstractVector{Tres}) where {Tsol, Tres, Tmsh}

  calcLFFlux(params, uL, uR, aux_vars, nrm, F)
  return nothing
end

"""
  calls [`calcEulerFlux_IR_revm`](@ref)
"""
mutable struct LFFlux_revm <: FluxType_revm
end

function (obj::LFFlux_revm)(params::ParamType,
                  qL::AbstractArray{Tsol,1}, qR::AbstractArray{Tsol, 1},
                  aux_vars::AbstractArray{Tres},
                  nrm::AbstractArray{Tmsh}, nrm_bar::AbstractArray{Tmsh},
                  F_bar::AbstractArray{Tres}) where {Tmsh, Tsol, Tres}

  calcLFFlux_revm(params, qL, qR, aux_vars, nrm, nrm_bar, F_bar)

  return nothing
end

mutable struct LFFlux_revq <: FluxType_revq
end

function (obj::LFFlux_revq)(params::ParamType,
                      qL::AbstractArray{Tsol,1}, qL_bar::AbstractArray{Tsol, 1},
                      qR::AbstractArray{Tsol, 1}, qR_bar::AbstractArray{Tsol, 1},
                      aux_vars::AbstractArray{Tres}, dir::AbstractArray{Tmsh},
                      F_bar::AbstractArray{Tres}) where {Tmsh, Tsol, Tres}

  calcLFFlux_revq(params, qL, qL_bar, qR, qR_bar, aux_vars, dir,
                        F_bar)

  return nothing
end

mutable struct StandardFlux <: FluxType
end

function (obj::StandardFlux)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractArray,
              F::AbstractArray{Tres}) where {Tsol, Tres}

  calcEulerFlux_standard(params, uL, uR, aux_vars, nrm, F)
  return nothing
end

"""
  Calls [`calcEulerFlux_Ducros`](@ref)
"""
mutable struct DucrosFlux <: FluxType
end

function (obj::DucrosFlux)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractVector,
              F::AbstractVector{Tres}) where {Tsol, Tres}

  calcEulerFlux_Ducros(params, uL, uR, aux_vars, nrm, F)
  return nothing
end

"""
  Calls [`calcEulerFlux_IR`](@ref)
"""
mutable struct IRFlux <: FluxType
end

function (obj::IRFlux)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractArray,
              F::AbstractArray{Tres}) where {Tsol, Tres}

  # this will dispatch to either the sinlge director or multi-dimension method
  calcEulerFlux_IR(params, uL, uR, aux_vars, nrm, F)
  return nothing
end

"""
  calls [`calcEulerFlux_IR_revm`](@ref)
"""
mutable struct IRFlux_revm <: FluxType_revm
end

function (obj::IRFlux_revm)(params::ParamType,
                  qL::AbstractArray{Tsol,1}, qR::AbstractArray{Tsol, 1},
                  aux_vars::AbstractArray{Tres},
                  nrm::AbstractArray{Tmsh}, nrm_bar::AbstractArray{Tmsh},
                  F_bar::AbstractArray{Tres}) where {Tmsh, Tsol, Tres}

  calcEulerFlux_IR_revm(params, qL, qR, aux_vars, nrm, nrm_bar, F_bar)

  return nothing
end

mutable struct IRFlux_revq <: FluxType_revq
end

function (obj::IRFlux_revq)(params::ParamType,
                      qL::AbstractArray{Tsol,1}, qL_bar::AbstractArray{Tsol, 1},
                      qR::AbstractArray{Tsol, 1}, qR_bar::AbstractArray{Tsol, 1},
                      aux_vars::AbstractArray{Tres}, dir::AbstractArray{Tmsh},
                      F_bar::AbstractArray{Tres}) where {Tmsh, Tsol, Tres}

  calcEulerFlux_IR_revq(params, qL, qL_bar, qR, qR_bar, aux_vars, dir,
                        F_bar)

  return nothing
end

"""
  Calls [`calcEulerFlux_IRSLF`](@ref)
"""
mutable struct IRSLFFlux <: FluxType
end

function (obj::IRSLFFlux)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractVector,
              F::AbstractVector{Tres}) where {Tsol, Tres}

  calcEulerFlux_IRSLF(params, uL, uR, aux_vars, nrm, F)
  return nothing
end

mutable struct IRSLFFlux_revm <: FluxType_revm
end

function (obj::IRSLFFlux_revm)(params::ParamType,
                  qL::AbstractArray{Tsol,1}, qR::AbstractArray{Tsol, 1},
                  aux_vars::AbstractArray{Tres},
                  nrm::AbstractArray{Tmsh}, nrm_bar::AbstractArray{Tmsh},
                  F_bar::AbstractArray{Tres}) where {Tmsh, Tsol, Tres}

  kernel = params.entropy_lf_kernel
  calcEulerFlux_IR_revm(params, qL, qR, aux_vars, nrm, nrm_bar, F_bar)
  applyEntropyKernel_diagE_revm(params, kernel, qL, qR, aux_vars, nrm, nrm_bar,
                                F_bar)

  return nothing
end

mutable struct IRSLFFlux_revq <: FluxType_revq
end

function (obj::IRSLFFlux_revq)(params::ParamType,
                      qL::AbstractArray{Tsol,1}, qL_bar::AbstractArray{Tsol, 1},
                      qR::AbstractArray{Tsol, 1}, qR_bar::AbstractArray{Tsol, 1},
                      aux_vars::AbstractArray{Tres}, nrm::AbstractArray{Tmsh},
                      F_bar::AbstractArray{Tres}) where {Tmsh, Tsol, Tres}

  kernel = params.entropy_lf_kernel
  calcEulerFlux_IR_revq(params, qL, qL_bar, qR, qR_bar, aux_vars, nrm,
                        F_bar)
  applyEntropyKernel_diagE_revq(params, kernel, qL, qL_bar, qR, qR_bar,
                                aux_vars, nrm, F_bar)
  return nothing
end

"""
  Computes only the penalty term from the [`IRSLFFlux`](@ref)
"""
mutable struct LFPenalty <: FluxType
end

function (obj::LFPenalty)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractVector,
              F::AbstractVector{Tres}) where {Tsol, Tres}

  kernel = params.entropy_lf_kernel
  fill!(F, 0)  # fluxes must *overwrite* the F array
  applyEntropyKernel_diagE(params, kernel, uL, uR, aux_vars, nrm, F)

  return nothing
end

mutable struct LFPenalty_revm <: FluxType_revm
end

function (obj::LFPenalty_revm)(params::ParamType,
                  qL::AbstractArray{Tsol,1}, qR::AbstractArray{Tsol, 1},
                  aux_vars::AbstractArray{Tres},
                  nrm::AbstractArray{Tmsh}, nrm_bar::AbstractArray{Tmsh},
                  F_bar::AbstractArray{Tres}) where {Tmsh, Tsol, Tres}

  kernel = params.entropy_lf_kernel
  applyEntropyKernel_diagE_revm(params, kernel, qL, qR, aux_vars, nrm, nrm_bar,
                                F_bar)
  return nothing
end

mutable struct LFPenalty_revq <: FluxType_revq
end

function (obj::LFPenalty_revq)(params::ParamType,
                      qL::AbstractArray{Tsol,1}, qL_bar::AbstractArray{Tsol, 1},
                      qR::AbstractArray{Tsol, 1}, qR_bar::AbstractArray{Tsol, 1},
                      aux_vars::AbstractArray{Tres}, nrm::AbstractArray{Tmsh},
                      F_bar::AbstractArray{Tres}) where {Tmsh, Tsol, Tres}

  kernel = params.entropy_lf_kernel
  applyEntropyKernel_diagE_revq(params, kernel, qL, qL_bar, qR, qR_bar,
                                aux_vars, nrm, F_bar)
  return nothing
end

"""
  Calls [`calcHLLFlux`](@ref)
"""
mutable struct HLLFlux <: FluxType
end

function (obj::HLLFlux)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractVector,
              F::AbstractVector{Tres}) where {Tsol, Tres}

  calcHLLFlux(params, uL, uR, aux_vars, nrm, F)

  return nothing
end

mutable struct HLLFlux_revq <: FluxType_revq
end

function (obj::HLLFlux_revq)(params::ParamType,
                      qL::AbstractArray{Tsol,1}, qL_bar::AbstractArray{Tsol, 1},
                      qR::AbstractArray{Tsol, 1}, qR_bar::AbstractArray{Tsol, 1},
                      aux_vars::AbstractArray{Tres}, nrm::AbstractArray{Tmsh},
                      F_bar::AbstractArray{Tres}) where {Tmsh, Tsol, Tres}

  calcHLLFlux_revq(params, qL, qL_bar, qR, qR_bar, aux_vars, nrm, F_bar)

  return nothing
end

mutable struct HLLFlux_revm <: FluxType_revm
end

function (obj::HLLFlux_revm)(params::ParamType{Tdim},
                  qL::AbstractArray{Tsol,1}, qR::AbstractArray{Tsol, 1},
                  aux_vars::AbstractArray{Tres},
                  nrm::AbstractArray{Tmsh}, nrm_bar::AbstractArray{Tmsh},
                  F_bar::AbstractArray{Tres}) where {Tdim, Tmsh, Tsol, Tres}

  calcHLLFlux_revm(params, qL, qR, aux_vars, nrm, nrm_bar, F_bar)

  return nothing
end

"""
  A special flux function for computing the face integral of (psi_L - psi_R),
  where psi is the potential flux.  F[1] = psiL - psiR, all other components
  zero.
"""
mutable struct PotentialFlux <: FluxType
end

function (obj::PotentialFlux)(params::ParamType{Tdim},
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractVector,
              F::AbstractVector{Tres}) where {Tdim, Tsol, Tres}

  # psi = rho*U_n
  psiL = zero(Tres); psiR = zero(Tres)
  for i=1:Tdim
    psiL += uL[i+1]*nrm[i]
    psiR += uR[i+1]*nrm[i]
  end

  fill!(F, 0)
  F[1] = psiL - psiR

  return nothing
end

mutable struct PotentialFlux_revq <: FluxType_revq
end

function (obj::PotentialFlux_revq)(params::ParamType{Tdim},
                      qL::AbstractArray{Tsol,1}, qL_bar::AbstractArray{Tsol, 1},
                      qR::AbstractArray{Tsol, 1}, qR_bar::AbstractArray{Tsol, 1},
                      aux_vars::AbstractArray{Tres}, nrm::AbstractArray{Tmsh},
                      F_bar::AbstractArray{Tres}) where {Tdim, Tmsh, Tsol, Tres}

  psiL_bar = F_bar[1]; psiR_bar = -F_bar[1]

  for i=1:Tdim
    qL_bar[i+1] += psiL_bar*nrm[i]
    qR_bar[i+1] += psiR_bar*nrm[i]
  end

  return nothing
end

mutable struct PotentialFlux_revm <: FluxType_revm
end

function (obj::PotentialFlux_revm)(params::ParamType{Tdim},
                  qL::AbstractArray{Tsol,1}, qR::AbstractArray{Tsol, 1},
                  aux_vars::AbstractArray{Tres},
                  nrm::AbstractArray{Tmsh}, nrm_bar::AbstractArray{Tmsh},
                  F_bar::AbstractArray{Tres}) where {Tdim, Tmsh, Tsol, Tres}

  psiL_bar = F_bar[1]; psiR_bar = -F_bar[1]

  for i=1:Tdim
    nrm_bar[i] += qL[i+1]*psiL_bar + qR[i+1]*psiR_bar
  end

  return nothing
end

@doc """
### EulerEquationMod.FluxDict

  This dictonary maps the names of the fluxes (Strings) to the
  functor object itself.  All flux functors should be added to the dictionary.

  All fluxes have one method that calculates the flux in a particular direction
  at a node.  Some fluxes have an additional method that computes the flux
  in several directions at a node in a single function call, which can be
  more efficient.  See calcEulerFlux_standard for an example.

  In general, these functors call similarly-named function in bc_solvers.jl.
  It is recommened to look at the documentation for those functions.

  Each functor must be callable as:

  ```
  function (obj::FunctorName)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractVector,
              F::AbstractVector{Tres}) where {Tsol, Tres}
  ```

  **Inputs**

   * params: ParamType
   * uL: vector containing the left state, length `numDofPerNode`
   * ur: vector containing the right state, length `numDofPerNode`
   * aux_vars: vector of auxiliary variables for the left state
   * nrm: scaled normal vector in xy space, length `mesh.dim`

  **Inputs/Outputs**

   * F: vector to overwrite with the flux, length `numDofPerNode`
  TODO: document signature of the functors here

"""->
global const FluxDict = Dict{String, FluxType}(
"ErrorFlux" => ErrorFlux(),
"ZeroFlux" => ZeroFlux(),
"IdentityFlux" => IdentityFlux(),
"RoeFlux" => RoeFlux(),
"LFFlux" => LFFlux(),
"StandardFlux" => StandardFlux(),
"DucrosFlux" => DucrosFlux(),
"IRFlux" => IRFlux(),
"IRSLFFlux" => IRSLFFlux(),
"LFPenalty" => LFPenalty(),
"HLLFlux" => HLLFlux(),
"PotentialFlux" => PotentialFlux(),
)

@doc """
### EulerEquationMod.getFluxFunctors

  This function retrieves the flux functors from the dictonary and
  stores them to eqn.flux_func.

  Inputs:
    mesh: an AbstractDGMesh
    sbp
    eqn
    opts
"""->
function getFluxFunctors(mesh::AbstractDGMesh, sbp, eqn, opts)

  name = opts["Flux_name"]
  eqn.flux_func = FluxDict[name]
  name = opts["Volume_flux_name"]
  eqn.volume_flux_func = FluxDict[name]

  assertFieldsConcrete(eqn.flux_func)
  assertFieldsConcrete(eqn.volume_flux_func)
  return nothing
end

"""
  Dictionary to put all the reverse mode wrt metric flux functors in.

  Each functor should be callable as:

  ```
  function (obj::FunctorName_revm)(params::ParamType,
                  qL::AbstractArray{Tsol,1}, qR::AbstractArray{Tsol, 1},
                  aux_vars::AbstractArray{Tres},
                  nrm::AbstractArray{Tmsh}, nrm_bar::AbstractArray{Tmsh},
                  F_bar::AbstractArray{Tres}) where {Tmsh, Tsol, Tres}
  ```

  The arguments are similar to the regular functor, the only differences are:

  * F_bar: input vector, same length as `F`
  * nrm_bar: vector to be updated with the back-propigation of `F_bar`, same
             length as `nrm`

"""
global const FluxDict_revm = Dict{String, FluxType_revm}(
"ErrorFlux" => ErrorFlux_revm(),
"ZeroFlux" => ZeroFlux_revm(),
"RoeFlux" => RoeFlux_revm(),
"LFFlux" => LFFlux_revm(),
"IRFlux" => IRFlux_revm(),
"IRSLFFlux" => IRSLFFlux_revm(),
"LFPenalty" => LFPenalty_revm(),
"HLLFlux" => HLLFlux_revm(),
"PotentialFlux" => PotentialFlux_revm(),
)

"""
  Populates the fields of the eqn object that need reverse mode flux functors.
  If opts["need_adjoint"] is false, the error flux functor is used instead.
"""
function getFluxFunctors_revm(mesh::AbstractDGMesh, sbp, eqn, opts)

  if opts["need_adjoint"]
    name = opts["Flux_name"]
  else
    name = "ErrorFlux"
  end
  eqn.flux_func_revm = FluxDict_revm[name]

  if opts["need_adjoint"]
    name = opts["Volume_flux_name"]
  else
    name = "ErrorFlux"
  end
  eqn.volume_flux_func_revm = FluxDict_revm[name]

  assertFieldsConcrete(eqn.flux_func_revm)
  assertFieldsConcrete(eqn.volume_flux_func_revm)

  return nothing
end # End function getFluxFunctors_revm

"""
  Dictionary to put all the reverse mode wrt the solution functors

  Each functor should be callable as

  ```
  function (obj::IRFlux_revq)(params::ParamType,
                      qL::AbstractArray{Tsol,1}, qL_bar::AbstractArray{Tsol, 1},
                      qR::AbstractArray{Tsol, 1}, qR_bar::AbstractArray{Tsol, 1},
                      aux_vars::AbstractArray{Tres}, dir::AbstractArray{Tmsh},
                      F_bar::AbstractArray{Tres}) where {Tmsh, Tsol, Tres}
  ```

  The arguments are similar to the regular functors, the only differences are:

   * F_bar: input vector, same length as `F`
   * qL_bar: vector corresponding to `qL`, to be updated with the back
             propigation of `F_bar`
   * qR_bar: vector corresponding to `qR`, to be updated with back
             propigation of `F_bar`
"""
global const FluxDict_revq = Dict{String, FluxType_revq}(
"ErrorFlux" => ErrorFlux_revq(),
"ZeroFlux" => ZeroFlux_revq(),
"RoeFlux" => RoeFlux_revq(),
"LFFlux" => LFFlux_revq(),
"IRFlux" => IRFlux_revq(),
"IRSLFFlux" => IRSLFFlux_revq(),
"LFPenalty" => LFPenalty_revq(),
"HLLFlux" => HLLFlux_revq(),
"PotentialFlux" => PotentialFlux_revq(),
)

"""
  Populates the fields of [`EulerData`](@ref) that take [`FluxType_revq`](@ref)s

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
"""
function getFluxFunctors_revq(mesh::AbstractDGMesh, sbp, eqn, opts)

  if opts["need_adjoint"]
    name = opts["Flux_name"]
  else
    name = "ErrorFlux"
  end
  eqn.flux_func_revq = FluxDict_revq[name]

  if opts["need_adjoint"]
    name = opts["Volume_flux_name"]
  else
    name = "ErrorFlux"
  end

  eqn.volume_flux_func_revq = FluxDict_revq[name]

  assertFieldsConcrete(eqn.flux_func_revq)
  assertFieldsConcrete(eqn.volume_flux_func_revq)

  return nothing
end # End function getFluxFunctors_revq

import PDESolver.setFluxFunction

function setFluxFunction(mesh::AbstractDGMesh, sbp::AbstractOperator,
                         eqn::EulerData, opts,
                         name::String=opts["Flux_name"])

  eqn.flux_func = FluxDict[name]
  if !opts["need_adjoint"]
    name = "ErrorFlux"
  end
  eqn.flux_func_revm = FluxDict_revm[name]
  eqn.flux_func_revq = FluxDict_revq[name]

  assertFieldsConcrete(eqn.flux_func)
  assertFieldsConcrete(eqn.flux_func_revm)
  assertFieldsConcrete(eqn.flux_func_revq)

  return nothing
end