text:
### A Pluto.jl notebook ###
# v0.18.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

# ╔═╡ 864e1180-f693-11ea-080e-a7d5aabc9ca5
begin
	using Colors, ColorVectorSpace, ImageShow, FileIO, ImageIO
	using ImageShow.ImageCore
	using ColorSchemes

	using InteractiveUtils, PlutoUI
	using LinearAlgebra, SparseArrays, Statistics
end

# ╔═╡ 0db6ee04-81b7-11eb-330c-11b578b72c90
html"""

Section 1.9

 Taking Advantage of Structure

Lecture Video

"""

# ╔═╡ ca1a1072-81b6-11eb-1fee-e7df687cc314
PlutoUI.TableOfContents(aside = true)

# ╔═╡ fe2028ba-f6dc-11ea-0228-938a81a91ace
myonehatvector = [0, 1, 0, 0, 0, 0]

# ╔═╡ 0a902426-f6dd-11ea-0ae4-fb0c47863fe7
# also one "cold"
1 .- myonehatvector

# ╔═╡ 4624cd26-f5d3-11ea-1cf8-7d6555eae5fa
struct OneHot <: AbstractVector{Int}
	n::Int
	k::Int
end

# ╔═╡ 397ac764-f5fe-11ea-20cc-8d7cab19d410
Base.size(x::OneHot) = (x.n,)

# ╔═╡ 82c7046c-f5d3-11ea-04e2-ef7c0f4db5da
Base.getindex(x::OneHot, i::Int) = Int(x.k == i)

# ╔═╡ b0ba5b8c-f5d1-11ea-1304-3f0e47f935fe
md"# Examples of structure"

# ╔═╡ 261c4df2-f5d2-11ea-2c72-7d4b09c46098
md"""
# One-hot vectors
This example comes from machine learning.
"""

# ╔═╡ 3cada3a0-81cc-11eb-04c8-bde26d36a84e
md"""
A **one-hot** vector has a single "hot" element, i.e. a single 1 in a sea of zeros. For example:
"""

# ╔═╡ 8d2c6910-f5d4-11ea-1928-1baf09815687
md"""How much "information" (numbers) do you need to represent a one-hot vector? Is it $n$ or is it two?
"""

# ╔═╡ 54649792-81cc-11eb-1038-9161a4037acf
md"""
(There are also "1-cold" vectors:)
"""

# ╔═╡ 4794e860-81b7-11eb-2c91-8561c20f308a
md"""
## Julia: `structs` (creating a new type in Julia)
"""

# ╔═╡ 67827da8-81cc-11eb-300e-278104d2d958
md"""
We can create our own new types in Julia. Let's create a new type to represent one-hot vectors. It will be a subtype of `AbstractVector`, meaning that it *behaves like* a vector.
"""

# ╔═╡ 9bdabef8-81cc-11eb-14a1-67a9a7d968c0
md"""
We need to specify how long the vector is:
"""

# ╔═╡ a22dcd2c-81cc-11eb-1252-13ace134192d
md"""
and how to extract the $i$th component:
"""

# ╔═╡ b024c318-81cc-11eb-018c-e1f7830ff51b
md"""
(Note that `x.k == i` returns a Boolean value, `true` or `false`, which we are then converting to an `Int`.)
"""

# ╔═╡ 93bfe3ac-f756-11ea-20fb-8f7d586b42f3
myonehotvector = OneHot(6,2)

# ╔═╡ 175039aa-f758-11ea-251a-5db57d7c4b32
myonehotvector[3]

# ╔═╡ c2a4b0a2-81cc-11eb-37a7-db601a6ddfdf
myonehotvector[2]

# ╔═╡ c5ed7d3e-81cc-11eb-3386-15b72db8155d
md"""
This behaves as if it were the original vector, but we are storing only 2 integers.
This is an example of taking advantage of structure.
"""

# ╔═╡ e2e354a8-81b7-11eb-311a-35151063c2a7
md"""
## Julia: dump and Dump
"""

# ╔═╡ dc5a96ba-81cc-11eb-3189-25920df48afa
md"""
`dump` shows the internal data stored inside a given object:
"""

# ╔═╡ af0d3c22-f756-11ea-37d6-11b630d2314a
with_terminal() do
	dump(myonehotvector)
end

# ╔═╡ 06e3a842-26b8-4417-9cf5-8a083ccdb264
md"""
Since `dump` writes to a terminal, you can't use it directly inside Pluto, which is why we wrap it inside a `with_terminal` block. You can also use the function `Dump` from `PlutoUI`:
"""

# ╔═╡ 91172a3e-81b7-11eb-0953-9f5e0207f863
Dump(myonehotvector)

# ╔═╡ 4bbf3f58-f788-11ea-0d24-6b0fb070829e
myonehotvector

# ╔═╡ fe70d104-81b7-11eb-14d0-eb5237d8ea6c
md"""
### Visualizing a one-hot vector
"""

# ╔═╡ ef8f44b2-f5fc-11ea-1e4d-bd873cd39d6c
md"""
n=$(@bind nn Slider(1:20, show_value=true))
"""

# ╔═╡ fd9211c0-f5fc-11ea-1745-7f2dae88af9e
md"""
k=$(@bind kk Slider(1:nn, default=1, show_value=true))
"""

# ╔═╡ f1154df8-f693-11ea-3b16-f32835fcc470
x = OneHot(nn, kk)

# ╔═╡ 81c35324-f5d4-11ea-2338-9f982d38732c
md"# Diagonal matrices"

# ╔═╡ 2cfda0dc-f5d5-11ea-16c4-b5a33b90e37f
md"Another example is diagonal matrices. Here's how you might see them in high school:"

# ╔═╡ 150432d4-f5d5-11ea-32b2-19a2a91d9637
denseD = [5 0 0
	    0 6 0
	    0 0 -10]

# ╔═╡ 44215aa4-f695-11ea-260e-b564c6fbcd4a
md"Julia has a better way of representing them:"

# ╔═╡ 21328d1c-f5d5-11ea-288e-4171ad35326d
D = Diagonal(denseD)

# ╔═╡ 75761cc0-81cd-11eb-1186-7d47debd68ca
md"""
It even displays nicely, with dots instead of zeros.

We can also create a diagonal matrix from the values on the diagonal:
"""

# ╔═╡ 6bd8a886-f758-11ea-2587-870a3fa9d710
Diagonal([5, 6, -10])

# ╔═╡ 4c533ac6-f695-11ea-3724-b955eaaeee49
md"How much information is stored for each representation? We can use Julia's `dump` function to find out:"

# ╔═╡ 466901ea-f5d5-11ea-1db5-abf82c96eabf
Dump(denseD)

# ╔═╡ b38c4aae-f5d5-11ea-39b6-7b0c7d529019
Dump(D)

# ╔═╡ 93e04ed8-81cd-11eb-214a-a761ef8c406f
md"""
We see that `Diagonal` stores only the diagonal entries, not the zeros!
"""

# ╔═╡ e90c55fc-f5d5-11ea-10f1-470ff772985d
md"""We should always look for *structure* where it exists!"""

# ╔═╡ 19775c3c-f5d6-11ea-15c2-89618e654a1e
md"# Sparse matrices"

# ╔═╡ 653792a8-f695-11ea-1ae0-43761c502583
md"A *sparse* matrix is a matrix that has many zeros, and is hence worth storing in a *sparse* representation:"

# ╔═╡ 79c94d2a-f75a-11ea-031d-09d70d229e15
denseM = [0 0 9; 0 0 0;12 0 4]

# ╔═╡ 10bc5d50-81b9-11eb-2ac7-354a6c6c826b
md"""
The above displays a sparse matrix in so-called `(i, j, value)` form.  We could
store sparse matrices in this way:
"""

# ╔═╡ 77d6a952-81ba-11eb-24e3-cb6510a59455
M = sparse(denseM)

# ╔═╡ 1f3ba55a-81b9-11eb-001f-593b9d8639ca
md"""
Although it looks like it's stored like this, in fact the actual storage format is different. In the Julia `SparseArrays.jl` package, the storage format
is **compressed sparse column** (CSC) format,
which is generally considered favorable for arithmetic, matrix-vector
	products and column slicing.  Of course, for specific matrices, other
		formats might be better.

* `nzval` contains the nonzero matrix entries
* `rowval` is the "i" or row entry for the corresponding value in nzval
  * length(rowval) == length(nzval)
* `colptr[j]` points into nzval and tells you the first nonzero in or after column j
* The last entry of colptr points beyond the end of nzval to indicate no more columns.
  * length(colptr) == number of columns + 1
"""

# ╔═╡ 3d4a702e-f75a-11ea-031c-333d591fc442
Dump(sparse(M))

# ╔═╡ 80ff4010-81bb-11eb-374e-215a57defb0b
md"""
 An example where CSC may not be a great choice is the following. The reason is that  `colptr` must have an entry in each column:
"""

# ╔═╡ 5de72b7c-f5d6-11ea-1b6f-35b830b5fb34
M2 = sparse([1, 2, 10^6], [4, 9, 10^6], [7, 8, 9])

# ╔═╡ 8b60629e-f5d6-11ea-27c8-d934460d3a57
with_terminal() do
	dump(M2)
end

# ╔═╡ 2fd7e52e-f5d7-11ea-3b5a-1f338e2451e0
M3 = [1 0 2 0 10; 0 3 4 0 9; 0 0 0 5 8; 0 0 0 0 7]

# ╔═╡ 2e87d4fe-81bc-11eb-0d16-b988bcedcc73
M4 = M3 .* 0

# ╔═╡ cde79f38-f5d6-11ea-3297-0b5b240f7b9e
Dump(sparse(M4))

# ╔═╡ aa09c008-f5d8-11ea-1bdc-b51ee6eb2478
sparse(M4)

# ╔═╡ 62a6ec62-f5d9-11ea-071e-ed33c5dea0cd
md"""# Random vectors
"""

# ╔═╡ 7f63daf6-f695-11ea-0b80-8702a83103a4
md"How much structure is there in a *random* vector?"

# ╔═╡ 67274c3c-f5d9-11ea-3475-c9d228e3bd5a
v = rand(1:9, 1_000_000)

# ╔═╡ 765c6552-f5d9-11ea-29d3-bfe7b4b04612
md"""You might guess that there is "no structure". But you can actually think of randomness itself as a structure.

For example, take the mean and standard deviation -- some would say that's the structure.

"""

# ╔═╡ 126fb3ea-f5da-11ea-2f7d-0b3259a296ce
mean(v), std(v), 5, sqrt(10 * 2/3)

# ╔═╡ ed0b2358-81ce-11eb-3339-93abcc06fd91
md"""
If you repeat the calculation, to 3 or 4 digits the mean and standard deviation don't change, and are approximately equal to the theoretical values on the right.
"""

# ╔═╡ 24ce92fa-81cf-11eb-30f0-b1e357d79d53
md"""
We can also count how many times each digit occurs in the data set:
"""

# ╔═╡ 2d4500e0-81cf-11eb-1699-d310074fddf5
[sum(v .== i) for i in 1:9]

# ╔═╡ 3546ff30-81cf-11eb-3afc-05c5db61366f
md"""
We see that each number occurs roughly the same number of times.
"""

# ╔═╡ 9b9e2c2a-f5da-11ea-369b-b513b196515b
md"Statisticians (and professors who've just graded exams) might say that under certain circumstances the mean and the variance give you the necessary structure, and the rest can be thrown away."

# ╔═╡ e68b98ea-f5da-11ea-1a9d-db45e4f80241
m = sum(v) / length(v)  # mean

# ╔═╡ f20ccac4-f5da-11ea-0e69-413b5e49f423
σ² = sum( (v .- m) .^ 2 ) / (length(v) - 1)

# ╔═╡ 12a2e96c-f5db-11ea-1c3e-494ae7446886
σ = sqrt(σ²)

# ╔═╡ 22487ce2-f5db-11ea-32e9-6f70ab2c0353
std(v)

# ╔═╡ 389ae62e-f5db-11ea-1557-c3adbbee0e5c
md"Sometimes the summary statistics are all you want. (But sometimes not.)"

# ╔═╡ 0c2b6408-f5d9-11ea-2b7f-7fece2eecc1f
md"# Multiplication tables"

# ╔═╡ 5d767290-f5dd-11ea-2189-81198fd216ce
outer(v, w) = [x * y for x ∈ v, y ∈ w]  # just a multiplication table

# ╔═╡ 587790ce-f6de-11ea-12d9-fde2a17ae314
outer(1:10, 1:10)

# ╔═╡ a39e8256-f6de-11ea-3170-c923b56609da
md"Did you memorize this in third grade?"

# ╔═╡ 8c84edd0-f6de-11ea-2180-61c6b81aac3b
@bind k Slider(1:14, show_value=true)

# ╔═╡ 22b73baa-f6df-11ea-197f-bbb4bd1a7ef5
outer(1:k, 1:k)

# ╔═╡ b2332814-f6e6-11ea-1c7d-556c7d4687f1
outer([2,4,6], [10,100,1000])

# ╔═╡ 9ab7a72e-81cf-11eb-2b78-073ff51cae58
md"""
A multiplication table is clearly a structure, but it's not sparse -- there are no zeros. Nonetheless you need much less information to reconstruct the matrix.
"""

# ╔═╡ fd8dd108-f6df-11ea-2f7c-3d99d054ac15
md"In the context of 1:k times 1:k, just one number k is needed."

# ╔═╡ 165788b2-f601-11ea-3e69-cdbbb6558e54
md"If you look at the following matrix? Does it have any structure? It's certainly more hidden."

# ╔═╡ 22941bb8-f601-11ea-1d6e-0d955297bc2e
outer( rand(3), rand(4) )  # but it's just a multiplication table

# ╔═╡ c33bf00e-81cf-11eb-1e1a-e5a879a45093
md"""
You might guess by visualizing the matrix that it is a multiplication table:
"""

# ╔═╡ 7ff664f0-f74b-11ea-0d2d-b53f19e4f4bf
md"We can factor out a multiplication table, if it's there:"

# ╔═╡ a0611eaa-81bc-11eb-1d23-c12ab14138b1
md"""
### Julia: Exceptions are thrown (generated) using `error`

An exception is anything that can interrupt a program, e.g. invalid input data.
"""

# ╔═╡ a4728944-f74b-11ea-03c3-9123908c1f8e
function factor( mult_table )
	v = mult_table[:, 1]
	w = mult_table[1, :]

	if v[1] ≠ 0
		w /= v[1]
	end

	# Good code has a check:
	if outer(v, w) ≈ mult_table
	   return v, w
	else
		error("Input is not a multiplication table")
	end
end

# ╔═╡ 05c35402-f752-11ea-2708-59cf5ef74fb4
factor( outer([1, 2, 3], [2, 2, 2] ) )

# ╔═╡ 8c11b19e-81bc-11eb-184b-bf6ffefe29de
md"""
A random 2x2 matrix is not a multiplication table. Most matrices are not given by multiplication tables.
"""

# ╔═╡ 8baaa994-f752-11ea-18d9-b3d0a6b9f7d9
factor( rand(2,2) )

# ╔═╡ d92bece4-f754-11ea-0242-99f198bb5b7b
md" Let's add two (or more) multiplication tables:"

# ╔═╡ e740999c-f754-11ea-2089-4b7a9aec6030
A = sum( outer(rand(3),rand(3)) for i=1:2 )

# ╔═╡ 0a79a7b4-f755-11ea-1b2d-21173567b257
md"Is it possible, given the matrix, to find the structure? E.g. to show that a matrix is a sum of outer products (multiplication table)."

# ╔═╡ 5adb98c2-f6e0-11ea-1fde-53b0fd6639c3
md"The answer is yes: The **Singular-Value Decomposition** (SVD) from algebra can find the structure!"

# ╔═╡ 487d6f9c-81d0-11eb-3bb0-336a4beb9b38
md"""
Let's take the SVD and calculate the sum of two outer products:
"""

# ╔═╡ 5a493052-f601-11ea-2f5f-f940412905f2
begin
	U, Σ, V = svd(A)

    outer( U[:, 1], V[:, 1] * Σ[1] ) + outer( U[:, 2], V[:, 2] * Σ[2] )
end

# ╔═╡ 55b76aee-81d0-11eb-0bcc-413f5bd14360
md"""
We see that we reconstruct the original matrix!"
"""

# ╔═╡ 709bf30a-f755-11ea-2e82-bd511e598c77
B = rand(3,3)

# ╔═╡ 782532b0-f755-11ea-1385-cd1a28c4b9d5
begin
	UU, ΣΣ, VV = svd( B )
    outer( UU[:,1], VV[:,1] * ΣΣ[1] ) + outer( UU[:,2], VV[:,2] * ΣΣ[2] )
end

# ╔═╡ 5bc4ab0a-f755-11ea-0fad-4987ad9fc02f
md"and it can approximate too!"

# ╔═╡ a5d637ea-f5de-11ea-3b70-877e876bc9c9
flag = outer([1,1,1,2,2,2,1,1,1], [1,1,1,1,1,1,1,1,1])

# ╔═╡ 21bbb60a-f5df-11ea-2c1b-dd716a657df8
cs = distinguishable_colors(100)

# ╔═╡ 2668e100-f5df-11ea-12b0-073a578a5edb
cs[flag]

# ╔═╡ e8d727f2-f5de-11ea-1456-f72602e81e0d
cs[flag + flag']

# ╔═╡ f5fcdeea-f75c-11ea-1fc3-731f0ef1ad14
outer([1,1,1,2,2,2,1,1,1], [1,1,1,1,1,1,1,1,1]) + outer([1,1,1,1,1,1,1,1,1], [1,1,1,2,2,2,1,1,1])

# ╔═╡ 0373fbf6-f75d-11ea-2a9e-cbb714d69cf4
cs[outer([1,1,1,2,2,2,1,1,1], [1,1,1,1,1,1,1,1,1]) + outer([1,1,1,1,1,1,1,1,1], [1,1,1,2,2,2,1,1,1])]

# ╔═╡ ebd72fb8-f5e0-11ea-0630-573337dff753
md"""
# Singular Value Decomposition (SVD): A tool to find structure
"""

# ╔═╡ b6478e1a-f5f6-11ea-3b92-6d4f067285f4
tree_url = "https://user-images.githubusercontent.com/6933510/110924885-d7f1b200-8322-11eb-9df7-7abf29c8db7d.png"

# ╔═╡ f2c11f88-f5f8-11ea-3e02-c1d4fa22031e
image = load(download(tree_url))

# ╔═╡ 29062f7a-f5f9-11ea-2682-1374e7694e32
picture = Float64.(channelview(image));

# ╔═╡ 5471fd30-f6e2-11ea-2cd7-7bd48c42db99
size(picture)

# ╔═╡ 6156fd1e-f5f9-11ea-06a9-211c7ab813a4
pr, pg, pb = eachslice(picture, dims=1)

# ╔═╡ a9766e68-f5f9-11ea-0019-6f9d02050521
[RGB.(pr, 0, 0) RGB.(0, pg, 0) RGB.(0, 0, pb)]

# ╔═╡ 0c0ee362-f5f9-11ea-0f75-2d2810c88d65
begin
	Ur, Σr, Vr = svd(pr)
	Ug, Σg, Vg = svd(pg)
	Ub, Σb, Vb = svd(pb)
end;

# ╔═╡ b95ce51a-f632-11ea-3a64-f7c218b9b3c9
@bind n Slider(1:200, show_value=true)

# ╔═╡ 7ba6e6a6-f5fa-11ea-2bcd-616d5a3c898b
RGB.(sum(outer(Ur[:,i], Vr[:,i]) .* Σr[i] for i in 1:n),
	 sum(outer(Ug[:,i], Vg[:,i]) .* Σg[i] for i in 1:n),
	 sum(outer(Ub[:,i], Vb[:,i]) .* Σb[i] for i in 1:n))

# ╔═╡ 8df84fcc-f5d5-11ea-312f-bf2a3b3ce2ce
md"# Appendix"

# ╔═╡ 0edd7cca-834f-11eb-0232-ff0850027f76
md"## Syntax Learned"

# ╔═╡ 69be8194-81b7-11eb-0452-0bc8b9f22286
md"""
Syntax to be learned:

* A `struct` is a great way to embody structure.
* `dump` and `Dump`: to see what's inside a data structure.
* `Diagonal`, `sparse`
* `error` (throws an exception)
* `svd` (Singular Value Decomposition)
"""

# ╔═╡ 1c462f68-834f-11eb-1447-85848814769b
[ dump, Diagonal, error, svd]

# ╔═╡ 5813e1b2-f5ff-11ea-2849-a1def74fc065
begin
	show_image(M) = get.([ColorSchemes.rainbow], M ./ maximum(M))
	show_image(x::AbstractVector) = show_image(x')
end

# ╔═╡ 982590d4-f5ff-11ea-3802-73292c75ad6c
show_image(x)

# ╔═╡ 2f75df7e-f601-11ea-2fc2-aff4f335af33
show_image( outer( rand(10), rand(10) ))

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
ColorVectorSpace = "c3611d14-8923-5661-9e6a-0046d554d3a4"
Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19"
ImageShow = "4e3cecfd-b093-5904-9786-8bbb286a6a31"
InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[compat]
ColorSchemes = "~3.17.1"
ColorVectorSpace = "~0.9.8"
Colors = "~0.12.8"
FileIO = "~1.13.0"
ImageIO = "~0.6.1"
ImageShow = "~0.3.3"
PlutoUI = "~0.7.37"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.AbstractFFTs]]
deps = ["ChainRulesCore", "LinearAlgebra"]
git-tree-sha1 = "6f1d9bc1c08f9f4a8fa92e3ea3cb50153a1b40d4"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.1.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.3"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.CEnum]]
git-tree-sha1 = "215a9aa4a1f23fbd05b92769fdd62559488d70e9"
uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
version = "0.4.1"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "c9a6160317d1abe9c44b3beb367fd448117679ca"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.13.0"

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

[[deps.ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "12fc73e5e0af68ad3137b886e3f7c1eacfca2640"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.17.1"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.ColorVectorSpace]]
deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "TensorCore"]
git-tree-sha1 = "3f1f500312161f1ae067abe07d13b40f78f32e07"
uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4"
version = "0.9.8"

[[deps.Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "96b0bc6c52df76506efc8a441c6cf1adcb1babc4"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.42.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "80ced645013a5dbdc52cf70329399c35ce007fae"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.13.0"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Graphics]]
deps = ["Colors", "LinearAlgebra", "NaNMath"]
git-tree-sha1 = "1c5a84319923bea76fa145d49e93aa4394c73fc2"
uuid = "a2bd30eb-e257-5431-a919-1863eab51364"
version = "1.1.1"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.ImageBase]]
deps = ["ImageCore", "Reexport"]
git-tree-sha1 = "b51bb8cae22c66d0f6357e3bcb6363145ef20835"
uuid = "c817782e-172a-44cc-b673-b171935fbb9e"
version = "0.1.5"

[[deps.ImageCore]]
deps = ["AbstractFFTs", "ColorVectorSpace", "Colors", "FixedPointNumbers", "Graphics", "MappedArrays", "MosaicViews", "OffsetArrays", "PaddedViews", "Reexport"]
git-tree-sha1 = "9a5c62f231e5bba35695a20988fc7cd6de7eeb5a"
uuid = "a09fc81d-aa75-5fe9-8630-4744c3626534"
version = "0.9.3"

[[deps.ImageIO]]
deps = ["FileIO", "JpegTurbo", "Netpbm", "OpenEXR", "PNGFiles", "QOI", "Sixel", "TiffImages", "UUIDs"]
git-tree-sha1 = "464bdef044df52e6436f8c018bea2d48c40bb27b"
uuid = "82e4d734-157c-48bb-816b-45c225c6df19"
version = "0.6.1"

[[deps.ImageShow]]
deps = ["Base64", "FileIO", "ImageBase", "ImageCore", "OffsetArrays", "StackViews"]
git-tree-sha1 = "d0ac64c9bee0aed6fdbb2bc0e5dfa9a3a78e3acc"
uuid = "4e3cecfd-b093-5904-9786-8bbb286a6a31"
version = "0.3.3"

[[deps.Imath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "87f7662e03a649cffa2e05bf19c303e168732d3e"
uuid = "905a6f67-0a94-5f89-b386-d35d92009cd1"
version = "3.1.2+0"

[[deps.IndirectArrays]]
git-tree-sha1 = "012e604e1c7458645cb8b436f8fba789a51b257f"
uuid = "9b13fd28-a010-5f03-acff-a1bbcff69959"
version = "1.0.0"

[[deps.Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "91b5dcf362c5add98049e6c29ee756910b03051d"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.3"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.JpegTurbo]]
deps = ["CEnum", "FileIO", "ImageCore", "JpegTurbo_jll", "TOML"]
git-tree-sha1 = "a77b273f1ddec645d1b7c4fd5fb98c8f90ad10a5"
uuid = "b835a17e-a41a-41e7-81f0-2f016b05efe0"
version = "0.1.1"

[[deps.JpegTurbo_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b53380851c6e6664204efb2e62cd24fa5c47e4ba"
uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
version = "2.1.2+0"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "3f7cb7157ef860c637f3f4929c8ed5d9716933c6"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.7"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MappedArrays]]
git-tree-sha1 = "e8b359ef06ec72e8c030463fe02efe5527ee5142"
uuid = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900"
version = "0.4.1"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MosaicViews]]
deps = ["MappedArrays", "OffsetArrays", "PaddedViews", "StackViews"]
git-tree-sha1 = "b34e3bc3ca7c94914418637cb10cc4d1d80d877d"
uuid = "e94cdb99-869f-56ef-bcf0-1ae2bcbe0389"
version = "0.3.3"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NaNMath]]
git-tree-sha1 = "b086b7ea07f8e38cf122f5016af580881ac914fe"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.7"

[[deps.Netpbm]]
deps = ["FileIO", "ImageCore"]
git-tree-sha1 = "18efc06f6ec36a8b801b23f076e3c6ac7c3bf153"
uuid = "f09324ee-3d7c-5217-9330-fc30815ba969"
version = "1.0.2"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OffsetArrays]]
deps = ["Adapt"]
git-tree-sha1 = "043017e0bdeff61cfbb7afeb558ab29536bbb5ed"
uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
version = "1.10.8"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OpenEXR]]
deps = ["Colors", "FileIO", "OpenEXR_jll"]
git-tree-sha1 = "327f53360fdb54df7ecd01e96ef1983536d1e633"
uuid = "52e1d378-f018-4a11-a4be-720524705ac7"
version = "0.3.2"

[[deps.OpenEXR_jll]]
deps = ["Artifacts", "Imath_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "923319661e9a22712f24596ce81c54fc0366f304"
uuid = "18a262bb-aa17-5467-a713-aee519bc75cb"
version = "3.1.1+0"

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.PNGFiles]]
deps = ["Base64", "CEnum", "ImageCore", "IndirectArrays", "OffsetArrays", "libpng_jll"]
git-tree-sha1 = "eb4dbb8139f6125471aa3da98fb70f02dc58e49c"
uuid = "f57f5aa1-a3ce-4bc8-8ab9-96f992907883"
version = "0.3.14"

[[deps.PaddedViews]]
deps = ["OffsetArrays"]
git-tree-sha1 = "03a7a85b76381a3d04c7a1656039197e70eda03d"
uuid = "5432bcbf-9aad-5242-b902-cca2824c8663"
version = "0.5.11"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.3"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PkgVersion]]
deps = ["Pkg"]
git-tree-sha1 = "a7a7e1a88853564e551e4eba8650f8c38df79b37"
uuid = "eebad327-c553-4316-9ea0-9fa01ccd7688"
version = "0.1.1"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "bf0a1121af131d9974241ba53f601211e9303a9e"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.37"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "d3538e7f8a790dc8903519090857ef8e1283eecd"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.5"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.ProgressMeter]]
deps = ["Distributed", "Printf"]
git-tree-sha1 = "afadeba63d90ff223a6a48d2009434ecee2ec9e8"
uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
version = "1.7.1"

[[deps.QOI]]
deps = ["ColorTypes", "FileIO", "FixedPointNumbers"]
git-tree-sha1 = "18e8f4d1426e965c7b532ddd260599e1510d26ce"
uuid = "4b34888f-f399-49d4-9bb3-47ed5cae4e65"
version = "1.0.0"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.Sixel]]
deps = ["Dates", "FileIO", "ImageCore", "IndirectArrays", "OffsetArrays", "REPL", "libsixel_jll"]
git-tree-sha1 = "8fb59825be681d451c246a795117f317ecbcaa28"
uuid = "45858cf5-a6b0-47a3-bbea-62219f50df47"
version = "0.1.2"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.1.4"

[[deps.StackViews]]
deps = ["OffsetArrays"]
git-tree-sha1 = "46e589465204cd0c08b4bd97385e4fa79a0c770c"
uuid = "cae243ae-269e-4f55-b966-ac2d0dc13c15"
version = "0.1.1"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.TensorCore]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "1feb45f88d133a655e001435632f019a9a1bcdb6"
uuid = "62fd8b95-f654-4bbd-a8a5-9c27f68ccd50"
version = "0.1.1"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.TiffImages]]
deps = ["ColorTypes", "DataStructures", "DocStringExtensions", "FileIO", "FixedPointNumbers", "IndirectArrays", "Inflate", "OffsetArrays", "PkgVersion", "ProgressMeter", "UUIDs"]
git-tree-sha1 = "aaa19086bc282630d82f818456bc40b4d314307d"
uuid = "731e570b-9d59-4bfa-96dc-6df516fadf69"
version = "0.5.4"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.libpng_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
version = "1.6.38+0"

[[deps.libsixel_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "78736dab31ae7a53540a6b752efc61f77b304c5b"
uuid = "075b6546-f08a-558a-be8f-8157d0f608a5"
version = "1.8.6+1"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─0db6ee04-81b7-11eb-330c-11b578b72c90
# ╟─ca1a1072-81b6-11eb-1fee-e7df687cc314
# ╟─b0ba5b8c-f5d1-11ea-1304-3f0e47f935fe
# ╠═864e1180-f693-11ea-080e-a7d5aabc9ca5
# ╟─261c4df2-f5d2-11ea-2c72-7d4b09c46098
# ╟─3cada3a0-81cc-11eb-04c8-bde26d36a84e
# ╠═fe2028ba-f6dc-11ea-0228-938a81a91ace
# ╟─8d2c6910-f5d4-11ea-1928-1baf09815687
# ╟─54649792-81cc-11eb-1038-9161a4037acf
# ╠═0a902426-f6dd-11ea-0ae4-fb0c47863fe7
# ╟─4794e860-81b7-11eb-2c91-8561c20f308a
# ╟─67827da8-81cc-11eb-300e-278104d2d958
# ╠═4624cd26-f5d3-11ea-1cf8-7d6555eae5fa
# ╟─9bdabef8-81cc-11eb-14a1-67a9a7d968c0
# ╠═397ac764-f5fe-11ea-20cc-8d7cab19d410
# ╟─a22dcd2c-81cc-11eb-1252-13ace134192d
# ╠═82c7046c-f5d3-11ea-04e2-ef7c0f4db5da
# ╟─b024c318-81cc-11eb-018c-e1f7830ff51b
# ╠═93bfe3ac-f756-11ea-20fb-8f7d586b42f3
# ╠═175039aa-f758-11ea-251a-5db57d7c4b32
# ╠═c2a4b0a2-81cc-11eb-37a7-db601a6ddfdf
# ╟─c5ed7d3e-81cc-11eb-3386-15b72db8155d
# ╟─e2e354a8-81b7-11eb-311a-35151063c2a7
# ╟─dc5a96ba-81cc-11eb-3189-25920df48afa
# ╠═af0d3c22-f756-11ea-37d6-11b630d2314a
# ╟─06e3a842-26b8-4417-9cf5-8a083ccdb264
# ╠═91172a3e-81b7-11eb-0953-9f5e0207f863
# ╠═4bbf3f58-f788-11ea-0d24-6b0fb070829e
# ╟─fe70d104-81b7-11eb-14d0-eb5237d8ea6c
# ╟─ef8f44b2-f5fc-11ea-1e4d-bd873cd39d6c
# ╟─fd9211c0-f5fc-11ea-1745-7f2dae88af9e
# ╠═f1154df8-f693-11ea-3b16-f32835fcc470
# ╠═982590d4-f5ff-11ea-3802-73292c75ad6c
# ╟─81c35324-f5d4-11ea-2338-9f982d38732c
# ╟─2cfda0dc-f5d5-11ea-16c4-b5a33b90e37f
# ╠═150432d4-f5d5-11ea-32b2-19a2a91d9637
# ╟─44215aa4-f695-11ea-260e-b564c6fbcd4a
# ╠═21328d1c-f5d5-11ea-288e-4171ad35326d
# ╟─75761cc0-81cd-11eb-1186-7d47debd68ca
# ╟─6bd8a886-f758-11ea-2587-870a3fa9d710
# ╟─4c533ac6-f695-11ea-3724-b955eaaeee49
# ╠═466901ea-f5d5-11ea-1db5-abf82c96eabf
# ╠═b38c4aae-f5d5-11ea-39b6-7b0c7d529019
# ╟─93e04ed8-81cd-11eb-214a-a761ef8c406f
# ╟─e90c55fc-f5d5-11ea-10f1-470ff772985d
# ╟─19775c3c-f5d6-11ea-15c2-89618e654a1e
# ╟─653792a8-f695-11ea-1ae0-43761c502583
# ╠═79c94d2a-f75a-11ea-031d-09d70d229e15
# ╟─10bc5d50-81b9-11eb-2ac7-354a6c6c826b
# ╠═77d6a952-81ba-11eb-24e3-cb6510a59455
# ╟─1f3ba55a-81b9-11eb-001f-593b9d8639ca
# ╠═3d4a702e-f75a-11ea-031c-333d591fc442
# ╟─80ff4010-81bb-11eb-374e-215a57defb0b
# ╠═5de72b7c-f5d6-11ea-1b6f-35b830b5fb34
# ╠═8b60629e-f5d6-11ea-27c8-d934460d3a57
# ╠═2fd7e52e-f5d7-11ea-3b5a-1f338e2451e0
# ╠═2e87d4fe-81bc-11eb-0d16-b988bcedcc73
# ╠═cde79f38-f5d6-11ea-3297-0b5b240f7b9e
# ╠═aa09c008-f5d8-11ea-1bdc-b51ee6eb2478
# ╟─62a6ec62-f5d9-11ea-071e-ed33c5dea0cd
# ╟─7f63daf6-f695-11ea-0b80-8702a83103a4
# ╠═67274c3c-f5d9-11ea-3475-c9d228e3bd5a
# ╟─765c6552-f5d9-11ea-29d3-bfe7b4b04612
# ╠═126fb3ea-f5da-11ea-2f7d-0b3259a296ce
# ╟─ed0b2358-81ce-11eb-3339-93abcc06fd91
# ╟─24ce92fa-81cf-11eb-30f0-b1e357d79d53
# ╠═2d4500e0-81cf-11eb-1699-d310074fddf5
# ╟─3546ff30-81cf-11eb-3afc-05c5db61366f
# ╟─9b9e2c2a-f5da-11ea-369b-b513b196515b
# ╠═e68b98ea-f5da-11ea-1a9d-db45e4f80241
# ╠═f20ccac4-f5da-11ea-0e69-413b5e49f423
# ╠═12a2e96c-f5db-11ea-1c3e-494ae7446886
# ╠═22487ce2-f5db-11ea-32e9-6f70ab2c0353
# ╟─389ae62e-f5db-11ea-1557-c3adbbee0e5c
# ╟─0c2b6408-f5d9-11ea-2b7f-7fece2eecc1f
# ╠═5d767290-f5dd-11ea-2189-81198fd216ce
# ╠═587790ce-f6de-11ea-12d9-fde2a17ae314
# ╟─a39e8256-f6de-11ea-3170-c923b56609da
# ╠═8c84edd0-f6de-11ea-2180-61c6b81aac3b
# ╠═22b73baa-f6df-11ea-197f-bbb4bd1a7ef5
# ╠═b2332814-f6e6-11ea-1c7d-556c7d4687f1
# ╟─9ab7a72e-81cf-11eb-2b78-073ff51cae58
# ╟─fd8dd108-f6df-11ea-2f7c-3d99d054ac15
# ╟─165788b2-f601-11ea-3e69-cdbbb6558e54
# ╠═22941bb8-f601-11ea-1d6e-0d955297bc2e
# ╟─c33bf00e-81cf-11eb-1e1a-e5a879a45093
# ╠═2f75df7e-f601-11ea-2fc2-aff4f335af33
# ╟─7ff664f0-f74b-11ea-0d2d-b53f19e4f4bf
# ╟─a0611eaa-81bc-11eb-1d23-c12ab14138b1
# ╠═a4728944-f74b-11ea-03c3-9123908c1f8e
# ╠═05c35402-f752-11ea-2708-59cf5ef74fb4
# ╟─8c11b19e-81bc-11eb-184b-bf6ffefe29de
# ╠═8baaa994-f752-11ea-18d9-b3d0a6b9f7d9
# ╟─d92bece4-f754-11ea-0242-99f198bb5b7b
# ╠═e740999c-f754-11ea-2089-4b7a9aec6030
# ╟─0a79a7b4-f755-11ea-1b2d-21173567b257
# ╟─5adb98c2-f6e0-11ea-1fde-53b0fd6639c3
# ╟─487d6f9c-81d0-11eb-3bb0-336a4beb9b38
# ╟─5a493052-f601-11ea-2f5f-f940412905f2
# ╟─55b76aee-81d0-11eb-0bcc-413f5bd14360
# ╠═709bf30a-f755-11ea-2e82-bd511e598c77
# ╠═782532b0-f755-11ea-1385-cd1a28c4b9d5
# ╟─5bc4ab0a-f755-11ea-0fad-4987ad9fc02f
# ╠═a5d637ea-f5de-11ea-3b70-877e876bc9c9
# ╠═21bbb60a-f5df-11ea-2c1b-dd716a657df8
# ╠═2668e100-f5df-11ea-12b0-073a578a5edb
# ╠═e8d727f2-f5de-11ea-1456-f72602e81e0d
# ╠═f5fcdeea-f75c-11ea-1fc3-731f0ef1ad14
# ╠═0373fbf6-f75d-11ea-2a9e-cbb714d69cf4
# ╟─ebd72fb8-f5e0-11ea-0630-573337dff753
# ╠═b6478e1a-f5f6-11ea-3b92-6d4f067285f4
# ╠═f2c11f88-f5f8-11ea-3e02-c1d4fa22031e
# ╠═29062f7a-f5f9-11ea-2682-1374e7694e32
# ╠═5471fd30-f6e2-11ea-2cd7-7bd48c42db99
# ╠═6156fd1e-f5f9-11ea-06a9-211c7ab813a4
# ╠═a9766e68-f5f9-11ea-0019-6f9d02050521
# ╠═0c0ee362-f5f9-11ea-0f75-2d2810c88d65
# ╠═b95ce51a-f632-11ea-3a64-f7c218b9b3c9
# ╠═7ba6e6a6-f5fa-11ea-2bcd-616d5a3c898b
# ╟─8df84fcc-f5d5-11ea-312f-bf2a3b3ce2ce
# ╟─0edd7cca-834f-11eb-0232-ff0850027f76
# ╟─69be8194-81b7-11eb-0452-0bc8b9f22286
# ╠═1c462f68-834f-11eb-1447-85848814769b
# ╟─5813e1b2-f5ff-11ea-2849-a1def74fc065
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using Statistics
using FillArrays
using Rotations
using CoordinateTransformations
using ImageTransformations
using Interpolations
using LinearAlgebra
using StaticArrays

"""
    inject(frame::AbstractMatrix, psf, [angle]; x, y, amp=1, center=center(frame), kwargs...)

Injects the `psf` into `frame` at the given position.

# Examples

`inject` works with two types of PSF models: matrices and synthetic models.

## Matrices

If you pass an `AbstractMatrix` to `inject`, you can optionally specify `degree` and `fill` as keyword arguments. By default, `degree=Interpolations.Linear()` and `fill=0`. These are used to create an `Interpolations.AbstractExtrapolation` type which can be arbitrarily transformed. From here, the `x`, `y`, and `amp` arguments will determine the position and will be optionally rotated `angle` degrees counter-clockwise around the `center`.

```jldoctest
julia> inject(zeros(5, 5), ones(1, 1); x=4, y=3, amp=2)
5×5 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  2.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> inject(zeros(5, 5), ones(3, 3), 90; x=4, y=3, amp=2)
5×5 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  2.0  2.0  2.0
 0.0  0.0  2.0  2.0  2.0
 0.0  0.0  2.0  2.0  2.0
 0.0  0.0  0.0  0.0  0.0
```

## Synthetic Models

The synthetic models from [`PSFModels.jl`](https://github.com/JuliaAstro/PSFModels.jl) can easily be used, too. In this case, any additional keyword arguments are directly passed to the given model-

```jldoctest
julia> using PSFModels

julia> inject(zeros(5, 5), gaussian; x=4, y=3, amp=2, fwhm=2)
5×5 Matrix{Float64}:
 0.000244141  0.00195312  0.00390625  0.00195312  0.000244141
 0.0078125    0.0625      0.125       0.0625      0.0078125
 0.0625       0.5         1.0         0.5         0.0625
 0.125        1.0         2.0         1.0         0.125
 0.0625       0.5         1.0         0.5         0.0625

julia> inject(zeros(5, 5), airydisk, 45; x=4, y=3, amp=2, fwhm=2, ratio=0.3)
5×5 Matrix{Float64}:
 0.042946   0.0944566   0.0741111  0.06456   0.0914308
 0.0944566  0.0191043   0.054543   0.108174  0.00031632
 0.0741111  0.054543    0.929442   1.29573   0.307196
 0.06456    0.108174    1.29573    1.76798   0.470669
 0.0914308  0.00031632  0.307196   0.470669  0.0620712
```

"""
inject(frame::AbstractMatrix, args...; kwargs...) = inject!(copy(frame), args...; kwargs...)

function inject!(frame::AbstractMatrix, kernel::AbstractMatrix{T}, args...; degree=Linear(), fill=zero(T), kwargs...) where T
    etp = ImageTransformations.box_extrapolation(kernel; method=degree, fillvalue=fill)
    return inject!(frame, etp, args...; kwargs...)
end

function inject!(frame::AbstractMatrix{T}, kernel::AbstractExtrapolation, angle=0; x, y, amp=one(T), center=center(frame), inds=CartesianIndices(frame), kwargs...) where T
    if iszero(angle)
        idxmap = identity
    else
        idxmap = recenter(RotMatrix{2}(deg2rad(angle)), center)
    end
    tform = Translation(ImageTransformations.center(kernel) .- (x, y)) ∘ idxmap
    @inbounds for idx in inds
        point = tform(SVector(idx.I))
        frame[idx] += amp * kernel(point...)
    end
    return frame
end

function inject!(frame::AbstractMatrix{T}, psfmodel, angle=0; center=center(frame), inds=CartesianIndices(frame), kwargs...) where T
    if iszero(angle)
        idxmap = identity
    else
        idxmap = recenter(RotMatrix{2}(deg2rad(angle)), center)
    end
    @inbounds for idx in inds
        point = idxmap(SVector(idx.I))
        frame[idx] += psfmodel(T, point; kwargs...)
    end
    return frame
end

################################

"""
    inject(cube, psf, [angles]; x, y, amp=1, kwargs...)

Injects the matrix `psf` into each frame of `cube` at the given position.

If `angles` are provided, each position will be rotated by the given parallactic angles in degrees, centered at `center`, which is the center of each frame, by default.
"""
inject(cube::AbstractArray{T,3}, args...; kwargs...) where T =
    inject!(copy(cube), args...; kwargs...)

"""
    inject!(cube, psf, [angles]; kwargs...)

In-place version of [`inject`](@ref) which modifies `cube`.
"""
function inject!(cube::AbstractArray{T,3}, kernel, angles=Zeros(size(cube, 3)); kwargs...) where T
    # All zeros position angles is actually 90° parallactic angle
    @inbounds for idx in axes(cube, 3)
        frame = @view cube[:, :, idx]
        inject!(frame, kernel, angles[idx]; kwargs...)
    end
    return cube
end

function inject!(cube::AnnulusView{T}, kernel, angles=Zeros(size(cube, 3)); kwargs...) where T
    # All zeros position angles is actually 90° parallactic angle
    @inbounds for idx in cube.indices[2]
        # relies on the fact that the underlying data can still be modified
        frame = @view cube[:, :, idx]
        inject!(frame, kernel, angles[idx]; inds=cube.indices[1], kwargs...)
    end
    return cube
end

function inject!(cube::MultiAnnulusView{T}, kernel, angles=Zeros(size(cube, 3)); kwargs...) where T
    # All zeros position angles is actually 90° parallactic angle
    @inbounds for annidx in cube.indices
        for idx in annidx[2]
            inds = annidx[1]
            # relies on the fact that the underlying data can still be modified
            frame = @view cube[:, :, idx]
            inject!(frame, kernel, angles[idx]; inds, kwargs...)
        end
    end
    return cube
end

ext: jl


text:
module EVE

const EVE_JL_VERSION = v"0.0.1"

include("tools.jl")
include("images.jl")
include("api.jl")
include("crest.jl")

export Crest

end # module


text:
using
    Oceananigans,
    Oceananigans.Operators,
    Oceananigans.TurbulenceClosures,
    Test,
    Random,
    JLD2,
    Printf,
    Statistics,
    OffsetArrays,
    FFTW

@hascuda using CuArrays

using Statistics: mean
using LinearAlgebra: norm
using GPUifyLoops: @launch, @loop

using Oceananigans: PoissonSolver, PPN, PNN, solve_poisson_3d!,
                    velocity_div!, compute_w_from_continuity!,
                    launch_config, datatuples, device,
                    parentdata, fill_halo_regions!, run_diagnostic,
                    TracerFields, buoyancy_frequency_squared, thermal_expansion, haline_contraction, ρ′,
                    RoquetIdealizedNonlinearEquationOfState

import Oceananigans: datatuple

using Oceananigans.TurbulenceClosures: ∂x_caa, ∂x_faa, ∂x²_caa, ∂x²_faa,
                                       ∂y_aca, ∂y_afa, ∂y²_aca, ∂y²_afa,
                                       ∂z_aac, ∂z_aaf, ∂z²_aac, ∂z²_aaf,
                                       ▶x_caa, ▶x_faa, ▶y_aca, ▶y_afa,
                                       ▶z_aac, ▶z_aaf

float_types = (Float32, Float64)

archs = (CPU(),)
@hascuda archs = (CPU(), GPU())

closures = (
            :ConstantIsotropicDiffusivity,
            :ConstantAnisotropicDiffusivity,
            :SmagorinskyLilly,
            :BlasiusSmagorinsky,
            :RozemaAnisotropicMinimumDissipation,
            :VerstappenAnisotropicMinimumDissipation
           )

EquationsOfState = (LinearEquationOfState, RoquetIdealizedNonlinearEquationOfState)

@testset "Oceananigans" begin
    include("test_grids.jl")
    include("test_fields.jl")
    include("test_halo_regions.jl")
    include("test_operators.jl")
    include("test_poisson_solvers.jl")
    include("test_coriolis.jl")
    include("test_buoyancy.jl")
    include("test_models.jl")
    include("test_time_stepping.jl")
    include("test_boundary_conditions.jl")
    include("test_forcings.jl")
    include("test_turbulence_closures.jl")
    include("test_dynamics.jl")
    include("test_diagnostics.jl")
    include("test_output_writers.jl")
    include("test_regression.jl")
    include("test_examples.jl")
end


text:
# Copyright (c) 2021, Joseph M. Rutherford
# All rights reserved.
#
# Code provided under the license contained in the LICENSE file.

"""
    resolve_subarray!(test)

Traverse a portion of a SolvablePuzzle to reduce logical options.
"""
function resolve_subarray!(test::Array{Sudoku.PuzzleEntry})
    rank_squared = length(test)

    # Visit each element to determine new mask
    mask = BitVector(undef,rank_squared)
    mask .= true
    for i = 1:rank_squared
        value = get_value(test[i])
        if value != 0
            mask[value] = false
        end
    end

    # Apply mask to all unknowns
    for i = 1:rank_squared
        if get_value(test[i]) == 0
            update_possibilities = BitVector(undef,rank_squared)
            update_possibilities.chunks[1] = test[i].possibilities
            update_possibilities .&= mask
            update_value = 0
            # if we have exactly 1 possibility, it's not unknown!
            possibility_sum = 0
            possible_value = 0
            for j = 1:rank_squared
                if update_possibilities[j]
                    possibility_sum += 1
                    possible_value = j
                end
            end
            test[i] = Sudoku.PuzzleEntry(update_possibilities)
        end
    end
    return nothing
end

"""
    resolve_puzzle(puzzle)

Traverse a puzzle by row, column, and block to apply logical rules.
"""
function resolve_puzzle!(puzzle::SolvablePuzzle)
    rank=get_rank(puzzle)
    rank_squared = rank*rank
    # use logical AND throughout loops
    for rowcol = 1:rank_squared
        # check full row and col
        update = puzzle.grid[rowcol,:]
        resolve_subarray!(update)
        puzzle.grid[rowcol,:] .= update
        update = puzzle.grid[:,rowcol]
        resolve_subarray!(update)
        puzzle.grid[:,rowcol] .= update
    end
    # Check rank x rank blocks
    for row_block = 1:rank
        start_row = rank*(row_block-1)+1
        stop_row = rank*(row_block-1)+rank
        for col_block = 1:rank
            start_col = rank*(col_block-1)+1
            stop_col = rank*(col_block-1)+rank
            update = puzzle.grid[start_row:stop_row, start_col:stop_col]
            resolve_subarray!(update)
            puzzle.grid[start_row:stop_row, start_col:stop_col] .= update
        end
    end
    return nothing
end

"""
    resolve_subarray!(test,n)

Traverse a portion of a SolvablePuzzle to reduce logical options for n cells claiming n values.
"""
function resolve_subarray!(test::Array{Sudoku.PuzzleEntry},n::Integer)
    rank_squared = length(test)
    for i = 1:rank_squared
        p1 = BitVector(undef,rank_squared)
        p1.chunks[1] = test[i].possibilities
        if sum(p1) != n
            continue # nothing to do
        end
        # Fill a vector of values this cell claims
        possible_values = Set{Integer}()
        for j = 1:rank_squared
            if p1[j]
                possible_values = union(possible_values,Set([j]))
            end
        end # secondary loop over possible values
        # Track what array indices claim a size-n match
        match_indices = Set([i])
        for j = 1:rank_squared
            # Compare integers holding state of possibilities
            if i != j && test[i].possibilities == test[j].possibilities
                match_indices = union(match_indices,Set([j]))
            end
        end # secondary loop to find size-n matches to p1
        count_matches_found = length(match_indices)
        # If < n entries must have those n values, do nothing
        # If exactly n entries must have those n values, they cannot be in any others
        # If more than n entries must have n values, this is impossible
        if count_matches_found < n
            continue
        elseif count_matches_found > n
            throw(DomainError("Cannot match more than n times"))
        end
        #  elseif count_matches_found == n
        # remove these entries from other cell's possibilities
        for j = 1:rank_squared
            if j in match_indices
                continue
            else
                # Force values to zero
                p2 = BitVector(undef,rank_squared)
                p2.chunks[1] = test[j].possibilities
                for v in possible_values
                    p2[v] = false
                end
                test[j] = PuzzleEntry(p2)
            end # branch 'remove possible values'
        end # secondary loop to remove entries
        # End of processing for matching possibilities for outermost chosen cell
    end # primary loop looking for n-size claims to possible value
    return nothing
end

"""
    resolve_puzzle!(puzzle,n)

Apply n-cell compound rule analyses to reduce uncertainty in puzzle.
"""
function resolve_puzzle!(puzzle::SolvablePuzzle, n::Integer)
    if n == 1
        resolve_puzzle!(puzzle)
        return nothing
    end
    # Handle shortcut for special case
    rank = get_rank(puzzle)
    rank_squared = rank*rank
    if n >= rank_squared
        throw(DomainError("Cannot evaluate compound resolve of size >= rank_squared"))
    end
    # use logical AND throughout loops
    for rowcol = 1:rank_squared
        # check full row and col
        update = puzzle.grid[rowcol,:]
        resolve_subarray!(update,n)
        puzzle.grid[rowcol,:] .= update
        update = puzzle.grid[:,rowcol]
        resolve_subarray!(update,n)
        puzzle.grid[:,rowcol] .= update
    end
    # Check rank x rank blocks
    for row_block = 1:rank
        start_row = rank*(row_block-1)+1
        stop_row = rank*(row_block-1)+rank
        for col_block = 1:rank
            start_col = rank*(col_block-1)+1
            stop_col = rank*(col_block-1)+rank
            update = puzzle.grid[start_row:stop_row, start_col:stop_col]
            resolve_subarray!(update,n)
            puzzle.grid[start_row:stop_row, start_col:stop_col] .= update
        end
    end
    # Compound resolve done for all rows, cols, and blocks
    return nothing
end

"""
    iterative_solve!(puzzle, n)

Iteratively resolve rows, columns, and blocks of a puzzle.
Compound rules evaluation of n >= 0 are used; n == 0 is a no-op.
"""
function iterative_solve!(puzzle::SolvablePuzzle, n::Integer)
    # Uncertainty is rank_squared*puzzle_size
    rank = get_rank(puzzle)
    rank_squared = rank*rank
    maximum_uncertainty = rank_squared*rank_squared*rank_squared

    iteration = 0
    previous_uncertainty = maximum_uncertainty
    current_uncertainty = uncertainty(puzzle)
    while current_uncertainty > 0
        if current_uncertainty > previous_uncertainty
            throw(DomainError("Solution diverging"))
        elseif current_uncertainty == previous_uncertainty
            break # Stalled -- stop
        else
            for i in 1:n
                resolve_puzzle!(puzzle,i)
            end
            previous_uncertainty = current_uncertainty
            current_uncertainty = uncertainty(puzzle)
        end
        iteration += 1
    end
    return [iteration,current_uncertainty]
end

"""
    guess_solutions(puzzle)

Return a vector of single-entry guesses for the unknowns.
All guesses correspond to possible values of a single entry.
Each guessed new puzzle may be attempted for further solution.
If valid after solution, the guess must satisfy puzzle.
"""
function guess_solutions(puzzle::Sudoku.SolvablePuzzle)
    # Find the point of greatest uncertainty
    peak_possibility = 0
    peak_index = 0
    rank = Sudoku.get_rank(puzzle.grid)
    rank_squared = rank*rank
    for i = 1:length(puzzle.grid)
        temp = BitVector(undef,rank_squared)
        temp.chunks[1] = puzzle.grid[i].possibilities # Convert int to BitVector
        possibility_count = sum(temp)
        if possibility_count > peak_possibility
            peak_possibility = possibility_count
            peak_index = i
        end
    end
    if peak_possibility < 1
        throw(DomainError("Cannot solve an overdetermined puzzle"))
    end
    # Now make guesses from this this one entry
    results = Array{Sudoku.SolvablePuzzle}(undef,peak_possibility)
    for i in 1:peak_possibility
        results[i] = SolvablePuzzle(rank)
        results[i].grid .= puzzle.grid
    end
    next_result_index = 1
    entry = BitVector(undef,rank_squared)
    entry.chunks[1] = puzzle.grid[peak_index].possibilities
    for i = 1:rank_squared
        if entry[i]
            # Attempt to solve with this guess
            guess_entry = BitVector(undef,rank_squared)
            guess_entry .= false
            guess_entry[i] = true
            results[next_result_index].grid[peak_index] = Sudoku.PuzzleEntry(guess_entry)
            next_result_index += 1
        end
    end
    return results
    if next_result_index != (peak_possibility+1)
        throw(ErrorException("Guessing solutions did not correctly traverse the uncertainty"))
    end
    return result
end

"""
    backtrack_solve(puzzle,n,r,c)

Attempt solution by guessing an instance of all unknown values.
Compound rule size n >= 0 are evaluated to reduce search space (n==0 means no reductions attempted).
Indeterminate solutions are recursed up to the given limit r.
Maximum solution count max at which evaluation will terminate.
"""
function backtrack_solve(puzzle::Sudoku.SolvablePuzzle,n::Integer,r::Integer,max::Integer)
    result = Array{SolvablePuzzle}(undef,0)
    if r < 1 || max < 1
        return result
    end
    # Gather a small volley of guesses from the single most unknown entry
    guesses = Sudoku.guess_solutions(puzzle)
    for guess in guesses
        if length(result) == max # Check for termination criteria
            break
        end
        try
            iterations, uncertainty = Sudoku.iterative_solve!(guess,n)
            if n > 0 && uncertainty == 0
                if Sudoku.valid_puzzle(Sudoku.as_values(guess))
                    push!(result,guess)
                end
            else # guaranteed above r > 0
                backtrack = backtrack_solve(guess,n,r-1,max-length(result))
                for b in backtrack
                    if Sudoku.valid_puzzle(Sudoku.as_values(b))
                        push!(result,b)
                    end
                    if length(result) == max # Check for termination criteria
                        break
                    end
                end
            end
        catch e
            # DomainErrors may occur with bad guesses
            if !isa(e,DomainError)
                throw(e)
            end
        end # try block
    end # loop over all guesses
    # Return whatever is found
    return result
end

max_stars_repo_name: joseph-m-rutherford/sudoku-julia

max_issues_repo_name: joseph-m-rutherford/sudoku-julia

max_forks_repo_name: joseph-m-rutherford/sudoku-julia


text:
@testset "Diagonal ambiguity methods" begin
    d = Diagonal([1,2,3])
    v = [2,3,4]
    tv = v.'

    @test (tv*d)::TransposedVector == [2,6,12].'
    @test_throws Exception d*tv

    @test (d*tv.')::Vector == [2,6,12]

    @test_throws Exception tv.'*d

    @test (d*tv')::Vector == [2,6,12]

    @test_throws Exception tv'*d

    @test (tv/d)::TransposedVector ≈ [2/1  3/2  4/3]

    @test_throws Exception d \ tv

end

@testset "Bidiagonal ambiguity methods" begin
    bd = Bidiagonal([1,2,3], [0,0], true)
    v = [2,3,4]
    tv = v.'

    @test (tv/bd)::TransposedVector ≈ [2/1  3/2  4/3]

    @test_throws Exception bd \ tv
end


text:
mutable struct TennisGame2
    P1point::Int
	P2point::Int

	P1res::String
	P2res::String
	player1Name::String
	player2Name::String

    TennisGame2(player1Name, player2Name) = new(0, 0, "", "", player1Name, player2Name)
end

function Tenniskata.get_score(game::TennisGame2)
    score = ""
	if game.P1point == game.P2point && game.P1point < 4
		if game.P1point == 0
			score = "Love"
        end
		if game.P1point == 1
			score = "Fifteen"
        end
		if game.P1point == 2
			score = "Thirty"
        end
		score *= "-All"
    end
	if game.P1point == game.P2point && game.P1point >= 3
		score = "Deuce"
    end

	if game.P1point > 0 && game.P2point == 0
		if game.P1point == 1
			game.P1res = "Fifteen"
        end
		if game.P1point == 2
			game.P1res = "Thirty"
        end
		if game.P1point == 3
			game.P1res = "Forty"
        end

		game.P2res = "Love"
		score = game.P1res * "-" * game.P2res
    end
	if game.P2point > 0 && game.P1point == 0
		if game.P2point == 1
			game.P2res = "Fifteen"
        end
		if game.P2point == 2
			game.P2res = "Thirty"
        end
		if game.P2point == 3
			game.P2res = "Forty"
        end

		game.P1res = "Love"
		score = game.P1res * "-" * game.P2res
    end

	if game.P1point > game.P2point && game.P1point < 4
		if game.P1point == 2
			game.P1res = "Thirty"
        end
		if game.P1point == 3
			game.P1res = "Forty"
        end
		if game.P2point == 1
			game.P2res = "Fifteen"
        end
		if game.P2point == 2
			game.P2res = "Thirty"
        end
		score = game.P1res * "-" * game.P2res
    end
	if game.P2point > game.P1point && game.P2point < 4
		if game.P2point == 2
			game.P2res = "Thirty"
        end
		if game.P2point == 3
			game.P2res = "Forty"
        end
		if game.P1point == 1
			game.P1res = "Fifteen"
        end
		if game.P1point == 2
			game.P1res = "Thirty"
        end
		score = game.P1res * "-" * game.P2res
    end

	if game.P1point > game.P2point && game.P2point >= 3
		score = "Advantage player1"
    end

	if game.P2point > game.P1point && game.P1point >= 3
		score = "Advantage player2"
    end

	if game.P1point >= 4 && game.P2point >= 0 && (game.P1point-game.P2point) >= 2
		score = "Win for player1"
    end
	if game.P2point >= 4 && game.P1point >= 0 && (game.P2point-game.P1point) >= 2
		score = "Win for player2"
    end
	return score
end

set_p1_score(game::TennisGame2, number::Int) = for _ = 1:number p1_score(game) end
set_p2_score(game::TennisGame2, number::Int) = for _ = 1:number p2_score(game) end

function p1_score(game::TennisGame2)
    game.P1point += 1
end

function p2_score(game::TennisGame2)
    game.P2point += 1
end

function Tenniskata.won_point(game::TennisGame2, player::String)
	if player == "player1"
		p1_score(game)
	else
		p2_score(game)
    end
end

max_stars_repo_name: jacob-alt-del/Tennis-Refactoring-Kata

max_issues_repo_name: jacob-alt-del/Tennis-Refactoring-Kata

max_forks_repo_name: jacob-alt-del/Tennis-Refactoring-Kata


text:
module VecchiaFactorization

import Base: size, getindex, permute!, invpermute!, parent, show, rand, randn

using LinearAlgebra # BLAS.set_num_threads(1)
import LinearAlgebra: mul!, lmul!, ldiv!, \, /, *, inv, pinv,
adjoint, transpose, Matrix, sqrt, Hermitian, Symmetric, cholesky

using SparseArrays: spdiagm
import SparseArrays: sparse

using ArrayLayouts # supposed to speed up mul! for Symmetric, etc...
				   # but for small matrices this doesn't appear to help

using BlockArrays: PseudoBlockArray,  Block,
blocks, blocksizes, blockedrange, findblockindex, blockindex, mortar

export Ridiagonal, Midiagonal, Inv, Piv, sparse

# VecchiaFactor{T}
# ===========================================

abstract type VecchiaFactor{T} <: AbstractMatrix{T} end

# Inv which compliments Adjoint
# ===========================================
include("lazy_inv.jl")

const Adj_VecchiaFactor{T}     = Adjoint{T,<:VecchiaFactor}
const Inv_VecchiaFactor{T}     = Inv{T,<:VecchiaFactor}
const Inv_Adj_VecchiaFactor{T} = Inv{T,<:Adj_VecchiaFactor}
const InvOrAdjOrVecc = Union{VecchiaFactor, Adj_VecchiaFactor, Inv_VecchiaFactor, Inv_Adj_VecchiaFactor}

# inv creates Inv generically (you can bypass these)
inv(A::VecchiaFactor)         = Inv(A)    # -> Inv_VecchiaFactor
inv(A::Inv_VecchiaFactor)     = A.parent  # -> VecchiaFactor
inv(A::Adj_VecchiaFactor)     = Inv(A)    # -> Inv_Adj_VecchiaFactor
inv(A::Inv_Adj_VecchiaFactor) = A.parent  # -> Adj_VecchiaFactor
pinv(A::InvOrAdjOrVecc)       = inv(A)

# adjoint automatically creats an Adjoint (you can bypass these)
# adjoint(A::VecchiaFactor)     # automatic -> Adj_VecchiaFactor
adjoint(A::Inv_VecchiaFactor)     = Inv(adjoint(A.parent))  # -> Inv_Adj_VecchiaFactor
adjoint(A::Adj_VecchiaFactor)     = A.parent               # -> VecchiaFactor
adjoint(A::Inv_Adj_VecchiaFactor) = Inv(adjoint(A.parent)) # -> Inv_VecchiaFactor

# for operating on vectors the base methods are these (define them for each new type)
# mul!(w, ::VecchiaFactor, v)
# mul!(w, ::Adj_VecchiaFactor, v)
# ldiv!(w, ::VecchiaFactor, v)
# ldiv!(w, ::Adj_VecchiaFactor, v)

# * calls out to mul! and ldiv!
*(VF::A,   w::AbstractVector) where {A<:VecchiaFactor}         = mul!(copy(w), VF, w)
*(VFᴴ::A,  w::AbstractVector) where {A<:Adj_VecchiaFactor}     = mul!(copy(w), VFᴴ, w)
*(iVF::A,  w::AbstractVector) where {A<:Inv_VecchiaFactor}     = ldiv!(iVF.parent, copy(w))
*(iVFᴴ::A, w::AbstractVector) where {A<:Inv_Adj_VecchiaFactor} = ldiv!(iVFᴴ.parent, copy(w))

# \ calls out to inv * ...
\(VF::A,   w::AbstractVector) where {A<:VecchiaFactor}         = inv(VF) * w
\(VFᴴ::A,  w::AbstractVector) where {A<:Adj_VecchiaFactor}     = inv(VFᴴ) * w
\(iVF::A,  w::AbstractVector) where {A<:Inv_VecchiaFactor}     = inv(iVF) * w
\(iVFᴴ::A, w::AbstractVector) where {A<:Inv_Adj_VecchiaFactor} = inv(iVFᴴ) * w

# Chain products by creating a tuple
# ===========================================
include("op_chain.jl")

# two specific VecchiaFactors
# ===========================================
include("mi_ri.jl")
include("pivot_type.jl") # TODO: standardize this as a Vecchia factor

# constructing sparse or matrix equivalents
# ===========================================
include("sparse_matrix_show.jl")

# constructor Vecchia factorization
# ===========================================
include("vecchia_approx.jl")

# ===========================================
include("util.jl")

end


text:
module SIMD

# The Julia SIMD vector type

export Vec
immutable Vec{N,T} <: DenseArray{T,1}
    elts::NTuple{N,T}
    Vec(elts::NTuple{N,T}) = new(elts)
end

# Type properties

# Note: eltype and ndims are provided by DenseArray

import Base: length
length{N,T}(::Type{Vec{N,T}}) = N
length{N,T}(::Vec{N,T}) = N

import Base: size
size{N,T}(::Vec{N,T}) = (N,)
size{N,T}(::Type{Vec{N,T}}) = (N,)
size{N,T}(::Vec{N,T}, n::Int) = size(Vec{N,T})[n]
size{N,T}(::Type{Vec{N,T}}, n::Int) = size(Vec{N,T})[n]

# Type conversion

import Base: convert
@generated function create{N,T}(::Type{Vec{N,T}}, x::T)
    quote
        $(Expr(:meta, :inline))
        Vec{N,T}($(Expr(:tuple, [:x for i in 1:N]...)))
    end
end
convert{N,T}(::Type{Vec{N,T}}, x::T) = create(Vec{N,T}, x)
convert{N,T}(::Type{Vec{N,T}}, x::Number) = create(Vec{N,T}, T(x))
convert{N,T}(::Type{Vec{N,T}}, xs::NTuple{N}) = Vec{N,T}(NTuple{N,T}(xs))

convert{N,T}(::Type{NTuple{N,T}}, v::Vec{N,T}) = v.elts

# Convert Julia types to LLVM types

llvmtype(::Type{Bool}) = "i8"   # Julia represents Tuple{Bool} as [1 x i8]

llvmtype(::Type{Int8}) = "i8"
llvmtype(::Type{Int16}) = "i16"
llvmtype(::Type{Int32}) = "i32"
llvmtype(::Type{Int64}) = "i64"
llvmtype(::Type{Int128}) = "i128"

llvmtype(::Type{UInt8}) = "i8"
llvmtype(::Type{UInt16}) = "i16"
llvmtype(::Type{UInt32}) = "i32"
llvmtype(::Type{UInt64}) = "i64"
llvmtype(::Type{UInt128}) = "i128"

llvmtype(::Type{Float16}) = "half"
llvmtype(::Type{Float32}) = "float"
llvmtype(::Type{Float64}) = "double"

# Type-dependent optimization flags
fastflags{T<:Signed}(::Type{T}) = "nsw"
fastflags{T<:Unsigned}(::Type{T}) = "nuw"
fastflags{T<:AbstractFloat}(::Type{T}) = "fast"

suffix{T}(N::Integer, ::Type{T}) = "v$(N)f$(8*sizeof(T))"

# Type-dependent LLVM constants
function llvmconst{T}(N::Integer, ::Type{T}, val)
    T(val) === T(0) && return "zeroinitializer"
    typ = llvmtype(T)
    "<" * join(["$typ $val" for i in 1:N], ", ") * ">"
end
function llvmconst(N::Integer, ::Type{Bool}, val)
    Bool(val) === false && return "zeroinitializer"
    typ = "i1"
    "<" * join(["$typ $(Int(val))" for i in 1:N], ", ") * ">"
end

# Type-dependent LLVM intrinsics
llvmins{T<:Integer}(::Type{Val{:+}}, N, ::Type{T}) = "add"
llvmins{T<:Integer}(::Type{Val{:-}}, N, ::Type{T}) = "sub"
llvmins{T<:Integer}(::Type{Val{:*}}, N, ::Type{T}) = "mul"
llvmins{T<:Signed}(::Type{Val{:div}}, N, ::Type{T}) = "sdiv"
llvmins{T<:Signed}(::Type{Val{:rem}}, N, ::Type{T}) = "srem"
llvmins{T<:Unsigned}(::Type{Val{:div}}, N, ::Type{T}) = "udiv"
llvmins{T<:Unsigned}(::Type{Val{:rem}}, N, ::Type{T}) = "urem"

llvmins{T<:Integer}(::Type{Val{:~}}, N, ::Type{T}) = "xor"
llvmins{T<:Integer}(::Type{Val{:&}}, N, ::Type{T}) = "and"
llvmins{T<:Integer}(::Type{Val{:|}}, N, ::Type{T}) = "or"
llvmins{T<:Integer}(::Type{Val{:$}}, N, ::Type{T}) = "xor"

llvmins{T<:Integer}(::Type{Val{:<<}}, N, ::Type{T}) = "shl"
llvmins{T<:Integer}(::Type{Val{:>>>}}, N, ::Type{T}) = "lshr"
llvmins{T<:Unsigned}(::Type{Val{:>>}}, N, ::Type{T}) = "lshr"
llvmins{T<:Signed}(::Type{Val{:>>}}, N, ::Type{T}) = "ashr"

llvmins{T<:Integer}(::Type{Val{:(==)}}, N, ::Type{T}) = "icmp eq"
llvmins{T<:Integer}(::Type{Val{:(!=)}}, N, ::Type{T}) = "icmp ne"
llvmins{T<:Signed}(::Type{Val{:(>)}}, N, ::Type{T}) = "icmp sgt"
llvmins{T<:Signed}(::Type{Val{:(>=)}}, N, ::Type{T}) = "icmp sge"
llvmins{T<:Signed}(::Type{Val{:(<)}}, N, ::Type{T}) = "icmp slt"
llvmins{T<:Signed}(::Type{Val{:(<=)}}, N, ::Type{T}) = "icmp sle"
llvmins{T<:Unsigned}(::Type{Val{:(>)}}, N, ::Type{T}) = "icmp ugt"
llvmins{T<:Unsigned}(::Type{Val{:(>=)}}, N, ::Type{T}) = "icmp uge"
llvmins{T<:Unsigned}(::Type{Val{:(<)}}, N, ::Type{T}) = "icmp ult"
llvmins{T<:Unsigned}(::Type{Val{:(<=)}}, N, ::Type{T}) = "icmp ule"

llvmins{T}(::Type{Val{:ifelse}}, N, ::Type{T}) = "select"

llvmins{T<:AbstractFloat}(::Type{Val{:+}}, N, ::Type{T}) = "fadd"
llvmins{T<:AbstractFloat}(::Type{Val{:-}}, N, ::Type{T}) = "fsub"
llvmins{T<:AbstractFloat}(::Type{Val{:*}}, N, ::Type{T}) = "fmul"
llvmins{T<:AbstractFloat}(::Type{Val{:/}}, N, ::Type{T}) = "fdiv"
llvmins{T<:AbstractFloat}(::Type{Val{:inv}}, N, ::Type{T}) = "fdiv"
llvmins{T<:AbstractFloat}(::Type{Val{:rem}}, N, ::Type{T}) = "frem"

llvmins{T<:AbstractFloat}(::Type{Val{:(==)}}, N, ::Type{T}) = "fcmp oeq"
llvmins{T<:AbstractFloat}(::Type{Val{:(!=)}}, N, ::Type{T}) = "fcmp une"
llvmins{T<:AbstractFloat}(::Type{Val{:(>)}}, N, ::Type{T}) = "fcmp ogt"
llvmins{T<:AbstractFloat}(::Type{Val{:(>=)}}, N, ::Type{T}) = "fcmp oge"
llvmins{T<:AbstractFloat}(::Type{Val{:(<)}}, N, ::Type{T}) = "fcmp olt"
llvmins{T<:AbstractFloat}(::Type{Val{:(<=)}}, N, ::Type{T}) = "fcmp ole"

llvmins{T<:AbstractFloat}(::Type{Val{:^}}, N, ::Type{T}) =
    "@llvm.pow.$(suffix(N,T))"
llvmins{T<:AbstractFloat}(::Type{Val{:abs}}, N, ::Type{T}) =
    "@llvm.fabs.$(suffix(N,T))"
llvmins{T<:AbstractFloat}(::Type{Val{:ceil}}, N, ::Type{T}) =
    "@llvm.ceil.$(suffix(N,T))"
llvmins{T<:AbstractFloat}(::Type{Val{:copysign}}, N, ::Type{T}) =
    "@llvm.copysign.$(suffix(N,T))"
llvmins{T<:AbstractFloat}(::Type{Val{:cos}}, N, ::Type{T}) =
    "@llvm.cos.$(suffix(N,T))"
llvmins{T<:AbstractFloat}(::Type{Val{:exp}}, N, ::Type{T}) =
    "@llvm.exp.$(suffix(N,T))"
llvmins{T<:AbstractFloat}(::Type{Val{:exp2}}, N, ::Type{T}) =
    "@llvm.exp2.$(suffix(N,T))"
llvmins{T<:AbstractFloat}(::Type{Val{:floor}}, N, ::Type{T}) =
    "@llvm.floor.$(suffix(N,T))"
llvmins{T<:AbstractFloat}(::Type{Val{:fma}}, N, ::Type{T}) =
    "@llvm.fma.$(suffix(N,T))"
llvmins{T<:AbstractFloat}(::Type{Val{:log}}, N, ::Type{T}) =
    "@llvm.log.$(suffix(N,T))"
llvmins{T<:AbstractFloat}(::Type{Val{:log10}}, N, ::Type{T}) =
    "@llvm.log10.$(suffix(N,T))"
llvmins{T<:AbstractFloat}(::Type{Val{:log2}}, N, ::Type{T}) =
    "@llvm.log2.$(suffix(N,T))"
llvmins{T<:AbstractFloat}(::Type{Val{:max}}, N, ::Type{T}) =
    "@llvm.maxnum.$(suffix(N,T))"
llvmins{T<:AbstractFloat}(::Type{Val{:min}}, N, ::Type{T}) =
    "@llvm.minnum.$(suffix(N,T))"
llvmins{T<:AbstractFloat}(::Type{Val{:muladd}}, N, ::Type{T}) =
    "@llvm.fmuladd.$(suffix(N,T))"
llvmins{T<:AbstractFloat}(::Type{Val{:powi}}, N, ::Type{T}) =
    "@llvm.powi.$(suffix(N,T))"
llvmins{T<:AbstractFloat}(::Type{Val{:round}}, N, ::Type{T}) =
    "@llvm.rint.$(suffix(N,T))"
llvmins{T<:AbstractFloat}(::Type{Val{:sin}}, N, ::Type{T}) =
    "@llvm.sin.$(suffix(N,T))"
llvmins{T<:AbstractFloat}(::Type{Val{:sqrt}}, N, ::Type{T}) =
    "@llvm.sqrt.$(suffix(N,T))"
llvmins{T<:AbstractFloat}(::Type{Val{:trunc}}, N, ::Type{T}) =
    "@llvm.trunc.$(suffix(N,T))"

# Convert between LLVM scalars, vectors, and arrays

function scalar2vector(vec, siz, typ, sca)
    instrs = []
    accum(nam, i) = i<0 ? "undef" : i==siz-1 ? nam : "$(nam)_$i"
    for i in 0:siz-1
        push!(instrs,
            "$(accum(vec,i)) = " *
                "insertelement <$siz x $typ> $(accum(vec,i-1)), " *
                "$typ $sca, i32 $i")
    end
    instrs
end

function scalar2array(varrec, siz, typ, sca)
    instrs = []
    accum(nam, i) = i<0 ? "undef" : i==siz-1 ? nam : "$(nam)_$i"
    for i in 0:siz-1
        push!(instrs,
            "$(accum(arr,i)) = " *
                "insertvalue [$siz x $typ] $(accum(arr,i-1)), $typ $sca, $i")
    end
    instrs
end

function array2vector(vec, siz, typ, arr, tmp=arr)
    instrs = []
    accum(nam, i) = i<0 ? "undef" : i==siz-1 ? nam : "$(nam)_$i"
    for i in 0:siz-1
        push!(instrs, "$(tmp)_$i = extractvalue [$siz x $typ] $arr, $i")
        push!(instrs,
            "$(accum(vec,i)) = " *
                "insertelement <$siz x $typ> $(accum(vec,i-1)), " *
                "$typ $(tmp)_$i, i32 $i")
    end
    instrs
end

function vector2array(arr, siz, typ, vec, tmp=vec)
    instrs = []
    accum(nam, i) = i<0 ? "undef" : i==siz-1 ? nam : "$(nam)_$i"
    for i in 0:siz-1
        push!(instrs, "$(tmp)_$i = extractelement <$siz x $typ> $vec, i32 $i")
        push!(instrs,
            "$(accum(arr,i)) = "*
                "insertvalue [$siz x $typ] $(accum(arr,i-1)), " *
                "$typ $(tmp)_$i, $i")
    end
    instrs
end

# Element-wise access

export setindex

@generated function setindex{N,T,I}(v::Vec{N,T}, ::Type{Val{I}}, x)
    @assert isa(I, Integer)
    1 <= I <= N || throw(BoundsError())
    quote
        $(Expr(:meta, :inline))
        Vec{N,T}($(Expr(:tuple,
            [i == I ? :(T(x)) : :(v.elts[$i]) for i in 1:N]...)))
    end
end

# Note: Julia has no equivalent of LLVM's "insertelement" function
@generated function setindex{N,T}(v::Vec{N,T}, i::Integer, x::Number)
    typ = llvmtype(T)
    ityp = llvmtype(Int)
    atyp = "[$N x $typ]"
    vtyp = "<$N x $typ>"
    decls = []
    instrs = []
    append!(instrs, array2vector("%arg1", N, typ, "%0", "%arg1arr"))
    push!(instrs, "%res = insertelement $vtyp %arg1, $typ %2, $ityp %1")
    append!(instrs, vector2array("%resarr", N, typ, "%res"))
    push!(instrs, "ret $atyp %resarr")
    quote
        $(Expr(:meta, :inline))
        let j = Int(i)
            @boundscheck 1 <= j <= N || throw(BoundsError())
            Vec{N,T}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
                NTuple{N,T}, Tuple{NTuple{N,T}, Int, T}, v.elts, j-1, T(x)))
        end
    end
end

import Base: getindex
getindex{N,T,I}(v::Vec{N,T}, ::Type{Val{I}}) = v.elts[I]
getindex{N,T}(v::Vec{N,T}, i::Integer) = v.elts[i]

# Type conversion

# Note: This can also change the number of vector elements if the size of the
# vector elements changes correspondingly; it is not possible to do this in
# plain Julia
@generated function Base.reinterpret{N,R,N1,T1}(::Type{Vec{N,R}},
        v1::Vec{N1,T1})
    @assert N*sizeof(R) == N1*sizeof(T1)
    typ1 = llvmtype(T1)
    atyp1 = "[$N1 x $typ1]"
    vtyp1 = "<$N1 x $typ1>"
    typr = llvmtype(R)
    atypr = "[$N x $typr]"
    vtypr = "<$N x $typr>"
    decls = []
    instrs = []
    append!(instrs, array2vector("%arg1", N1, typ1, "%0", "%arg1arr"))
    push!(instrs, "%res = bitcast $vtyp1 %arg1 to $vtypr")
    append!(instrs, vector2array("%resarr", N, typr, "%res"))
    push!(instrs, "ret $atypr %resarr")
    quote
        $(Expr(:meta, :inline))
        Vec{N,R}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            NTuple{N,R}, Tuple{NTuple{N1,T1}}, v1.elts))
    end
end

# Conditionals

for op in (:(==), :(!=), :(<), :(<=), :(>), :(>=))
    @eval begin
        import Base: $op
        @generated function $op{N,T}(v1::Vec{N,T}, v2::Vec{N,T})
            op = $op
            quote
                $(Expr(:meta, :inline))
                Vec{N,Bool}($(Expr(:tuple,
                    [:($op(v1.elts[$i], v2.elts[$i])) for i in 1:N]...)))
            end
        end
    end
end

import Base: ifelse
@generated function ifelse{N,T}(v1::Vec{N,Bool}, v2::Vec{N,T}, v3::Vec{N,T})
    quote
        $(Expr(:meta, :inline))
        Vec{N,T}($(Expr(:tuple,
            [:(ifelse(v1.elts[$i], v2.elts[$i], v3.elts[$i]))
                for i in 1:N]...)))
    end
end

# Boolean functions

for op in (:!,)
    @eval begin
        import Base: $op
        @generated function $op{N}(v1::Vec{N,Bool})
            op = $op
            quote
                $(Expr(:meta, :inline))
                Vec{N,Bool}($(Expr(:tuple,
                    [:($op(v1.elts[$i])) for i in 1:N]...)))
            end
        end
    end
end

# Arithmetic functions

for op in (
        :~, :+, :-,
        :abs, :ceil, :cos, :exp, :exp2, :floor, :inv, :log, :log10, :log2,
        :round, :sin, :sqrt, :trunc)
    @eval begin
        import Base: $op
        @generated function $op{N,T}(v1::Vec{N,T})
            op = $op
            quote
                $(Expr(:meta, :inline))
                Vec{N,T}($(Expr(:tuple, [:($op(v1.elts[$i])) for i in 1:N]...)))
            end
        end
    end
end
import Base: exp10
@inline exp10{N,T}(v1::Vec{N,T}) = Vec{N,T}(10)^v1

for op in (
        :&, :|, :$, :+, :-, :*, :/, :^,
        :cld, :copysign, :div, :fld, :max, :min, :mod, :rem)
    @eval begin
        import Base: $op
        @generated function $op{N,T}(v1::Vec{N,T}, v2::Vec{N,T})
            op = $op
            quote
                $(Expr(:meta, :inline))
                Vec{N,T}($(Expr(:tuple,
                    [:($op(v1.elts[$i], v2.elts[$i])) for i in 1:N]...)))
            end
        end
    end
end

for op in (:fma, :muladd)
    @eval begin
        import Base: $op
        @generated function $op{N,T}(v1::Vec{N,T}, v2::Vec{N,T}, v3::Vec{N,T})
            op = $op
            quote
                $(Expr(:meta, :inline))
                Vec{N,T}($(Expr(:tuple,
                    [:($op(v1.elts[$i], v2.elts[$i], v3.elts[$i]))
                        for i in 1:N]...)))
            end
        end
    end
end

# TODO: add missing functions
# TODO: add load/store functions

end


text:
using Ipopt
using MathProgBase
using SpecialFunctions

mutable struct GenTDist <: ContinuousUnivariateDistribution
    mu::Float64
    sigma::Float64
    t::TDist
end

GenTDist( mu::Float64, sigma::Float64, nu::Float64 ) = GenTDist( mu, sigma, TDist( nu ) )

Base.rand( t::GenTDist ) = t.sigma .* rand( t.t ) .+ t.mu

Base.rand( t::GenTDist, n::Int ) = t.sigma .* rand( t.t, n ) .+ t.mu

Distributions.pdf( t::GenTDist, x::Float64 ) = pdf( t.t, (x - t.mu)/t.sigma)/t.sigma

Distributions.cdf( t::GenTDist, x::Float64 ) = cdf( t.t, (x - t.mu)/t.sigma)

Distributions.mean( t::GenTDist ) = t.mu

Distributions.std( t::GenTDist ) = t.sigma * sqrt(t.t.ν/(t.t.ν-2))

randomparameters( ::Type{GenTDist} ) = [randn(), rand( Exponential() ), 1 + rand(Exponential())]

randomparameters( ::Type{GenTDist}, n::Int ) =
    [randn( 1, n ); rand( Exponential(), 1, n ); 1 .+ rand( Exponential(), 1, n )]

mutable struct GenTDistOptimizer <: MathProgBase.AbstractNLPEvaluator
    y::Vector{Float64}
    w::Vector{Float64}
    debug::Int
end

GenTDistOptimizer( y::Vector, w::Vector ) = GenTDistOptimizer( y, w, 0 )

function MathProgBase.initialize( ::GenTDistOptimizer, requested_features::Vector{Symbol} )
    unimplemented = setdiff( requested_features, [:Grad] )
    if !isempty( unimplemented )
        error( "The following features aren't implemented: " * join( string.(unimplemented), ", " ) )
    end
end

MathProgBase.features_available( ::GenTDistOptimizer ) = [:Grad]

function MathProgBase.eval_f( t::GenTDistOptimizer, x; debug=0 )
    if t.debug > 0
        println( "eval_f called with $x" )
    end
    (mu, sigma, nu) = x

    sigma <= 0.0 && return -Inf

    y = t.y
    w = t.w
    n = sum(w)
    normalysq = ((y .- mu)/sigma).^2
    constant = n*(lgamma((nu+1)/2) - log(nu*pi)/2 - lgamma(nu/2) - log(sigma))
    return constant - (nu+1)/2*sum(w .* log.(1 .+ normalysq/nu))
end

function MathProgBase.eval_grad_f( t::GenTDistOptimizer, g, x )
    (mu, sigma, nu) = x
    y = t.y
    w = t.w
    n = sum(w)
    normaly = (y .- mu)/sigma
    normalysq = normaly.^2
    g[1] = (nu+1)*sum(w .* (y .- mu)./(nu*sigma^2 .+ (y .- mu).^2))
    g[2] = -n/sigma + (nu+1)*sum(w .* normalysq ./ (sigma*(nu .+ normalysq)))
    g[3] = n*(digamma((nu+1)/2)/2 - 1/(2*nu) - digamma(nu/2)/2)
    g[3] += -sum(w .* log.(1 .+ normalysq/nu))/2 + (nu+1)/(2*nu)*sum(w .* normalysq ./ (nu .+ normalysq))
    if t.debug > 0
        println( "grad_f called with $x, returning $g" )
    end
end

function fit_mle!(
    ::Type{GenTDist},
    parameters::AbstractVector{Out},
    x::Vector{Out},
    w::Vector{Calc},
    scratch::Dict{Symbol,Any};
    max_iter::Int = 3000,
    print_level::Int = 0,
) where {Calc,Out}
    t = GenTDistOptimizer( x, w )

    solver = IpoptSolver(print_level=print_level, max_iter=max_iter)

    model = MathProgBase.NonlinearModel(solver)
    MathProgBase.loadproblem!(model, 3, 0, [-Inf; inequalityconstraintvector(GenTDist)], fill(Inf,3), Float64[], Float64[], :Max, t)
    MathProgBase.setwarmstart!( model, parameters )
    MathProgBase.optimize!(model)

    status = MathProgBase.status(model)
    if status != :Optimal
        @warn( "Non-optimal optimization result: $status" )
    end
    parameters[:] = MathProgBase.getsolution(model)
end


text:
using BinDeps

@BinDeps.setup

const OPENJPEG_VERSION = "v2.3.0"
privides(Sources, URI("https://github.com/uclouvain/openjpeg/releases/download/v2.3.0/openjpeg-$(OPENJPEG_VERSION)-linux-x86_64.tar.gz"))


text:
function readModelAndGenerateMeshMref(m,pad::Int64,newSize::Vector,domain::Vector)
########################## m,mref are in Velocity here. ###################################

mref = copy(m);

mref = getSimilarLinearModel(m,1.75,2.9);

sea = abs.(m[:] .- minimum(m)) .< 1e-1;
mref[sea] = m[sea];
if newSize!=[]
	m    = expandModelNearest(m,   collect(size(m)),newSize);
	mref = expandModelNearest(mref,collect(size(mref)),newSize);
end

Minv = getRegularMesh(domain,collect(size(m))-1);

(mPadded,MinvPadded) = addAbsorbingLayer(m,Minv,pad);
(mrefPadded,MinvPadded) = addAbsorbingLayer(mref,Minv,pad);

N = prod(MinvPadded.n+1);
boundsLow  = minimum(mPadded);
boundsHigh = maximum(mPadded);

boundsLow  = ones(N)*boundsLow;
boundsLow = convert(Array{Float32},boundsLow);
boundsHigh = ones(N)*boundsHigh;
boundsHigh = convert(Array{Float32},boundsHigh);

return (mPadded,MinvPadded,mrefPadded,boundsHigh,boundsLow);
end


text:
struct UserQFunction{F,K}
    f::F
    fptr::Ptr{Nothing}
    kf::K
    cuf::Union{Nothing,Ptr{Nothing}}
end

@inline function extract_context(ptr, ::Type{T}) where {T}
    unsafe_load(Ptr{T}(ptr))
end

@inline function extract_array(ptr, idx, dims)
    UnsafeArray(Ptr{CeedScalar}(unsafe_load(ptr, idx)), dims)
end

function generate_user_qfunction(
    ceed,
    def_module,
    qf_name,
    constants,
    array_names,
    ctx,
    dims_in,
    dims_out,
    body,
)
    idx = gensym(:i)
    Q = gensym(:Q)
    ctx_ptr = gensym(:ctx_ptr)
    in_ptr = gensym(:in_ptr)
    out_ptr = gensym(:out_ptr)

    const_assignments = Vector{Expr}(undef, length(constants))
    for (i, c) ∈ enumerate(constants)
        const_assignments[i] = :($(c[1]) = $(c[2]))
    end

    narrays = length(array_names)
    arrays = Vector{Expr}(undef, narrays)
    array_views = Vector{Expr}(undef, narrays)
    n_in = length(dims_in)
    for (i, arr_name) ∈ enumerate(array_names)
        i_inout = (i <= n_in) ? i : i - n_in
        dims = (i <= n_in) ? dims_in[i] : dims_out[i-n_in]
        ptr = (i <= n_in) ? in_ptr : out_ptr
        arr_name_gen = gensym(arr_name)
        arrays[i] = :($arr_name_gen = extract_array($ptr, $i_inout, (Int($Q), $(dims...))))
        ndims = length(dims)
        slice = Expr(:ref, arr_name_gen, idx, (:(:) for i = 1:ndims)...)
        if i <= n_in
            if ndims == 0
                array_views[i] = :($arr_name = $slice)
            else
                S = Tuple{dims...}
                array_views[i] = :($arr_name = LibCEED.SArray{$S}(@view $slice))
            end
        else
            array_views[i] = :($arr_name = @view $slice)
        end
    end

    if isnothing(ctx)
        ctx_assignment = nothing
    else
        ctx_assignment = :($(ctx.name) = extract_context($ctx_ptr, $(ctx.type)))
    end

    qf1 = gensym(qf_name)
    f = Core.eval(
        def_module,
        quote
            @inline function $qf1(
                $ctx_ptr::Ptr{Cvoid},
                $Q::CeedInt,
                $in_ptr::Ptr{Ptr{CeedScalar}},
                $out_ptr::Ptr{Ptr{CeedScalar}},
            )
                $(const_assignments...)
                $ctx_assignment
                $(arrays...)
                @inbounds @simd for $idx = 1:$Q
                    $(array_views...)
                    $body
                end
                CeedInt(0)
            end
        end,
    )
    f_qn = QuoteNode(f)
    rt = :CeedInt
    at = :(Core.svec(Ptr{Cvoid}, CeedInt, Ptr{Ptr{CeedScalar}}, Ptr{Ptr{CeedScalar}}))
    fptr = eval(Expr(:cfunction, Ptr{Cvoid}, f_qn, rt, at, QuoteNode(:ccall)))

    # COV_EXCL_START
    if iscuda(ceed)
        getresource(ceed) == "/gpu/cuda/gen" && error(
            string(
                "/gpu/cuda/gen is not compatible with user Q-functions defined with ",
                "libCEED.jl.\nPlease use a different backend, for example: /gpu/cuda/shared ",
                "or /gpu/cuda/ref",
            ),
        )
        if isdefined(@__MODULE__, :CUDA)
            !has_cuda() && error("No valid CUDA installation found")
            qf2 = gensym(qf_name)
            kf = Core.eval(
                def_module,
                quote
                    @inline function $qf2($ctx_ptr::Ptr{Cvoid}, $(array_names...))
                        $(const_assignments...)
                        $ctx_assignment
                        $body
                        nothing
                    end
                end,
            )
            cuf = mk_cufunction(ceed, def_module, qf_name, kf, dims_in, dims_out)
        else
            error(
                string(
                    "User Q-functions with CUDA backends require the CUDA.jl package to be ",
                    "loaded.\nThe libCEED backend is: $(getresource(ceed))\n",
                    "Please ensure that the CUDA.jl package is installed and loaded.",
                ),
            )
        end
    else
        kf = nothing
        cuf = nothing
    end
    # COV_EXCL_STOP

    UserQFunction(f, fptr, kf, cuf)
end

function meta_user_qfunction(ceed, def_module, qf, args)
    qf_name = Meta.quot(qf)

    ctx = nothing
    constants = Expr[]
    dims_in = Expr[]
    dims_out = Expr[]
    names_in = Symbol[]
    names_out = Symbol[]

    for a ∈ args[1:end-1]
        if Meta.isexpr(a, :(=))
            a1 = Meta.quot(a.args[1])
            a2 = esc(a.args[2])
            push!(constants, :(($a1, $a2)))
        elseif Meta.isexpr(a, :tuple)
            arr_name = a.args[1]
            inout = a.args[2].value
            ndim = length(a.args) - 3
            dims = Vector{Expr}(undef, ndim)
            for d = 1:ndim
                dims[d] = :(Int($(a.args[d+3])))
            end
            dims_expr = :(Int[$(esc.(a.args[4:end])...)])
            if inout == :in
                push!(dims_in, dims_expr)
                push!(names_in, arr_name)
            elseif inout == :out
                push!(dims_out, dims_expr)
                push!(names_out, arr_name)
            else
                error("Array specification must be either :in or :out. Given $inout.")
            end
        elseif Meta.isexpr(a, :(::))
            ctx = (name=a.args[1], type=a.args[2])
        else
            error("Bad argument to @user_qfunction")
        end
    end

    body = Meta.quot(args[end])

    return :(generate_user_qfunction(
        $ceed,
        $def_module,
        $qf_name,
        [$(constants...)],
        $([names_in; names_out]),
        $ctx,
        [$(dims_in...)],
        [$(dims_out...)],
        $body,
    ))
end

"""
    @interior_qf name=def

Creates a user-defined interior (volumetric) Q-function, and assigns it to a variable named
`name`. The definition of the Q-function is given as:
```
@interior_qf user_qf=(
    ceed::CEED,
    [const1=val1, const2=val2, ...],
    [ctx::ContextType],
    (I1, :in, EvalMode, dims...),
    (I2, :in, EvalMode, dims...),
    (O1, :out, EvalMode, dims...),
    body
)
```
The definitions of form `const=val` are used for definitions which will be compile-time
constants in the Q-function. For example, if `dim` is a variable set to the dimension of the
problem, then `dim=dim` will make `dim` available in the body of the Q-function as a
compile-time constant.

If the user wants to provide a context struct to the Q-function, that can be achieved by
optionally including `ctx::ContextType`, where `ContextType` is the type of the context
struct, and `ctx` is the name to which is will be bound in the body of the Q-function.

This is followed by the definition of the input and output arrays, which take the form
`(arr_name, (:in|:out), EvalMode, dims...)`. Each array will be bound to a variable named
`arr_name`. Input arrays should be tagged with :in, and output arrays with :out. An
`EvalMode` should be specified, followed by the dimensions of the array. If the array
consists of scalars (one number per Q-point) then `dims` should be omitted.

# Examples

- Q-function to compute the "Q-data" for the mass operator, which is given by the quadrature
  weight times the Jacobian determinant. The mesh Jacobian (the gradient of the nodal mesh
  points) and the quadrature weights are given as input arrays, and the Q-data is the output
  array. `dim` is given as a compile-time constant, and so the array `J` is statically
  sized, and therefore `det(J)` will automatically dispatch to an optimized implementation
  for the given dimension.
```
@interior_qf build_qfunc = (
    ceed, dim=dim,
    (J, :in, EVAL_GRAD, dim, dim),
    (w, :in, EVAL_WEIGHT),
    (qdata, :out, EVAL_NONE),
    qdata[] = w*det(J)
)
```
"""
macro interior_qf(args)
    if !Meta.isexpr(args, :(=))
        error("@interior_qf must be of form `qf = (body)`") # COV_EXCL_LINE
    end

    qf = args.args[1]
    user_qf = esc(qf)
    args = args.args[2].args
    ceed = esc(args[1])

    # Calculate field sizes
    fields_in = Expr[]
    fields_out = Expr[]
    for a ∈ args
        if Meta.isexpr(a, :tuple)
            field_name = String(a.args[1])
            inout = a.args[2].value
            evalmode = a.args[3]
            ndim = length(a.args) - 3
            dims = Vector{Expr}(undef, ndim)
            for d = 1:ndim
                dims[d] = esc(:(Int($(a.args[d+3]))))
            end
            sz_expr = :(prod(($(dims...),)))
            if inout == :in
                push!(fields_in, :(add_input!($user_qf, $field_name, $sz_expr, $evalmode)))
            elseif inout == :out
                push!(
                    fields_out,
                    :(add_output!($user_qf, $field_name, $sz_expr, $evalmode)),
                )
            end
        end
    end

    gen_user_qf = meta_user_qfunction(ceed, __module__, qf, args[2:end])

    quote
        $user_qf = create_interior_qfunction($ceed, $gen_user_qf)
        $(fields_in...)
        $(fields_out...)
    end
end

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
# Syntax for mat"" string interpolation

# A really basic parser intended only to handle checking whether
# a variable is on the left or right hand side of an expression
type DumbParserState
    paren_depth::Int
    in_string::Bool
end
DumbParserState() = DumbParserState(0, false)

# Returns true if an = is encountered and updates pstate
function dumb_parse!(pstate::DumbParserState, str::AbstractString)
    paren_depth = pstate.paren_depth
    in_string = pstate.in_string
    x = '\0'
    s = start(str)
    while !done(str, s)
        lastx = x
        (x, s) = next(str, s)
        if in_string
            if x == '\''
                if !done(str, s) && next(str, s)[1] == '\''
                    (x, s) = next(str, s)
                else
                    in_string = false
                end
            end
        else
            if x == '('
                paren_depth += 1
            elseif x == ')'
                paren_depth -= 1
            elseif x == '\'' && lastx in ",( \t\0;"
                in_string = true
            elseif x == '=' && !(lastx in "<>~")
                if !done(str, s) && next(str, s)[1] == '='
                    (x, s) = next(str, s)
                else
                    return true
                end
            elseif x == '%'
                break
            end
        end
    end
    pstate.paren_depth = paren_depth
    pstate.in_string = in_string
    return false
end

# Check if a given variable is assigned, used, or both. Returns the#
# assignment and use status
function check_assignment(interp, i)
    # Go back to the last newline
    before = AbstractString[]
    for j = i-1:-1:1
        if isa(interp[j], AbstractString)
            sp = split(interp[j], "\n")
            unshift!(before, sp[end])
            for k = length(sp)-1:-1:1
                match(r"\.\.\.[ \t]*\r?$", sp[k]) == nothing && @goto done_before
                unshift!(before, sp[k])
            end
        end
    end
    @label done_before

    # Check if this reference is inside parens at the start, or on the rhs of an assignment
    pstate = DumbParserState()
    (dumb_parse!(pstate, join(before)) || pstate.paren_depth > 1) && return (false, true)

    # Go until the next newline or comment
    after = AbstractString[]
    both_sides = false
    for j = i+1:length(interp)
        if isa(interp[j], AbstractString)
            sp = split(interp[j], "\n")
            push!(after, sp[1])
            for k = 2:length(sp)
                match(r"\.\.\.[ \t]*\r?$", sp[k-1]) == nothing && @goto done_after
                push!(after, sp[k])
            end
        elseif interp[j] == interp[i]
            both_sides = true
        end
    end
    @label done_after

    assigned = dumb_parse!(pstate, join(after))
    used = !assigned || both_sides || (i < length(interp) && match(r"^[ \t]*\(", interp[i+1]) != nothing)
    return (assigned, used)
end

function do_mat_str(ex)
    # Hack to do interpolation
    interp = parse(string("\"\"\"", replace(ex, "\"\"\"", "\\\"\"\""), "\"\"\""))
    if isa(interp, AbstractString)
        interp = [interp]
    elseif interp.head == :string
        interp = interp.args
    elseif interp.head == :macrocall
        interp = interp.args[2:end]
    else
        throw(ArgumentError("unexpected input"))
    end

    # Handle interpolated variables
    putblock = Expr(:block)
    getblock = Expr(:block)
    usedvars = Set{Symbol}()
    assignedvars = Set{Symbol}()
    varmap = Dict{Symbol,Symbol}()
    for i = 1:length(interp)
        if !isa(interp[i], AbstractString)
            # Don't put the same symbol to MATLAB twice
            if haskey(varmap, interp[i])
                var = varmap[interp[i]]
            else
                var = Symbol("matlab_jl_", i)
                if isa(interp[i], Symbol)
                    varmap[interp[i]] = var
                end
            end

            # Try to determine if variable is being used in an assignment
            (assigned, used) = check_assignment(interp, i)

            if used && !(var in usedvars)
                push!(usedvars, var)
                (var in assignedvars) || push!(putblock.args, :(put_variable($(Meta.quot(var)), $(esc(interp[i])))))
            end
            if assigned && !(var in assignedvars)
                push!(assignedvars, var)
                push!(getblock.args, Expr(:(=), esc(interp[i]), :(get_variable($(Meta.quot(var))))))
            end

            interp[i] = var
        end
    end

    # Clear `ans` and set `matlab_jl_has_ans` before we run the code
    unshift!(interp, "clear ans;\nmatlab_jl_has_ans = 0;\n")

    # Add a semicolon to the end of the last statement to suppress output
    isa(interp[end], AbstractString) && (interp[end] = rstrip(interp[end]))
    push!(interp, ";")

    # Figure out if `ans` exists in code to avoid an error if it doesn't
    push!(interp, "\nmatlab_jl_has_ans = exist('ans', 'var');")

    quote
        $(putblock)
        eval_string($(join(interp)))
        $(getblock)
        $(if !isempty(usedvars) || !isempty(assignedvars)
            # Clear variables we created
            :(eval_string($(string("clear ", join(union(usedvars, assignedvars), " "), ";"))))
        end)
        if get_variable(:matlab_jl_has_ans) != 0
            # Return ans if it was set
            get_variable(:ans)
        end
    end
end

macro mat_str(ex)
    do_mat_str(ex)
end

# Only needed for Julia 0.3
macro mat_mstr(ex)
    do_mat_str(ex)
end