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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 4c1ebac8-81b7-11eb-19fa-f704b4d84a21
html"""

Section 2.1

 Principal Component Analysis

Lecture Video

"""

# ╔═╡ cf82077a-81c2-11eb-1de2-09ed6c35d810
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add([
			Pkg.PackageSpec(name="Images", version="0.22.4"),
			Pkg.PackageSpec(name="ImageMagick", version="0.7"),
			Pkg.PackageSpec(name="PlutoUI", version="0.7"),
			Pkg.PackageSpec(name="Plots", version="1.10"),
			Pkg.PackageSpec(name="Colors", version="0.12"),
			Pkg.PackageSpec(name="ColorSchemes", version="3.10"),
			Pkg.PackageSpec(name="ForwardDiff"),
			Pkg.PackageSpec(name="LaTeXStrings"),

			])

	using PlutoUI
	using Colors, ColorSchemes, Images
	using Plots
	using LaTeXStrings

	using Statistics, LinearAlgebra  # standard libraries
end

# ╔═╡ c593a748-81b6-11eb-295a-a9800f9dec6d
PlutoUI.TableOfContents(aside=true)

# ╔═╡ deb2af50-8524-11eb-0dd4-9d799ff6d3e2
md"""
# Introduction: Understanding data

In this notebook we will start looking at more general kinds of **data**, not only images, and we'll try to extract some information from the image using statistical methods, namely [**principal component analysis**](https://en.wikipedia.org/wiki/Principal_component_analysis).

This method tries to answer the questions "which 'directions' are the most important  in the data" and "can we [reduce the dimensionality](https://en.wikipedia.org/wiki/Dimensionality_reduction) (number of useful variables) of the data"?

This can be viewed as another take on the question of finding and exploiting **structure** in data. It also leads towards ideas of **machine learning**.
"""

# ╔═╡ 2e50a070-853f-11eb-2045-b1cc43c29768
md"""
# Rank of a matrix
"""

# ╔═╡ ed7ff6b2-f863-11ea-1a59-eb242a8674e3
md"## Flags"

# ╔═╡ fed5845e-f863-11ea-2f95-c331d3c62647
md"Let's start off by recalling ideas about  **multiplication tables** or **outer products**:"

# ╔═╡ 0e1a6d80-f864-11ea-074a-5f7890180114
outer(v, w) = [x * y for x in v, y in w]

# ╔═╡ 2e497e30-f895-11ea-09f1-d7f2c1f61193
outer(1:10, 1:12)

# ╔═╡ cfdd04f8-815a-11eb-0409-79a2599c29ab
md"""
Each column is a multiple (in general not an integer multiple) of any other column; and each row is a multiple of any other row.

This is an example of a **structured matrix**, i.e. one in which we need to store *less information* to store the matrix than the full $(m \times n)$ table of numbers (even though no element is 0). For example, for an outer product we only need $m + n$ numbers, which is usually *much* less.
"""

# ╔═╡ ab3d55cc-f905-11ea-2f22-5398f3aca803
md"Some flags give simple examples of outer products, for example"

# ╔═╡ 13b6c108-f864-11ea-2447-2b0741f15c7b
flag = outer([1, 0.1, 2], ones(6))

# ╔═╡ e66b30a6-f914-11ea-2c0f-35282d45a30a
ones(6)

# ╔═╡ 43bff19e-f864-11ea-2315-0f85b532a325
show_image(flag)

# ╔═╡ 71d1b12e-f895-11ea-39df-f5c18a7766c3
flag2 = outer([1, 0.1, 2], [1, 1, 1, 3, 3, 3])

# ╔═╡ 79d2c6f4-f895-11ea-30c4-9d1102c99482
show_image(flag2)

# ╔═╡ 356267fa-815b-11eb-1c57-ad14fd6e91a7
md"""
Note that outer products are not always immediate to recognise just by looking at an image! But you should be able to recognise that there is some kind of structure.
"""

# ╔═╡ cdbe1d8e-f905-11ea-3884-efeeef386dda
md"## Matrix rank"

# ╔═╡ d9aa9af0-f865-11ea-379e-f16b452bd94c
md"""
If a matrix can be written exactly as a *single* multiplication table / outer product, we say that its **rank** is 1, and we call it a **rank-1** matrix. Similarly, if it can be written as the *sum* of *two* outer products, it has **rank 2**, etc.
"""

# ╔═╡ 2e8ae92a-f867-11ea-0219-1bdd9627c1ea
md"Let's see what a random rank-1 matrix looks like:"

# ╔═╡ 9ad13804-815c-11eb-0253-8f8baf15eee3
w = 300

# ╔═╡ 38adc490-f867-11ea-1de5-3b633aff7c97
image = outer([1; 0.4; rand(50)], rand(w));

# ╔═╡ b183b6ca-f864-11ea-0b34-4dd3f4f5e69d
show_image(image)

# ╔═╡ 946fde3c-815b-11eb-3039-db4105bc43ab
md"""
It has a characteristic checkerboard or patchwork look.
"""

# ╔═╡ ab924210-815b-11eb-07fe-411db58fbc3a
md"""
Here's a random rank-2 matrix:
"""

# ╔═╡ 74c04322-815b-11eb-2308-7b3d571cf613
begin

	image2 = outer([1; 0.4; rand(50)], rand(w)) +
	         outer(rand(52), rand(w))

	show_image(image2)
end

# ╔═╡ b5094384-815b-11eb-06fd-1f40134c6fd8
md"""
We see that it starts to look less regular.
"""

# ╔═╡ cc4f3fee-815b-11eb-2982-9b797b806b45
md"""
#### Exercise:

Make an interactive visualisation of a random rank-$n$ matrix where you can vary $n$.
"""

# ╔═╡ dc55775a-815b-11eb-15b7-7993190bffab
md"""
## Effect of noise
"""

# ╔═╡ 9cf23f9a-f864-11ea-3a08-af448aceefd8
md"""
Now what happens if we add a bit of **noise**, i.e. randomness, to a rank-1 matrix?
"""

# ╔═╡ a5b62530-f864-11ea-21e8-71ccfed487f8
noisy_image = image .+ 0.03 .* randn.();

# ╔═╡ f6713bec-815b-11eb-2fc4-6b0326a64b16
show_image(image)

# ╔═╡ 5471ddce-f867-11ea-2519-21981f5ea68b
show_image(noisy_image)

# ╔═╡ c41df86c-f865-11ea-1253-4942bbdbe9d2
md"""The noisy image now has a rank larger than 1. But visually we can see that it is "close to" the original rank-1 matrix.

Given this matrix, how could we discover that it is close to a structured, rank-1 matrix? We would like to be able to find this out and say that the matrix is close to a simple one."""

# ╔═╡ 7fca33ac-f864-11ea-2a8b-933eb382c172
md"## Images as data"

# ╔═╡ 283f5da4-f866-11ea-27d4-957ca2551b92
md"""

Now let's treat the image as a **data matrix**, so that each column of the image / matrix is a **vector** representing one observation of data. (In data science it is often the rows that correspond to observations.)

Let's try to visualize those vectors, taking just the first two rows of the image as the $x$ and $y$ coordinates of our data points:
"""

# ╔═╡ 1957f71c-f8eb-11ea-0dcf-339bfa7f96fc
show_image(image[1:2, 1:20])

# ╔═╡ 54977286-f908-11ea-166d-d1df33f38454
image[1:2, 1:20]

# ╔═╡ 7b4e90b4-f866-11ea-26b3-95efde6c650b
begin
	xx = image[1, :]
	yy = image[2, :]
end

# ╔═╡ f574ad7c-f866-11ea-0efa-d9d0602aa63b
md"# From images to data"

# ╔═╡ 8775b3fe-f866-11ea-3e6f-9732e39a3525
md"""
We would like to **visualise** this data with the `Plots.jl` package, passing it the $x$ and $y$ coordinates as arguments and plotting a point at each $(x_i, y_i)$ pair.
We obtain the following plot of the original rank-1 matrix and the noisy version:
"""

# ╔═╡ cb7e85e4-815c-11eb-3923-7b9537801bae
gr()

# ╔═╡ 1147cbda-f867-11ea-08fa-ef6ed2ae1e93
begin
	xs = noisy_image[1, :]
	ys = noisy_image[2, :]

	scatter(xs, ys, label="noisy", m=:., alpha=0.3, ms=4, ratio=1)

	scatter!(xx, yy,
			leg=:topleft, label="rank-1", ms=3, alpha=0.3,
			size=(500, 400), m=:square, c=:red,
			framestyle=:origin)

	title!("Plotting a rank-1 matrix gives a straight line!")
end

# ╔═╡ 8a611e36-f867-11ea-121f-317b7c145fe3
md"We see that the exact rank-1 matrix has columns that **lie along a line** through the origin in this representation, since they are just multiples of one another.

E.g. If $(x_1, y_1)$ and $(x_2, y_2)$ are two columns with $x_2 = cx_1$ and $y_2 = cy_1$, then $y_2 / x_2 = y_1 / x_1$, so they lie along the same line through the origin.

The approximate rank-1 matrix has columns that **lie *close to* the line**!"

# ╔═╡ f7371934-f867-11ea-3b53-d1566684585c
md"So, given the data, we want to look at it do see if it lies close to a line or not.
How can we do so in an *automatic* way?
"

# ╔═╡ 987c1f2e-f868-11ea-1125-0d8c02843ae4
md"""
## Measuring data cloud "size" -- using statistics
"""

# ╔═╡ 9e78b048-f868-11ea-192e-d903265d1eb5
md"Looking at this cloud of data points, a natural thing to do is to try to *measure* it: How wide is it, and how tall?"

# ╔═╡ 24df1f32-ec90-11ea-1f6d-03c1bfa5df8e
md"""For example, let's think about calculating the width of the cloud, i.e. the range of possible $x$-values of the data. For this purpose the $y$-values are actually irrelevant.
"""

# ╔═╡ b264f724-81bf-11eb-1052-295b81cde5fb
md"""
A natural idea would be to just scan the data and take the maximum and minimum values. However, real data often contains anomalously large values called **outliers**, which would dramatically affect this calculation. Instead we need to use a **statistical** method where we weight the data and average over all the data points. This process will hopefully be affected less by outliers, and will give a more representative idea of the size of the **bulk** of the data.
"""

# ╔═╡ 09991230-8526-11eb-04aa-fb904bd2036c
md"""
A first step in analysing data is often to **centre** it data around 0 by subtracting the mean, sometimes called "de-meaning":
"""

# ╔═╡ aec46a9b-f743-4cbd-97a7-3ef3cac78b12
begin
	xs_centered = xs .- mean(xs)
	ys_centered = ys .- mean(ys)
end

# ╔═╡ 1b8c743e-ec90-11ea-10aa-e3b94f768f82
scatter(xs_centered, ys_centered, ms=5, alpha=0.5, ratio=1, leg=false, framestyle=:origin)

# ╔═╡ eb867f18-852e-11eb-005f-e15b6d0d0d95
md"""
## Measuring a "width" of a data set
"""

# ╔═╡ f7079016-852e-11eb-3bc9-53fa0846276f
md"""
A natural way to measure the width of a data set could be to measure some kind of width *separately* in both the $x$ and $y$ directions, in other words by **projecting** the data onto one of the axes, while ignoring the other one. Let's start with the $x$ coordinates of the centred data and try to average them.

If we literally average them we will get $0$ (since we have subtracted the mean). We could ask "on average how far away from the origin is the data". This would give the [mean absolute deviation](https://en.wikipedia.org/wiki/Average_absolute_deviation):
"""

# ╔═╡ 870d3efa-f8fc-11ea-1593-1552511dcf86
begin
	scatter(xs_centered, ys_centered, ms=5, alpha=0.5, ratio=1, leg=false, framestyle=:origin)

	scatter!(xs_centered, zeros(size(xs_centered)), ms=5, alpha=0.1, ratio=1, leg=false, framestyle=:origin)

	for i in 1:length(xs_centered)
		plot!([(xs_centered[i], ys_centered[i]), (xs_centered[i], 0)], ls=:dash, c=:black, alpha=0.1)
	end

	plot!()

end

# ╔═╡ 4fb82f18-852f-11eb-278d-cf93571f4adc
mean(abs.(xs_centered))

# ╔═╡ 5fcf832c-852f-11eb-1354-792933a891a5
md"""
This is a perfectly good computational measure of distance. However, there is another measure which is easier to reason about theoretically / analytically:
"""

# ╔═╡ ef4a2a54-81bf-11eb-358b-0da2072f20c8
md"""
### Root-mean-square distance: Standard deviation
"""

# ╔═╡ f5358ce4-f86a-11ea-2989-b1f37be89183
md"""
The **standard deviation** is the **root-mean-square** distance of the centered data from the origin.

In other words, we first *square* the distances (or displacements) from the origin, then take the mean of those, giving the **variance**. However, since we have squared the original distances, this gives a quantity with units "distance$^2$", so we need to take the square root to get back to a measurable length:
"""

# ╔═╡ 2c3721da-f86b-11ea-36cf-3fe4c6622dc6
begin
	σ_x = √(mean(xs_centered.^2))   # root-mean-square distance from 0
	σ_y = √(mean(ys_centered.^2))
end

# ╔═╡ 03ab44c0-f8fd-11ea-2243-1f3580f98a65
md"This gives the following approximate extents (standard deviations) of the cloud:"

# ╔═╡ 6dec0db8-ec93-11ea-24ad-e17870ee64c2
begin
	scatter(xs_centered, ys_centered, ms=5, alpha=0.5, ratio=1, leg=false,
			framestyle=:origin)

	vline!([-2*σ_x, 2*σ_x], ls=:dash, lw=2, c=:green)
	hline!([-2*σ_y, 2*σ_y], ls=:dash, lw=2, c=:blue)

	annotate!( 2σ_x * 0.93, 0.03, text(L"2\sigma_x",  14, :green))
	annotate!(-2σ_x * 0.88, 0.03, text(L"-2\sigma_x", 14, :green))

	annotate!(0.05,  2σ_y * 1.13, text(L"2\sigma_y",  14, :blue))
	annotate!(0.06, -2σ_y * 1.14, text(L"-2\sigma_y", 14, :blue))

end

# ╔═╡ 5fab2c32-f86b-11ea-2f27-ed5feaac1fa5
md"""
We expect most (around 95%) of the data to be contained within the interval $\mu \pm 2 \sigma$, where $\mu$ is the mean and $\sigma$ is the standard deviation. (This assumes that the data is **normally distributed**, which is not actually the case for the data generated above.)
"""

# ╔═╡ ae9a2900-ec93-11ea-1ae5-0748221328fc
md"## Correlated data"

# ╔═╡ b81c9db2-ec93-11ea-0dbd-4bd0951cb2cc
md"""
However, from the figure it is clear that $x$ and $y$ are not the "correct directions" to use for this data set. It would be more natural to think about other directions: the direction in which the data set is mainly pointing (roughly, the direction in which it's longest), together with the approximately perpendicular direction in which it is  narrowest.

We need to find *from the data* which directions these are, and the extent (width) of the data cloud in those directions.

However, we cannot obtain any information about those directions by looking *separately* at $x$-coordinates and $y$-coordinates, since within the same bounding boxes that we just calculated the data can be distributed in many different ways.

Rather, the information that we need is encoded in the *relationship* between the values of $x_i$ and $y_i$ *for the points in the data set*.

For our data set, when $x$ is large and negative, $y$ is also rather negative; when $x$ is 0, $y$ is near $0$, and when $x$ is large and positive, so is $y$. We say that $x$ and $y$ are **correlated** -- literally they are mutually ("co") related, such that knowing some information about one of them allows us to predict something about the other.

For example, if I measure a new data point from the same process and I find that $x$ is around $0.25$ then I would expect $y$ to be within the range $0.05$ to $0.2$, and it would be very surprising if $y$ were -0.5.

"""

# ╔═╡ 80722856-f86d-11ea-363d-53fc5f6b8152
md"""
Although there are standard methods to calculate this correlation, we prefer to hone our **intuition** using computational thinking instead!
"""

# ╔═╡ b8fa6a1c-f86d-11ea-3d6b-2959d737254b
md"""
We want to think about different *directions*, so let's introduce an angle $\theta$ to describe the direction along which we are looking. We want to calculate the width of the cloud *along that direction*.

Effectively we are *changing coordinates* to a new coordinate, oriented along the line. To do this  requires more linear algebra than we are assuming in this course, but let's see what it looks like:
"""

# ╔═╡ 3547f296-f86f-11ea-1698-53d3c1a0bc30
md"## Rotating the axes"

# ╔═╡ 7a83101e-f871-11ea-1d87-4946162777b5
md"""By rotating the axes we can "look in different directions" and calculate the width of the data set "along that direction". What we are really doing is a perpendicular **projection** of the data onto that direction."""

# ╔═╡ e8276b4e-f86f-11ea-38be-218a72452b10
M = [xs_centered ys_centered]'

# ╔═╡ 1f373bd0-853f-11eb-0f8e-19cb7f376182
eigvals(cov(M')) .* 199

# ╔═╡ 31e4b138-84e8-11eb-36a8-8b90746fbb0f
variances = σs.^2 ./ 199

# ╔═╡ d71fdaea-f86f-11ea-1a1f-45e4d50926d3
imax = argmax(M[1, :])

# ╔═╡ 757c6808-f8fe-11ea-39bb-47e4da65113a
svdvals(M)

# ╔═╡ 1232e848-8540-11eb-089b-2185cc06f23a
M

# ╔═╡ 3b71142c-f86f-11ea-0d43-47011d00786c
p1 = begin

	scatter(M[1, :], M[2, :], ratio=1, leg=false, ms=2.5, alpha=0.5,
			framestyle=:origin)

	projected = ([cos(θ) sin(θ)] * M) .* [cos(θ) sin(θ)]'
	scatter!(projected[1, :], projected[2, :], m=:3, alpha=0.1, c=:green)

	lines_x = reduce(vcat, [M[1, i], projected[1, i], NaN] for i in 1:size(M, 2))
	lines_y = reduce(vcat, [M[2, i], projected[2, i], NaN] for i in 1:size(M, 2))

# 	for i in 1:size(M, 2)
# 		plot!([M[1, i], projected[1, i]], [M[2, i], projected[2, i]], ls=:dash, c=:black, alpha=0.1)
# 	end

	plot!(lines_x, lines_y, ls=:dash, c=:black, alpha=0.1)

	plot!([0.7 .* (-cos(θ), -sin(θ)), 0.7 .* (cos(θ), sin(θ))], lw=1, arrow=true, c=:red, alpha=0.3)
	xlims!(-0.7, 0.7)
	ylims!(-0.7, 0.7)

	scatter!([M[1, imax]], [M[2, imax]],  ms=3, alpha=1, c=:yellow)

	annotate!(0, 1.2, text("align arrow with cloud", :red, 10))
end;

# ╔═╡ 88bbe1bc-f86f-11ea-3b6b-29175ddbea04
p2 = begin
	M2 = R(θ) * M

	scatter(M2[1, :], M2[2, :],ratio=1, leg=false, ms=2.5, alpha=0.3, framestyle=:origin, size=(500, 500))

	#plot!([(-0.6, 0), (0.6, 0)], lw=3, arrow=true, c=:red, xaxis=false, yaxis=false, xticks=[], yticks=[])

	#scatter!([M2[1, imax]], [M2[2, imax]], ms=3, alpha=1, c=:yellow)

	xlims!(-0.7, 0.7)
	ylims!(-0.7, 0.7)

	scatter!(M2[1, :], zeros(size(xs_centered)), ms=3, alpha=0.1, ratio=1, leg=false, framestyle=:origin, c=:green)

	lines2_x = reduce(vcat, [M2[1, i], M2[1, i], NaN] for i in 1:size(M2, 2))
	lines2_y = reduce(vcat, [M2[2, i], 0, NaN] for i in 1:size(M2, 2))

	# for i in 1:size(M2, 2)
	# 	plot!([(M2[1, i], M2[2, i]), (M2[1, i], 0)], ls=:dash, c=:black, alpha=0.1)
	# end

	plot!(lines2_x, lines2_y, ls=:dash, c=:black, alpha=0.1)

	σ = std(M2[1, :])
	vline!([-2σ, 2σ], ls=:dash, lw=2)

# 	plot!(0.5 * [-cos(θ), cos(θ)], 0.5 * [sin(θ), -sin(θ)], c=:black, alpha=0.5, lw=1, arrow=true)
# 	plot!(0.5 * [-sin(θ), sin(θ)], 0.5 * [-cos(θ), cos(θ)], c=:black, alpha=0.5, lw=1, arrow=true)

# 	plot!(0.5 * [cos(θ), -cos(θ)], 0.5 * [-sin(θ), sin(θ)], c=:black, alpha=0.5, lw=1, arrow=true)
# 	plot!(0.5 * [sin(θ), -sin(θ)], 0.5 * [cos(θ), -cos(θ)], c=:black, alpha=0.5, lw=1, arrow=true)

	title!("σ = $(round(σ, digits=4))")

	annotate!(2σ+0.05, 0.05, text("2σ", 10, :green))
	annotate!(-2σ-0.05, 0.05, text("-2σ", 10, :green))

end;

# ╔═╡ 7cb04c9a-8358-11eb-1255-8d8c90916c37
gr()

# ╔═╡ cd9e05ee-f86f-11ea-0422-25f8329c7ef2
R(θ)= [cos(θ) sin(θ)
	  -sin(θ) cos(θ)]

# ╔═╡ 7eb51908-f906-11ea-19d2-e947d81cb743
md"In the following figure, we are rotating the axis (red arrow) around in the left panel. In the right panel we are viewing the data from the point of view of that new coordinate direction, in other words projecting onto that direction, effectively as if we rotated our head so the red vector was horizontal:"

# ╔═╡ 8b8e6b2e-8531-11eb-1ea6-637db25b28d5
p1

# ╔═╡ 4f1980ea-f86f-11ea-3df2-35cca6c961f3
md"""
degrees = $(@bind degrees Slider(0:360, default=28, show_value=true))
"""

# ╔═╡ c9da6e64-8540-11eb-3984-47fdf8be0dac
md"""
## Rotating the data
"""

# ╔═╡ f70065aa-835a-11eb-00cb-ffa27bcb486e
θ = π * degrees / 180   # radians

# ╔═╡ 2ffe7ed0-f870-11ea-06aa-390581500ca1
plot(p2)

# ╔═╡ a5cdad52-f906-11ea-0486-755a6403a367
md"""
We see that the extent of the data in the direction $\theta$ varies as $\theta$ changes. Let's plot the variance in direction $\theta$ as a function of $\theta$:
"""

# ╔═╡ 0115c974-f871-11ea-1204-054510848849
begin
	variance(θ) = var( (R(θ) * M)[1, :] )
	variance(θ::AbstractArray) = variance(θ[1])
end

# ╔═╡ 0935c870-f871-11ea-2a0b-b1b824379350
p3 = begin
	plot(0:360, variance.(range(0, 2π, length=361)), leg=false, size=(400, 200))
	scatter!([degrees], [σ^2])
	xlabel!("θ")
	ylabel!("variance in direction θ")
end

# ╔═╡ 6646abe0-835b-11eb-328a-55ca22f89c7d
σs = svdvals(M)

# ╔═╡ ef850e8e-84e7-11eb-1cb0-870c3000841d
1 ./ σs

# ╔═╡ 031a2894-84e8-11eb-1381-9b9e86f2fa0a
M

# ╔═╡ e4af4d26-f877-11ea-1de3-a9f8d389138e
md"""The direction in which the variance is **maximised** gives the most important direction: It is the direction along which the data "points", or the direction which best distinguishes different data points. This is often called the first **principal component** in statistics, or the first **singular vector** in linear algebra.

We can also now quantify *how close* the data is to lying along that single line, using the width in the *perpendicular* direction: if that width is "much smaller" than  the width in the first principal direction then the data is close to being rank 1.
"""

# ╔═╡ bf57f674-f906-11ea-08eb-9b50818a025b
md"The simplest way to maximise this function is to evaluate it everywhere and find one of the places where it takes the maximum value:"

# ╔═╡ 17e015fe-f8ff-11ea-17b4-a3aa072cd7b3
begin
	θs = 0:0.01:2π
	fs = variance.(θs)

	θmax = θs[argmax(fs)]
	θmin = θs[argmin(fs)]

	fmax = variance(θmax)
	fmin = variance(θmin)
end

# ╔═╡ 045b9b98-f8ff-11ea-0d49-5b209319e951
begin
	scatter(xs_centered, ys_centered, ms=5, alpha=0.3, ratio=1, leg=false,
			framestyle=:origin)

	plot!([(0, 0), 2*sqrt(fmax) .* (cos(θmax), sin(θmax))], arrow=true, lw=3, c=:red)
	plot!([(0, 0), 2*sqrt(fmin) .* (cos(θmin), sin(θmin))], arrow=true, lw=3, c=:red)

end

# ╔═╡ cfec1ec4-f8ff-11ea-265d-ab4844f0f739
md"""
Note that the directions that maximise and minimise variance are perpendicular. This is always the case, as shown using the Singular-Value Decomposition (SVD) in linear algebra.

There are different ways to think about this procedure. We can think of it as effectively "fitting an ellipse" to the data, where the widths of the ellipse axes show the relative importance of each direction in the data.

Alternatively we can think of it as fitting a multivariate normal distribution by finding a suitable covariance matrix.
"""

# ╔═╡ e6e900b8-f904-11ea-2a0d-953b99785553
begin
	circle = [cos.(θs) sin.(θs)]'
	stretch = [2 * sqrt(fmax) 0
				0   		2 * sqrt(fmin)]
	ellipse = R(-θmax) * stretch * circle

	plot!(ellipse[1, :], ellipse[2, :], series=:shape, alpha=0.4, fill=true, c=:orange)
end

# ╔═╡ 301f4d06-8162-11eb-1cd6-31dd8da164b6
md"""
Note also that an ellipse is the image of the unit circle under a *linear* transformation. We are effectively learning what the best linear transformation is that transforms the unit circle into our data cloud.
"""

# ╔═╡ aaff88e8-f877-11ea-1527-ff4d3db663db
md"## Higher dimensions"

# ╔═╡ aefa84de-f877-11ea-3e26-678008e9739e
md"""
Can we generalise this to dimensions higher than 2?
Let's think about 3D. If we take columns of the first *three* rows of the original image, we have vectors in 3D.

If we plot these vectors in 3D, a rank-1 matrix will give a straight line in 3D, while a rank-2 matrix will give a **plane** in 3D. Rank-2 + noise gives a noisy cloud lying close to a plane.

Similarly to what we did above, we need to calculate the ellipsoid that best fits the data. The widths of the axes of the ellipsoid tell us how close to being a line or a plane (rank-1 or rank-2) the data is.
"""

# ╔═╡ 0bd9358e-f879-11ea-2c83-ed4e7bf9d903
md"In more than 3D we can no longer visualise the data, but the same idea applies. The calculations are done using the SVD.

If the widths of the ellipsoid in some directions are very small, we can ignore those directions and hence reduce the dimensionality of the data, by changing coordinates to the principal components."

# ╔═╡ 9960e1c2-85a0-11eb-3d68-cd3a07a9c0b5
md"""
## What is the Singular-Value Decomposition (SVD)?
"""

# ╔═╡ dc4cca88-85a0-11eb-2791-d7610a610e36
md"""
The **Singular-Value Decomposition (SVD)** is a way of writing any matrix in terms of simpler matrices. Thinking in terms of linear transformations, *any* linear transformation $T$ has the same effect as a sequence of three simple transformations:

T = rotation₂ ∘ stretch ∘ rotation₁

In terms of matrices, for *any* matrix $M$ it is possible to write

$$M = U \, \Sigma \, V^\text{T}$$

where $U$ and $V$ are **orthogonal** matrices, i.e. they satisfy $U U^\text{T} = I$ and $V V^\text{T} = I$, and $\Sigma$ is a diagonal matrix.

Orthogonal matrices have determinant $\pm 1$, so they leave areas unchanged. They are rotations, possibly combined with reflections.

There are algorithms from numerical linear algebra to calculate this decomposition. In Julia we call the `svd` function, e.g.
"""

# ╔═╡ 453689c2-85a2-11eb-2cbc-7d6476b42f2f
let
	M = [2 1
		 1 1]

	svd(M)
end

# ╔═╡ 91b77acc-85a2-11eb-19bb-bd2bdd0c1a68
md"""
Let's look at the action of the matrix on the unit disc. To generate points in the unit disc we generate points in $[-1, 1]^2$ and *reject* those lying outside:
"""

# ╔═╡ f92c75f6-85a3-11eb-1689-23aeaa3daeb7
begin
	unit_disc = [ (-1.0 .+ 2.0 .* rand(2)) for i in 1:2000 ]
	unit_disc = reduce(hcat, [x for x in unit_disc if x[1]^2 + x[2]^2 < 1])
end

# ╔═╡ 03069da6-85a4-11eb-2ac5-87b767846550
scatter(unit_disc[1, :], unit_disc[2, :], ratio=1, leg=false, alpha=0.5, ms=3)

# ╔═╡ 40b87cbe-85a4-11eb-30f8-cf7b5e79c19a
pp1 = begin
	scatter(unit_disc[1, :], unit_disc[2, :], ratio=1, leg=false, alpha=0.5, title="stretch + rotate")
	result =  [1 + tt  tt; tt  1] * unit_disc
	scatter!(result[1, :], result[2, :], alpha=0.2)

	ylims!(-3, 3)
	xlims!(-3, 3)
end;

# ╔═╡ 28a7d6dc-85a5-11eb-0e4b-b7e9b4c592ed
pp2 = begin
	UU, Sigma, VV = svd([1 + tt  tt; tt  1])
	scatter(unit_disc[1, :], unit_disc[2, :], ratio=1, leg=false, alpha=0.5, title="stretch")

	result2 = Diagonal(Sigma) * unit_disc

	scatter!(result2[1, :], result2[2, :], alpha=0.2)

	ylims!(-3, 3)
	xlims!(-3, 3)

end;

# ╔═╡ 1647a126-85a4-11eb-3923-5f5a6f703403
md"""
t = $(@bind tt Slider(0:0.01:1, show_value=true))
"""

# ╔═╡ 6ec7f980-85a5-11eb-12fc-cb132db28d83
plot(pp2, pp1)

# ╔═╡ 92a2827e-84e9-11eb-1e85-6f49b1da7277
md"""
## Rotations in 300 dimensions
"""

# ╔═╡ e84adec2-84e8-11eb-2157-dd491588ccf0
md"""
We have been thinking of 300 points in 2 dimensions. Instead, we could think about 2 points **in 300 dimensions**. In some sense, this is what
"transpose really does"!
"""

# ╔═╡ 571b88ac-85a6-11eb-2887-a368a19bce4d
md"""
The $V$ in the SVD is another rotation that we "cannot see" in the above pictures. In our case it is a rotation in 300 dimensions! But in fact we can visualise it as follows, where we multiply our data by a $300 \times 300$ orthogonal matrix!

We first take a random *anti-symmetric* matrix, i.e. one for which $M = -M^\text{t}$.
We then turn that into an orthogonal matrix by taking the so-called **matrix exponential**.
"""

# ╔═╡ 12010a58-84eb-11eb-106f-cb4e3e0c879b
begin
	dim = size(M, 2)
	anti_symmetric = randn(dim, dim)
	anti_symmetric -= anti_symmetric'
end

# ╔═╡ 696d2768-84eb-11eb-39e0-612e074a2c27
@bind t Slider(0:0.0002:1, show_value=true, default=0.0)

# ╔═╡ 7b7b5128-84eb-11eb-3974-9b4c08fab8bb
begin
	M_rotated = M * exp(t * anti_symmetric)
	scatter(M_rotated[1, :], M_rotated[2, :], leg=false, alpha=0.5, )
	scatter!(M[1, :], M[2, :], alpha=0.5)

	ylims!(-0.3, 0.3)
	xlims!(-0.6, 0.6)
end

# ╔═╡ 805b9616-85a7-11eb-22e8-db8ee67071ae
md"""
We see that the data rotates around in 300 dimensions, but always is projected to the *same* ellipse.
"""

# ╔═╡ 90656ce6-84fb-11eb-1aac-4bd7747613db
U, Σ, V = svd(M, full=true)

# ╔═╡ aec542a2-84fb-11eb-322c-27fc2c45f6ef
M18 = M * V

# ╔═╡ b55dcfd2-84fb-11eb-1766-17dc8b7a17d0
scatter(M18[1, :], M18[2, :], alpha=0.5, leg=false, ratio=1, xlim=(-5,5))

# ╔═╡ 1cf3e098-f864-11ea-3f3a-c53017b73490
md"#### Appendix"

# ╔═╡ 2917943c-f864-11ea-3ee6-db952ca7cd67
begin
	show_image(M) = get.(Ref(ColorSchemes.rainbow), M ./ maximum(M))
	show_image(x::AbstractVector) = show_image(x')
end

# ╔═╡ 72bb11b0-f88f-11ea-0e55-b1108300f854
loss(M1, M2) = sum( (M1[i] - M2[i])^2 for i in 1:length(M1) if !ismissing(M2[i]) )

# ╔═╡ feeeb24a-f88f-11ea-287f-219e53615f32
function split_up(v, m, n)
	return v[1:m], v[m+1:m+n], v[m+n+1:2m+n], v[2m+n+1:2m+2n]
end

# ╔═╡ 0bcc8852-f890-11ea-3715-11cbead7f636
function ff(v, m, n)
	v1, w1, v2, w2 = split_up(v, m, n)

	loss(outer(v1, w1) + outer(v2, w2), M3)
end

# ╔═╡ 7040dc72-f893-11ea-3d22-4fbd452faa41
ff2(v) = ff(v, m, n)

# ╔═╡ 1dbcf15a-f890-11ea-008c-8935edfbdb1c
ff(rand(total), m, n)

# ╔═╡ Cell order:
# ╟─4c1ebac8-81b7-11eb-19fa-f704b4d84a21
# ╠═cf82077a-81c2-11eb-1de2-09ed6c35d810
# ╠═c593a748-81b6-11eb-295a-a9800f9dec6d
# ╟─deb2af50-8524-11eb-0dd4-9d799ff6d3e2
# ╟─2e50a070-853f-11eb-2045-b1cc43c29768
# ╟─ed7ff6b2-f863-11ea-1a59-eb242a8674e3
# ╟─fed5845e-f863-11ea-2f95-c331d3c62647
# ╠═0e1a6d80-f864-11ea-074a-5f7890180114
# ╠═2e497e30-f895-11ea-09f1-d7f2c1f61193
# ╟─cfdd04f8-815a-11eb-0409-79a2599c29ab
# ╟─ab3d55cc-f905-11ea-2f22-5398f3aca803
# ╠═13b6c108-f864-11ea-2447-2b0741f15c7b
# ╠═e66b30a6-f914-11ea-2c0f-35282d45a30a
# ╠═43bff19e-f864-11ea-2315-0f85b532a325
# ╠═71d1b12e-f895-11ea-39df-f5c18a7766c3
# ╠═79d2c6f4-f895-11ea-30c4-9d1102c99482
# ╟─356267fa-815b-11eb-1c57-ad14fd6e91a7
# ╟─cdbe1d8e-f905-11ea-3884-efeeef386dda
# ╟─d9aa9af0-f865-11ea-379e-f16b452bd94c
# ╟─2e8ae92a-f867-11ea-0219-1bdd9627c1ea
# ╠═9ad13804-815c-11eb-0253-8f8baf15eee3
# ╠═38adc490-f867-11ea-1de5-3b633aff7c97
# ╠═b183b6ca-f864-11ea-0b34-4dd3f4f5e69d
# ╟─946fde3c-815b-11eb-3039-db4105bc43ab
# ╟─ab924210-815b-11eb-07fe-411db58fbc3a
# ╠═74c04322-815b-11eb-2308-7b3d571cf613
# ╟─b5094384-815b-11eb-06fd-1f40134c6fd8
# ╟─cc4f3fee-815b-11eb-2982-9b797b806b45
# ╟─dc55775a-815b-11eb-15b7-7993190bffab
# ╟─9cf23f9a-f864-11ea-3a08-af448aceefd8
# ╠═a5b62530-f864-11ea-21e8-71ccfed487f8
# ╠═f6713bec-815b-11eb-2fc4-6b0326a64b16
# ╠═5471ddce-f867-11ea-2519-21981f5ea68b
# ╟─c41df86c-f865-11ea-1253-4942bbdbe9d2
# ╟─7fca33ac-f864-11ea-2a8b-933eb382c172
# ╟─283f5da4-f866-11ea-27d4-957ca2551b92
# ╠═1957f71c-f8eb-11ea-0dcf-339bfa7f96fc
# ╠═54977286-f908-11ea-166d-d1df33f38454
# ╠═7b4e90b4-f866-11ea-26b3-95efde6c650b
# ╟─f574ad7c-f866-11ea-0efa-d9d0602aa63b
# ╟─8775b3fe-f866-11ea-3e6f-9732e39a3525
# ╠═cb7e85e4-815c-11eb-3923-7b9537801bae
# ╠═1147cbda-f867-11ea-08fa-ef6ed2ae1e93
# ╟─8a611e36-f867-11ea-121f-317b7c145fe3
# ╟─f7371934-f867-11ea-3b53-d1566684585c
# ╟─987c1f2e-f868-11ea-1125-0d8c02843ae4
# ╟─9e78b048-f868-11ea-192e-d903265d1eb5
# ╟─24df1f32-ec90-11ea-1f6d-03c1bfa5df8e
# ╟─b264f724-81bf-11eb-1052-295b81cde5fb
# ╟─09991230-8526-11eb-04aa-fb904bd2036c
# ╠═aec46a9b-f743-4cbd-97a7-3ef3cac78b12
# ╠═1b8c743e-ec90-11ea-10aa-e3b94f768f82
# ╟─eb867f18-852e-11eb-005f-e15b6d0d0d95
# ╟─f7079016-852e-11eb-3bc9-53fa0846276f
# ╟─870d3efa-f8fc-11ea-1593-1552511dcf86
# ╠═4fb82f18-852f-11eb-278d-cf93571f4adc
# ╟─5fcf832c-852f-11eb-1354-792933a891a5
# ╟─ef4a2a54-81bf-11eb-358b-0da2072f20c8
# ╟─f5358ce4-f86a-11ea-2989-b1f37be89183
# ╠═2c3721da-f86b-11ea-36cf-3fe4c6622dc6
# ╟─03ab44c0-f8fd-11ea-2243-1f3580f98a65
# ╟─6dec0db8-ec93-11ea-24ad-e17870ee64c2
# ╟─5fab2c32-f86b-11ea-2f27-ed5feaac1fa5
# ╟─ae9a2900-ec93-11ea-1ae5-0748221328fc
# ╟─b81c9db2-ec93-11ea-0dbd-4bd0951cb2cc
# ╟─80722856-f86d-11ea-363d-53fc5f6b8152
# ╟─b8fa6a1c-f86d-11ea-3d6b-2959d737254b
# ╟─3547f296-f86f-11ea-1698-53d3c1a0bc30
# ╟─7a83101e-f871-11ea-1d87-4946162777b5
# ╠═e8276b4e-f86f-11ea-38be-218a72452b10
# ╠═1f373bd0-853f-11eb-0f8e-19cb7f376182
# ╠═31e4b138-84e8-11eb-36a8-8b90746fbb0f
# ╟─d71fdaea-f86f-11ea-1a1f-45e4d50926d3
# ╠═757c6808-f8fe-11ea-39bb-47e4da65113a
# ╠═1232e848-8540-11eb-089b-2185cc06f23a
# ╠═3b71142c-f86f-11ea-0d43-47011d00786c
# ╠═88bbe1bc-f86f-11ea-3b6b-29175ddbea04
# ╠═7cb04c9a-8358-11eb-1255-8d8c90916c37
# ╟─cd9e05ee-f86f-11ea-0422-25f8329c7ef2
# ╟─7eb51908-f906-11ea-19d2-e947d81cb743
# ╠═8b8e6b2e-8531-11eb-1ea6-637db25b28d5
# ╟─4f1980ea-f86f-11ea-3df2-35cca6c961f3
# ╟─c9da6e64-8540-11eb-3984-47fdf8be0dac
# ╠═f70065aa-835a-11eb-00cb-ffa27bcb486e
# ╠═2ffe7ed0-f870-11ea-06aa-390581500ca1
# ╟─a5cdad52-f906-11ea-0486-755a6403a367
# ╟─0115c974-f871-11ea-1204-054510848849
# ╠═0935c870-f871-11ea-2a0b-b1b824379350
# ╠═6646abe0-835b-11eb-328a-55ca22f89c7d
# ╠═ef850e8e-84e7-11eb-1cb0-870c3000841d
# ╠═031a2894-84e8-11eb-1381-9b9e86f2fa0a
# ╟─e4af4d26-f877-11ea-1de3-a9f8d389138e
# ╟─bf57f674-f906-11ea-08eb-9b50818a025b
# ╠═17e015fe-f8ff-11ea-17b4-a3aa072cd7b3
# ╟─045b9b98-f8ff-11ea-0d49-5b209319e951
# ╟─cfec1ec4-f8ff-11ea-265d-ab4844f0f739
# ╠═e6e900b8-f904-11ea-2a0d-953b99785553
# ╟─301f4d06-8162-11eb-1cd6-31dd8da164b6
# ╟─aaff88e8-f877-11ea-1527-ff4d3db663db
# ╟─aefa84de-f877-11ea-3e26-678008e9739e
# ╟─0bd9358e-f879-11ea-2c83-ed4e7bf9d903
# ╟─9960e1c2-85a0-11eb-3d68-cd3a07a9c0b5
# ╟─dc4cca88-85a0-11eb-2791-d7610a610e36
# ╠═453689c2-85a2-11eb-2cbc-7d6476b42f2f
# ╟─91b77acc-85a2-11eb-19bb-bd2bdd0c1a68
# ╠═f92c75f6-85a3-11eb-1689-23aeaa3daeb7
# ╠═03069da6-85a4-11eb-2ac5-87b767846550
# ╠═40b87cbe-85a4-11eb-30f8-cf7b5e79c19a
# ╠═28a7d6dc-85a5-11eb-0e4b-b7e9b4c592ed
# ╟─1647a126-85a4-11eb-3923-5f5a6f703403
# ╠═6ec7f980-85a5-11eb-12fc-cb132db28d83
# ╟─92a2827e-84e9-11eb-1e85-6f49b1da7277
# ╟─e84adec2-84e8-11eb-2157-dd491588ccf0
# ╟─571b88ac-85a6-11eb-2887-a368a19bce4d
# ╠═12010a58-84eb-11eb-106f-cb4e3e0c879b
# ╠═696d2768-84eb-11eb-39e0-612e074a2c27
# ╠═7b7b5128-84eb-11eb-3974-9b4c08fab8bb
# ╟─805b9616-85a7-11eb-22e8-db8ee67071ae
# ╠═90656ce6-84fb-11eb-1aac-4bd7747613db
# ╠═aec542a2-84fb-11eb-322c-27fc2c45f6ef
# ╠═b55dcfd2-84fb-11eb-1766-17dc8b7a17d0
# ╟─1cf3e098-f864-11ea-3f3a-c53017b73490
# ╠═2917943c-f864-11ea-3ee6-db952ca7cd67
# ╠═72bb11b0-f88f-11ea-0e55-b1108300f854
# ╠═feeeb24a-f88f-11ea-287f-219e53615f32
# ╠═0bcc8852-f890-11ea-3715-11cbead7f636
# ╠═7040dc72-f893-11ea-3d22-4fbd452faa41
# ╠═1dbcf15a-f890-11ea-008c-8935edfbdb1c


text:
# singleton for indicating if no default arguments are present
struct NoDefault end
const NO_DEFAULT = NoDefault()

"""
    @addlogprob!(ex)

Add the result of the evaluation of `ex` to the joint log probability.

# Examples

This macro allows you to [include arbitrary terms in the likelihood](https://github.com/TuringLang/Turing.jl/issues/1332)

```jldoctest; setup = :(using Distributions)
julia> myloglikelihood(x, μ) = loglikelihood(Normal(μ, 1), x);

julia> @model function demo(x)
           μ ~ Normal()
           @addlogprob! myloglikelihood(x, μ)
       end;

julia> x = [1.3, -2.1];

julia> loglikelihood(demo(x), (μ=0.2,)) ≈ myloglikelihood(x, 0.2)
true
```

and to [reject samples](https://github.com/TuringLang/Turing.jl/issues/1328):

```jldoctest; setup = :(using Distributions, LinearAlgebra)
julia> @model function demo(x)
           m ~ MvNormal(zero(x), I)
           if dot(m, x) < 0
               @addlogprob! -Inf
               # Exit the model evaluation early
               return
           end
           x ~ MvNormal(m, I)
           return
       end;

julia> logjoint(demo([-2.1]), (m=[0.2],)) == -Inf
true
```

!!! note
    The `@addlogprob!` macro increases the accumulated log probability regardless of the evaluation context,
    i.e., regardless of whether you evaluate the log prior, the log likelihood or the log joint density.
    If you would like to avoid this behaviour you should check the evaluation context.
    It can be accessed with the internal variable `__context__`.
    For instance, in the following example the log density is not accumulated when only the log prior is computed:
    ```jldoctest; setup = :(using Distributions)
    julia> myloglikelihood(x, μ) = loglikelihood(Normal(μ, 1), x);

    julia> @model function demo(x)
               μ ~ Normal()
               if DynamicPPL.leafcontext(__context__) !== PriorContext()
                   @addlogprob! myloglikelihood(x, μ)
               end
           end;

    julia> x = [1.3, -2.1];

    julia> logprior(demo(x), (μ=0.2,)) ≈ logpdf(Normal(), 0.2)
    true

    julia> loglikelihood(demo(x), (μ=0.2,)) ≈ myloglikelihood(x, 0.2)
    true
    ```
"""
macro addlogprob!(ex)
    return quote
        $(esc(:(__varinfo__))) = acclogp!!($(esc(:(__varinfo__))), $(esc(ex)))
    end
end

"""
    addargnames!(args)

Adds names to unnamed arguments in `args`.

The names are generated with `gensym(:arg)` to avoid conflicts with other variable names.

# Examples

```jldoctest; filter = r"var\\"##arg#[0-9]+\\""
julia> args = :(f(x::Int, y, ::Type{T}=Float64)).args[2:end]
3-element Vector{Any}:
 :(x::Int)
 :y
 :($(Expr(:kw, :(::Type{T}), :Float64)))

julia> DynamicPPL.addargnames!(args)

julia> args
3-element Vector{Any}:
 :(x::Int)
 :y
 :($(Expr(:kw, :(var"##arg#301"::Type{T}), :Float64)))
```
"""
function addargnames!(args)
    if isempty(args)
        return nothing
    end

    @inbounds for i in eachindex(args)
        arg = args[i]
        if MacroTools.@capture(arg, ::T_)
            args[i] = Expr(:(::), gensym(:arg), T)
        elseif MacroTools.@capture(arg, ::T_ = val_)
            args[i] = Expr(:kw, Expr(:(::), gensym(:arg), T), val)
        end
    end

    return nothing
end

"""
    getargs_dottilde(x)

Return the arguments `L` and `R`, if `x` is an expression of the form `L .~ R` or
`(~).(L, R)`, or `nothing` otherwise.
"""
getargs_dottilde(x) = nothing
function getargs_dottilde(expr::Expr)
    return MacroTools.@match expr begin
        (.~)(L_, R_) => (L, R)
        (~).(L_, R_) => (L, R)
        # Julia 1.6: see https://github.com/TuringLang/Turing.jl/issues/1525
        (L_ .~ R_) => (L, R)
        x_ => nothing
    end
end

"""
    getargs_tilde(x)

Return the arguments `L` and `R`, if `x` is an expression of the form `L ~ R`, or `nothing`
otherwise.
"""
getargs_tilde(x) = nothing
function getargs_tilde(expr::Expr)
    return MacroTools.@match expr begin
        (~)(L_, R_) => (L, R)
        x_ => nothing
    end
end

"""
    getargs_assignment(x)

Return the arguments `L` and `R`, if `x` is an expression of the form `L = R`, or `nothing`
otherwise.
"""
getargs_assignment(x) = nothing
function getargs_assignment(expr::Expr)
    return MacroTools.@match expr begin
        (L_ = R_) => (L, R)
        x_ => nothing
    end
end

function to_namedtuple_expr(syms, vals=syms)
    length(syms) == 0 && return :(NamedTuple())

    names_expr = Expr(:tuple, QuoteNode.(syms)...)
    vals_expr = Expr(:tuple, vals...)
    return :(NamedTuple{$names_expr}($vals_expr))
end

#####################################################
# Helper functions for vectorize/reconstruct values #
#####################################################

vectorize(d::UnivariateDistribution, r::Real) = [r]
vectorize(d::MultivariateDistribution, r::AbstractVector{<:Real}) = copy(r)
vectorize(d::MatrixDistribution, r::AbstractMatrix{<:Real}) = copy(vec(r))

# NOTE:
# We cannot use reconstruct{T} because val is always Vector{Real} then T will be Real.
# However here we would like the result to be specifric type, e.g. Array{Dual{4,Float64}, 2},
# otherwise we will have error for MatrixDistribution.
# Note this is not the case for MultivariateDistribution so I guess this might be lack of
# support for some types related to matrices (like PDMat).
reconstruct(d::Distribution, val::AbstractVector) = reconstruct(size(d), val)
reconstruct(::Tuple{}, val::AbstractVector) = val[1]
reconstruct(s::NTuple{1}, val::AbstractVector) = copy(val)
reconstruct(s::NTuple{2}, val::AbstractVector) = reshape(copy(val), s)
function reconstruct!(r, d::Distribution, val::AbstractVector)
    return reconstruct!(r, d, val)
end
function reconstruct!(r, d::MultivariateDistribution, val::AbstractVector)
    r .= val
    return r
end
function reconstruct(d::Distribution, val::AbstractVector, n::Int)
    return reconstruct(size(d), val, n)
end
function reconstruct(::Tuple{}, val::AbstractVector, n::Int)
    return copy(val)
end
function reconstruct(s::NTuple{1}, val::AbstractVector, n::Int)
    return copy(reshape(val, s[1], n))
end
function reconstruct(s::NTuple{2}, val::AbstractVector, n::Int)
    tmp = reshape(val, s..., n)
    orig = [tmp[:, :, i] for i in 1:n]
    return orig
end
function reconstruct!(r, d::Distribution, val::AbstractVector, n::Int)
    return reconstruct!(r, d, val, n)
end
function reconstruct!(r, d::MultivariateDistribution, val::AbstractVector, n::Int)
    r .= val
    return r
end

# Uniform random numbers with range 4 for robust initializations
# Reference: https://mc-stan.org/docs/2_19/reference-manual/initialization.html
randrealuni(rng::Random.AbstractRNG) = 4 * rand(rng) - 2
randrealuni(rng::Random.AbstractRNG, args...) = 4 .* rand(rng, args...) .- 2

const Transformable = Union{
    PositiveDistribution,
    UnitDistribution,
    TransformDistribution,
    SimplexDistribution,
    PDMatDistribution,
}
istransformable(dist) = false
istransformable(::Transformable) = true

#################################
# Single-sample initialisations #
#################################

inittrans(rng, dist::UnivariateDistribution) = invlink(dist, randrealuni(rng))
function inittrans(rng, dist::MultivariateDistribution)
    return invlink(dist, randrealuni(rng, size(dist)[1]))
end
inittrans(rng, dist::MatrixDistribution) = invlink(dist, randrealuni(rng, size(dist)...))

################################
# Multi-sample initialisations #
################################

inittrans(rng, dist::UnivariateDistribution, n::Int) = invlink(dist, randrealuni(rng, n))
function inittrans(rng, dist::MultivariateDistribution, n::Int)
    return invlink(dist, randrealuni(rng, size(dist)[1], n))
end
function inittrans(rng, dist::MatrixDistribution, n::Int)
    return invlink(dist, [randrealuni(rng, size(dist)...) for _ in 1:n])
end

#######################
# Convenience methods #
#######################
collectmaybe(x) = x
collectmaybe(x::Base.AbstractSet) = collect(x)

#######################
# BangBang.jl related #
#######################
function set!!(obj, lens::Setfield.Lens, value)
    lensmut = BangBang.prefermutation(lens)
    return Setfield.set(obj, lensmut, value)
end
function set!!(obj, vn::VarName{sym}, value) where {sym}
    lens = BangBang.prefermutation(Setfield.PropertyLens{sym}() ∘ AbstractPPL.getlens(vn))
    return Setfield.set(obj, lens, value)
end

#############################
# AbstractPPL.jl extensions #
#############################
# This is preferable to `haskey` because the order of arguments is different, and
# we're more likely to specialize on the key in these settings rather than the container.
# TODO: I'm not sure about this name.
"""
    canview(lens, container)

Return `true` if `lens` can be used to view `container`, and `false` otherwise.

# Examples
```jldoctest; setup=:(using Setfield; using DynamicPPL: canview)
julia> canview(@lens(_.a), (a = 1.0, ))
true

julia> canview(@lens(_.a), (b = 1.0, )) # property `a` does not exist
false

julia> canview(@lens(_.a[1]), (a = [1.0, 2.0], ))
true

julia> canview(@lens(_.a[3]), (a = [1.0, 2.0], )) # out of bounds
false
```
"""
canview(lens, container) = false
canview(::Setfield.IdentityLens, _) = true
function canview(lens::Setfield.PropertyLens{field}, x) where {field}
    return hasproperty(x, field)
end

# `IndexLens`: only relevant if `x` supports indexing.
canview(lens::Setfield.IndexLens, x) = false
canview(lens::Setfield.IndexLens, x::AbstractArray) = checkbounds(Bool, x, lens.indices...)

# `ComposedLens`: check that we can view `.outer` and `.inner`, but using
# value extracted using `.outer`.
function canview(lens::Setfield.ComposedLens, x)
    return canview(lens.outer, x) && canview(lens.inner, get(x, lens.outer))
end

"""
    parent(vn::VarName)

Return the parent `VarName`.

# Examples
```julia-repl; setup=:(using DynamicPPL: parent)
julia> parent(@varname(x.a[1]))
x.a

julia> (parent ∘ parent)(@varname(x.a[1]))
x

julia> (parent ∘ parent ∘ parent)(@varname(x.a[1]))
x
```
"""
function parent(vn::VarName)
    p = parent(getlens(vn))
    return p === nothing ? VarName(vn, Setfield.IdentityLens()) : VarName(vn, p)
end

"""
    parent(lens::Setfield.Lens)

Return the parent lens. If `lens` doesn't have a parent,
`nothing` is returned.

See also: [`parent_and_child`].

# Examples
```jldoctest; setup=:(using Setfield; using DynamicPPL: parent)
julia> parent(@lens(_.a[1]))
(@lens _.a)

julia> # Parent of lens without parents results in `nothing`.
       (parent ∘ parent)(@lens(_.a[1])) === nothing
true
```
"""
parent(lens::Setfield.Lens) = first(parent_and_child(lens))

"""
    parent_and_child(lens::Setfield.Lens)

Return a 2-tuple of lenses `(parent, child)` where `parent` is the
parent lens of `lens` and `child` is the child lens of `lens`.

If `lens` does not have a parent, we return `(nothing, lens)`.

See also: [`parent`].

# Examples
```jldoctest; setup=:(using Setfield; using DynamicPPL: parent_and_child)
julia> parent_and_child(@lens(_.a[1]))
((@lens _.a), (@lens _[1]))

julia> parent_and_child(@lens(_.a))
(nothing, (@lens _.a))
```
"""
parent_and_child(lens::Setfield.Lens) = (nothing, lens)
function parent_and_child(lens::Setfield.ComposedLens)
    p, child = parent_and_child(lens.inner)
    parent = p === nothing ? lens.outer : lens.outer ∘ p
    return parent, child
end

"""
    splitlens(condition, lens)

Return a 3-tuple `(parent, child, issuccess)` where, if `issuccess` is `true`,
`parent` is a lens such that `condition(parent)` is `true` and `parent ∘ child == lens`.

If `issuccess` is `false`, then no such split could be found.

# Examples
```jldoctest; setup=:(using Setfield; using DynamicPPL: splitlens)
julia> p, c, issucesss = splitlens(@lens(_.a[1])) do parent
           # Succeeds!
           parent == @lens(_.a)
       end
((@lens _.a), (@lens _[1]), true)

julia> p ∘ c
(@lens _.a[1])

julia> splitlens(@lens(_.a[1])) do parent
           # Fails!
           parent == @lens(_.b)
       end
(nothing, (@lens _.a[1]), false)
```
"""
function splitlens(condition, lens)
    current_parent, current_child = parent_and_child(lens)
    # We stop if either a) `condition` is satisfied, or b) we reached the root.
    while !condition(current_parent) && current_parent !== nothing
        current_parent, c = parent_and_child(current_parent)
        current_child = c ∘ current_child
    end

    return current_parent, current_child, condition(current_parent)
end

ext: jl

max_stars_repo_name: TuringLang/MicroPPL

max_issues_repo_name: TuringLang/MicroPPL

max_forks_repo_name: TuringLang/MicroPPL


text:
import Base: ==

# Rudimentary fixed width string. Nothing fancy just bare minimum, for
# lookup fast string manipulations. Unfortunately `nextind` and other string functions
# are notoriously slow.
struct String32 <: AbstractString
    s::Vector{Char}
    l::Int
end

function String32(s::S) where {S <: AbstractString}
    x = Vector{Char}(s)
    return String32(x, length(x))
end

String32(s::String32) = s

Base.:length(s::String32) = s.l
==(s1::String32, s2::String32) = (s1.s == s2.s) & (s1.l == s2.l)

function Base.:iterate(s::String32, state::Integer = 1)
    state > length(s.s) && return nothing
    return (s.s[state], state + 1)
end

Base.:ncodeunits(s::String32) = s.l
Base.:isvalid(s::String32, i::Integer) = true
Base.:lowercase(x::String32) = String32(lowercase.(x.s), x.l)
Base.:getindex(x::String32, i::Integer) = x.s[i]

struct DictWord{T <: Integer}
    word::String32
    count::T
end

DictWord(word::String, count::T) where {T <: Integer} = DictWord{T}(String32(word), count)

==(dw1::DictWord{T}, dw2::DictWord{T}) where {T <: Integer} = (dw1.word == dw2.word) & (dw1.count == dw2.count)
Base.:length(dw::DictWord) = length(dw.word)

mutable struct LookupOptions
    include_unknown::Bool
    ignore_token::Union{Nothing,AbstractString,Regex,AbstractChar}
    transfer_casing::Bool
    verbosity::Verbosity
    max_edit_distance::Int
end

mutable struct SymSpell{T <: Integer, K <: Unsigned}
    words::Vector{DictWord{T}}
    below_threshold_words::Dict{String, T}
    deletes::Dict{String, Vector{K}}
    words_idx::Dict{String, K}

    max_dictionary_edit_distance::Int
    prefix_length::Int
    count_threshold::Int
    max_length::Int

    lookup::LookupOptions
end

SymSpell(; max_dictionary_edit_distance = 2, prefix_length = 7, count_threshold = 1) =
    SymSpell{Int, UInt32}(Vector(), Dict(), Dict(), Dict(),
        max_dictionary_edit_distance, prefix_length, count_threshold, 0,
        LookupOptions(false, nothing, false, VerbosityTOP, max_dictionary_edit_distance))

function SymSpell(path; sep = " ", max_dictionary_edit_distance = 2, prefix_length = 7, count_threshold = 1)
    d = SymSpell(max_dictionary_edit_distance = max_dictionary_edit_distance,
                    prefix_length = prefix_length, count_threshold = count_threshold)
    update!(d, path, sep = sep)
end

# TODO: add support for DataFrames and CSV
function update!(dict::SymSpell{T, K}, path; sep = " ") where {T, K}
    for line in eachline(path)
        line_parts = strip.(split(line, sep))
        if length(line_parts) >= 2
            push!(dict, string(line_parts[1]), parse(T, strip(line_parts[2])))
        end
    end

    # this sorting increases chance of early stop on TOP/CLOSEST verbosity level
    for v in values(dict.deletes)
        sort!(v, by = idx -> length(dict.words[idx]))
    end

    dict
end

"""
    edits!(delete_words, word, edit_distance, max_dictionary_distance)

Inexpensive and language independent: only deletes, no transposes + replaces + inserts
replaces and inserts are expensive and language dependent
"""
function edits!(delete_words, word, edit_distance, max_dictionary_edit_distance)
    edit_distance += 1
    word_len = length(word)
    if word_len > 1
        let i = 0
            while i < lastindex(word)
                delete = word[1:i] * word[nextind(word, nextind(word, i)):end]
                i = nextind(word, i)
                if !(delete in delete_words)
                    push!(delete_words, delete)
                    # recursion, if maximum edit distance not yet reached
                    if edit_distance < max_dictionary_edit_distance
                        edits!(delete_words, delete, edit_distance, max_dictionary_edit_distance)
                    end
                end
            end
        end
    end

    return delete_words
end

function edits_prefix(key::S, max_dictionary_edit_distance, prefix_length) where {S <: AbstractString}
    hash_set = Set{S}()

    if length(key) <= max_dictionary_edit_distance
        push!(hash_set, "")
    end
    if length(key) > prefix_length
        key = key[1:nextind(key, 0, prefix_length)]
    end
    push!(hash_set, key)
    edits!(hash_set, key, 0, max_dictionary_edit_distance)
end

function Base.:push!(dict::SymSpell{T, K}, key, cnt = 1) where {T, K}
    if cnt < 0
        if dict.count_threshold > 0 return false end
        cnt = 0
    end

    # look first in below threshold words, update count, and allow
    # promotion to correct spelling word if count reaches threshold
    # threshold must be >1 for there to be the possibility of low
    # threshold words
    if dict.count_threshold > 1 && key in keys(dict.below_threshold_words)
        cnt_prev = dict.below_threshold_words[key]
        # calculate new count for below threshold word
        cnt = typemax(T) - cnt_prev > cnt ? cnt_prev + cnt : typemax(T)
        # has reached threshold - remove from below threshold
        # collection (it will be added to correct words below)
        if cnt >= dict.count_threshold
            delete!(dict.below_threshold_words, key)
        else
            dict.below_threshold_words[key] = cnt
            return false
        end
    elseif key in keys(dict.words_idx)
        # just update count if it's an already added above
        # threshold word
        idx = dict.words_idx[key]
        cnt_prev = dict.words[idx].count
        dict.words[idx] = DictWord(dict.words[idx].word, typemax(T) - cnt_prev > cnt ? cnt_prev + cnt : typemax(T))
        return false
    elseif cnt < dict.count_threshold
        # new below threshold word
        dict.below_threshold_words[key] = cnt
        return false
    end

    # what we have at this point is a new, above threshold word
    push!(dict.words, DictWord(key, cnt))
    idx = K(length(dict.words))
    dict.words_idx[key] = idx

    # edits/suggestions are created only once, no matter how often
    # word occurs. edits/suggestions are created as soon as the
    # word occurs in the corpus, even if the same term existed
    # before in the dictionary as an edit from another word
    if length(key) > dict.max_length
        dict.max_length = length(key)
    end

    # create deletes
    edits = edits_prefix(key, dict.max_dictionary_edit_distance, dict.prefix_length)
    for delete in edits
        push!(get!(dict.deletes, delete, []), idx)
    end

    return true
end


text:
struct CuDeviceArray{T,N} <: AbstractArray{T,N}
    shape::NTuple{N,Int}
    ptr::Ptr{T}
end

@inline function Base.setindex!(A::CuDeviceArray{T}, x, index::Integer) where {T}
    @boundscheck checkbounds(A, index)
    Base.unsafe_store!(A.ptr, x, index)
end

@inline function Base.getindex(A::CuDeviceArray{T}, index::Integer) where {T}
    @boundscheck checkbounds(A, index)
    Base.unsafe_load(A.ptr, index)::T
end

Base.size(g::CuDeviceArray) = g.shape

Base.checkbounds(::CuDeviceArray, I...) = nothing

@inline blockIdx_x() = Int((ccall("llvm.nvvm.read.ptx.sreg.ctaid.x", llvmcall, UInt32, ()))+UInt32(1))
@inline blockIdx_y() = Int((ccall("llvm.nvvm.read.ptx.sreg.ctaid.y", llvmcall, UInt32, ()))+UInt32(1))
@inline blockIdx_z() = Int((ccall("llvm.nvvm.read.ptx.sreg.ctaid.z", llvmcall, UInt32, ()))+UInt32(1))

@inline blockDim_x() = Int((ccall("llvm.nvvm.read.ptx.sreg.ntid.x", llvmcall, UInt32, ()))+UInt32(1))
@inline blockDim_y() = Int((ccall("llvm.nvvm.read.ptx.sreg.ntid.y", llvmcall, UInt32, ()))+UInt32(1))
@inline blockDim_z() = Int((ccall("llvm.nvvm.read.ptx.sreg.ntid.z", llvmcall, UInt32, ()))+UInt32(1))

@inline threadIdx_x() = Int((ccall("llvm.nvvm.read.ptx.sreg.tid.x", llvmcall, UInt32, ()))+UInt32(1))
@inline threadIdx_y() = Int((ccall("llvm.nvvm.read.ptx.sreg.tid.y", llvmcall, UInt32, ()))+UInt32(1))
@inline threadIdx_z() = Int((ccall("llvm.nvvm.read.ptx.sreg.tid.z", llvmcall, UInt32, ()))+UInt32(1))

@inline blockIdx() = (blockIdx_x(), blockIdx_y(), blockIdx_z())
@inline blockDim() = (blockDim_x(), blockDim_y(), blockDim_z())
@inline threadIdx() = (x=threadIdx_x(), y=threadIdx_y(), z=threadIdx_z())

function kernel_vadd(a, b, c)
    i = threadIdx().x
    c[i] = a[i] + b[i]

    return nothing
end

code_llvm(stdout, kernel_vadd, Tuple{CuDeviceArray{Float32,1},
                                     CuDeviceArray{Float32,1},
                                     CuDeviceArray{Float32,1}})

using TinyCassette
struct Ctx end

code_llvm(stdout, TinyCassette.execute,
          Tuple{Ctx,
                typeof(kernel_vadd),
                CuDeviceArray{Float32,1},
                CuDeviceArray{Float32,1},
                CuDeviceArray{Float32,1}})


text:
using ClimateModels, MITgcmTools, OceanStateEstimation
using Pluto, PlutoUI, PlutoSliderServer, Downloads
import Plots
import CairoMakie

#MC=MITgcm_config(configuration="global_with_exf")
#setup(MC)
# build(MC,"--allow-skip")
# launch(MC)

tmp=ModelConfig(model=ClimateModels.RandomWalker)
setup(tmp)
launch(tmp)

notebook_url="https://raw.githubusercontent.com/gaelforget/OceanStateEstimation.jl/master/examples/ECCO/ECCO_standard_plots.jl"
path_to_notebook = Downloads.download(notebook_url)
include(path_to_notebook)

max_stars_repo_name: gaelforget/Notebooks

max_issues_repo_name: gaelforget/Notebooks

max_forks_repo_name: gaelforget/Notebooks


text: largez04() = largez(4)


text:
using Documenter, ModelingToolkit

makedocs(
    sitename="ModelingToolkit.jl",
    authors="Chris Rackauckas",
    modules=[ModelingToolkit],
    clean=true,doctest=false,
    format=Documenter.HTML(# analytics = "UA-90474609-3",
                             assets=["assets/favicon.ico"],
                             canonical="https://mtk.sciml.ai/stable/"),
    pages=[
        "Home" => "index.md",
        "Symbolic Modeling Tutorials" => Any[
            "tutorials/ode_modeling.md",
            "tutorials/spring_mass.md",
            "tutorials/acausal_components.md",
            "tutorials/higher_order.md",
            "tutorials/tearing_parallelism.md",
            "tutorials/nonlinear.md",
            "tutorials/optimization.md",
            "tutorials/stochastic_diffeq.md",
            "tutorials/nonlinear_optimal_control.md",
            "tutorials/parameter_identifiability.md"
        ],
        "ModelingToolkitize Tutorials" => Any[
            "mtkitize_tutorials/modelingtoolkitize.md",
            "mtkitize_tutorials/modelingtoolkitize_index_reduction.md",
            # "mtkitize_tutorials/sparse_jacobians",
        ],
        "Basics" => Any[
            "basics/AbstractSystem.md",
            "basics/ContextualVariables.md",
            "basics/Composition.md",
            "basics/Validation.md",
            "basics/DependencyGraphs.md",
            "basics/FAQ.md"
        ],
        "System Types" => Any[
            "systems/ODESystem.md",
            "systems/SDESystem.md",
            "systems/JumpSystem.md",
            "systems/NonlinearSystem.md",
            "systems/OptimizationSystem.md",
            "systems/ControlSystem.md",
            "systems/PDESystem.md",
        ],
        "comparison.md",
        "internals.md",
    ]
)

deploydocs(
   repo="github.com/SciML/ModelingToolkit.jl.git";
   push_preview=true
)

max_stars_repo_name: paulrschrater/ModelingToolsJulia

max_issues_repo_name: paulrschrater/ModelingToolsJulia

max_forks_repo_name: paulrschrater/ModelingToolsJulia


text:
"""
$(SIGNATURES)

Generalized surrogate function of `Symbol`
"""
symbolize(x::AbstractString) = Symbol(x)
symbolize(x::AA{T}) where {T<:AbstractString} = [symbolize(y) for y in x]

"""
$(SIGNATURES)

Transform to static array
"""
function static_array(x::AV)
    y = MVector{length(x)}(collect(x))

    if x isa OffsetArray
        idx0 = eachindex(x) |> first
        idx1 = eachindex(x) |> last

        y = OffsetArray(y, idx0:idx1)
    end

    return y
end

function static_array(x::AM)
    y = MMatrix{size(x, 1),size(x, 2)}(collect(x))

    if x isa OffsetArray
        idx0 = axes(x, 1) |> first
        idx1 = axes(x, 1) |> last
        idy0 = axes(x, 2) |> first
        idy1 = axes(x, 2) |> last

        y = OffsetArray(y, idx0:idx1, idy0:idy1)
    end

    return y
end

function static_array(x::AA{T,3}) where {T}
    y = MArray{Tuple{size(x, 1),size(x, 2),size(x, 3)}}(collect(x))

    if x isa OffsetArray
        idx0 = axes(x, 1) |> first
        idx1 = axes(x, 1) |> last
        idy0 = axes(x, 2) |> first
        idy1 = axes(x, 2) |> last
        idz0 = axes(x, 3) |> first
        idz1 = axes(x, 3) |> last

        y = OffsetArray(y, idx0:idx1, idy0:idy1, idz0:idz1)
    end

    return y
end

function static_array(x::AA{T,4}) where {T}
    y = MArray{Tuple{size(x, 1),size(x, 2),size(x, 3),size(x, 4)}}(collect(x))

    if x isa OffsetArray
        ida0 = axes(x, 1) |> first
        ida1 = axes(x, 1) |> last
        idb0 = axes(x, 2) |> first
        idb1 = axes(x, 2) |> last
        idc0 = axes(x, 3) |> first
        idc1 = axes(x, 3) |> last
        idd0 = axes(x, 4) |> first
        idd1 = axes(x, 4) |> last

        y = OffsetArray(y, ida0:ida1, idb0:idb1, idc0:idc1, idd0:idd1)
    end

    return y
end

"""
$(SIGNATURES)

Transform to dynamic array
"""
dynamic_array(x) = x
dynamic_array(x::StaticArray) = Array(x)


text:
import Flux
import Flux.Tracker

struct TrackerGradAlg <: GradAlg end
const TrackerGrad = TrackerGradAlg()

function back!(rv, ω, ::TrackerGradAlg)
  l = rv(ω)
  Flux.back!(l)
end

# # FIXME: Remove this
# function gradient(::TrackerGradAlg, U, sω::SimpleΩ{I, V}) where {I, V <: AbstractArray}
#   l = U(sω)
#   Flux.back!(l)
# end
# # The New #

struct GradView{G, O}
  grads::G
  vals::O
end
Base.getindex(gv::GradView, i) = value(gv.grads[gv.vals[i]])

function grad(rv, ω::Ω, ::TrackerGradAlg)
  v = collect(values(ω))
  grad_ = Tracker.gradient(() -> rv(ω), Flux.params(v...))
  GradView(grad_, v)
end

function gradarray(rv, ω::Ω, ::TrackerGradAlg)
  vs = values(ω)
  grads_ = grad(rv, ω, vs, TrackerGrad)
  map(v -> grads_[v], vs)
end

lineargradient(rv, ω, ::TrackerGradAlg) = (back!(rv, ω, TrackerGrad); linearize(ω))

function value(x::Union{Tracker.Tracked,
                        Tracker.TrackedReal,
                        Tracker.TrackedArray,
                        Tracker.TrackedTuple})
  Tracker.data(x)
end


text:
global reps = 10

@testset "Construction and Conversions" begin
    @test DNACodon(DNA_A, DNA_G, DNA_T) === DNAMer("AGT")
    @test RNACodon(RNA_A, RNA_G, RNA_U) === RNAMer("AGU")

    # Check that kmers in strings survive round trip conversion:
    #   String → Kmer → String
    function check_string_construction(::Type{T}, seq::AbstractString) where {T<:AbstractMer}
        return String(T(seq)) == uppercase(seq)
    end

    # Check that DNAMers can be constructed from a LongDNASeq
    #   LongDNASeq → Kmer → LongDNASeq
    function check_dnasequence_construction(::Type{T}, seq::LongDNASeq) where {T<:AbstractMer}
        return LongDNASeq(T(seq)) == seq
    end

    # Check that RNAMers can be constructed from a LongRNASeq
    #   LongRNASeq → Kmer → LongRNASeq
    function check_rnasequence_construction(::Type{T}, seq::LongRNASeq) where {T<:AbstractMer}
        return LongRNASeq(T(seq)) == seq
    end

    # Check that kmers can be constructed from a BioSequence
    #   BioSequence → Kmer → BioSequence
    function check_biosequence_construction(::Type{T}, seq::LongSequence) where {T<:AbstractMer}
        return LongSequence(T(seq)) == seq
    end

    # Check that kmers can be constructed from an array of nucleotides
    #   Vector{NucleicAcid} → Kmer → Vector{NucleicAcid}
    function check_nucarray_kmer(::Type{M}, seq::Vector{T}) where {T<:NucleicAcid,M<:AbstractMer}
        return String([convert(Char, c) for c in seq]) == String(M(seq))
    end

    # Check that kmers in strings survive round trip conversion:
    #   String → BioSequence → Kmer → BioSequence → String
    function check_roundabout_construction(::Type{T}, A2, seq::AbstractString) where {T<:AbstractMer}
        return String(LongSequence{A2}(T(LongSequence{A2}(seq)))) == uppercase(seq)
    end

    function check_uint_conversion(::Type{T}) where {T<:AbstractMer}
        U = BioSequences.encoded_data_type(T)
        uint = rand(typemin(U):U(one(U) << 2BioSequences.ksize(T) - 1))
        return convert(U, T(uint)) === uint
    end

    @testset "Skipmer conversion" begin
        for len in [1, 16, 32, 64]

            if len <= 32
                # UInt64 conversions
                @test all(Bool[check_uint_conversion(DNAMer{len}) for _ in 1:reps])
                @test all(Bool[check_uint_conversion(RNAMer{len}) for _ in 1:reps])
                # String construction
                @test all(Bool[check_string_construction(DNAMer{len}, random_dna_kmer(len)) for _ in 1:reps])
                @test all(Bool[check_string_construction(RNAMer{len}, random_rna_kmer(len)) for _ in 1:reps])
                # DNA/LongRNASeq Constructions
                @test all(Bool[check_dnasequence_construction(DNAMer{len}, LongDNASeq(random_dna_kmer(len))) for _ in 1:reps])
                @test all(Bool[check_rnasequence_construction(RNAMer{len}, LongRNASeq(random_rna_kmer(len))) for _ in 1:reps])
                # BioSequence Construction
                @test all(Bool[check_biosequence_construction(DNAMer{len}, LongDNASeq(random_dna_kmer(len))) for _ in 1:reps])
                @test all(Bool[check_biosequence_construction(RNAMer{len}, LongRNASeq(random_rna_kmer(len))) for _ in 1:reps])
                # Construction from nucleotide arrays
                @test all(Bool[check_nucarray_kmer(DNAMer{len}, random_dna_kmer_nucleotides(len)) for _ in 1:reps])
                @test all(Bool[check_nucarray_kmer(RNAMer{len}, random_rna_kmer_nucleotides(len)) for _ in 1:reps])
                # Roundabout conversions
                @test all(Bool[check_roundabout_construction(DNAMer{len}, DNAAlphabet{2}, random_dna_kmer(len)) for _ in 1:reps])
                @test all(Bool[check_roundabout_construction(DNAMer{len}, DNAAlphabet{4}, random_dna_kmer(len)) for _ in 1:reps])
                @test all(Bool[check_roundabout_construction(RNAMer{len}, RNAAlphabet{2}, random_rna_kmer(len)) for _ in 1:reps])
                @test all(Bool[check_roundabout_construction(RNAMer{len}, RNAAlphabet{4}, random_rna_kmer(len)) for _ in 1:reps])
            end

            # UInt64 conversions
            @test all(Bool[check_uint_conversion(BigDNAMer{len}) for _ in 1:reps])
            @test all(Bool[check_uint_conversion(BigDNAMer{len}) for _ in 1:reps])
            # String construction
            @test all(Bool[check_string_construction(BigDNAMer{len}, random_dna_kmer(len)) for _ in 1:reps])
            @test all(Bool[check_string_construction(BigRNAMer{len}, random_rna_kmer(len)) for _ in 1:reps])
            # DNA/LongRNASeq Constructions
            @test all(Bool[check_dnasequence_construction(BigDNAMer{len}, LongDNASeq(random_dna_kmer(len))) for _ in 1:reps])
            @test all(Bool[check_rnasequence_construction(BigRNAMer{len}, LongRNASeq(random_rna_kmer(len))) for _ in 1:reps])
            # BioSequence Construction
            @test all(Bool[check_biosequence_construction(BigDNAMer{len}, LongDNASeq(random_dna_kmer(len))) for _ in 1:reps])
            @test all(Bool[check_biosequence_construction(BigRNAMer{len}, LongRNASeq(random_rna_kmer(len))) for _ in 1:reps])
            # Construction from nucleotide arrays
            @test all(Bool[check_biosequence_construction(BigDNAMer{len}, LongDNASeq(random_dna_kmer(len))) for _ in 1:reps])
            @test all(Bool[check_biosequence_construction(BigRNAMer{len}, LongRNASeq(random_rna_kmer(len))) for _ in 1:reps])
            # Roundabout conversions
            @test all(Bool[check_roundabout_construction(BigDNAMer{len}, DNAAlphabet{2}, random_dna_kmer(len)) for _ in 1:reps])
            @test all(Bool[check_roundabout_construction(BigDNAMer{len}, DNAAlphabet{4}, random_dna_kmer(len)) for _ in 1:reps])
            @test all(Bool[check_roundabout_construction(BigRNAMer{len}, RNAAlphabet{2}, random_rna_kmer(len)) for _ in 1:reps])
            @test all(Bool[check_roundabout_construction(BigRNAMer{len}, RNAAlphabet{4}, random_rna_kmer(len)) for _ in 1:reps])
        end
    end

    @test_throws MethodError Mer() # can't construct 0-mer using `Kmer()`
    @test_throws MethodError BigMer() # can't construct 0-mer using `BigKmer()`
    @test_throws ArgumentError DNAMer(dna"") # 0-mers not allowed
    @test_throws ArgumentError BigDNAMer(dna"") # 0-mers not allowed
    @test_throws ArgumentError DNAMer{0}(UInt64(0)) # 0-mers not allowed
    @test_throws ArgumentError RNAMer{0}(UInt64(0)) # 0-mers not allowed
    @test_throws ArgumentError RNAMer(RNA_A, RNA_C, RNA_G, RNA_N, RNA_U) # no Ns in kmers
    @test_throws ArgumentError DNAMer(DNA_A, DNA_C, DNA_G, DNA_N, DNA_T) # no Ns in kmers
    @test_throws ArgumentError RNAMer(rna"ACGNU")# no Ns in kmers
    @test_throws ArgumentError DNAMer(dna"ACGNT") # no Ns in kmers
    @test_throws MethodError DNAMer(RNA_A, DNA_A) # no mixing of RNA and DNA
    @test_throws ArgumentError RNAMer(random_rna(33)) # no kmer larger than 32nt
    @test_throws ArgumentError DNAMer(random_dna(33)) # no kmer larger than 32nt
    @test_throws ArgumentError RNAMer(
                      RNA_A, RNA_C, RNA_G, RNA_U, # no kmer larger than 32nt
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U)
    @test_throws ArgumentError DNAMer(
                      DNA_A, DNA_C, DNA_G, DNA_T, # no kmer larger than 32nt
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T)

    @test_throws ArgumentError BigDNAMer{0}(UInt128(0)) # 0-mers not allowed
    @test_throws ArgumentError BigRNAMer{0}(UInt128(0)) # 0-mers not allowed
    @test_throws ArgumentError BigRNAMer(RNA_A, RNA_C, RNA_G, RNA_N, RNA_U) # no Ns in kmers
    @test_throws ArgumentError BigDNAMer(DNA_A, DNA_C, DNA_G, DNA_N, DNA_T) # no Ns in kmers
    @test_throws ArgumentError BigRNAMer(rna"ACGNU")# no Ns in kmers
    @test_throws ArgumentError BigDNAMer(dna"ACGNT") # no Ns in kmers
    @test_throws MethodError BigDNAMer(RNA_A, DNA_A) # no mixing of RNA and DNA
    @test_throws ArgumentError BigRNAMer(random_rna(65)) # no kmer larger than 64nt
    @test_throws ArgumentError BigDNAMer(random_dna(65)) # no kmer larger than 64nt
    @test_throws ArgumentError BigRNAMer(
                      RNA_A, RNA_C, RNA_G, RNA_U, # no kmer larger than 64nt
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U,
                      RNA_A, RNA_C, RNA_G, RNA_U)
    @test_throws ArgumentError BigDNAMer(
                      DNA_A, DNA_C, DNA_G, DNA_T, # no kmer larger than 64nt
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T,
                      DNA_A, DNA_C, DNA_G, DNA_T)

    @testset "From strings" begin
        @test DNAMer("ACTG") == DNAMer(LongDNASeq("ACTG"))
        @test RNAMer("ACUG") == RNAMer(LongRNASeq("ACUG"))
        @test BigDNAMer("ACTG") == BigDNAMer(LongDNASeq("ACTG"))
        @test BigRNAMer("ACUG") == BigRNAMer(LongRNASeq("ACUG"))

        # N is not allowed in Kmers
        @test_throws Exception DNAMmer("ACGTNACGT")
        @test_throws Exception RNAMer("ACGUNACGU")
        @test_throws Exception BigDNAMmer("ACGTNACGT")
        @test_throws Exception BigRNAMer("ACGUNACGU")

        # Test string literals
        @test mer"ACTG"dna == DNAMer(LongDNASeq("ACTG"))
        @test isa(mer"ACGT"dna, DNAMer{4})
        @test_throws LoadError eval(:(mer"ACGN"dna))
        @test_throws LoadError eval(:(mer"ACG-"dna))

        @test bigmer"ACTG"dna == BigDNAMer(LongDNASeq("ACTG"))
        @test isa(bigmer"ACGT"dna, BigDNAMer{4})
        @test_throws LoadError eval(:(bigmer"ACGN"dna))
        @test_throws LoadError eval(:(bigmer"ACG-"dna))
    end

    @testset "Capacity" begin
        @test BioSequences.capacity(DNAMer(random_dna_kmer(10))) == 32
        @test BioSequences.capacity(RNAMer(random_rna_kmer(10))) == 32
        @test BioSequences.capacity(DNAMer(random_dna_kmer(32))) == 32
        @test BioSequences.capacity(RNAMer(random_rna_kmer(32))) == 32

        @test BioSequences.capacity(BigDNAMer(random_dna_kmer(10))) == 64
        @test BioSequences.capacity(BigRNAMer(random_rna_kmer(10))) == 64
        @test BioSequences.capacity(BigDNAMer(random_dna_kmer(32))) == 64
        @test BioSequences.capacity(BigRNAMer(random_rna_kmer(32))) == 64
    end

    @testset "N unused" begin
        @test BioSequences.n_unused(DNAMer(random_dna_kmer(10))) == 22
        @test BioSequences.n_unused(RNAMer(random_rna_kmer(10))) == 22
        @test BioSequences.n_unused(DNAMer(random_dna_kmer(32))) == 0
        @test BioSequences.n_unused(RNAMer(random_rna_kmer(32))) == 0

        @test BioSequences.n_unused(BigDNAMer(random_dna_kmer(40))) == 24
        @test BioSequences.n_unused(BigRNAMer(random_rna_kmer(40))) == 24
        @test BioSequences.n_unused(BigDNAMer(random_dna_kmer(64))) == 0
        @test BioSequences.n_unused(BigRNAMer(random_rna_kmer(64))) == 0
    end
end


text:
using Test
using Ipaper

# println(dirname(@__FILE__))
# println(pwd())

# cd(dirname(@__FILE__)) do

include("test-Ipaper.jl")
include("test-missing.jl")
include("test-statistics.jl")
include("test-Pipe.jl")
include("test-dt.jl")
include("test-dt_pipe.jl")

# include("test-smooth_whit.jl")
# include("test-smooth_SG.jl")
# include("test_wTSM.jl")
# include("test_whittaker.jl")
# include("test-lambda_init.jl")
# end


text:
using ComputationalBiologyCourse, Test

@test 1==1

@testset "Builders" begin
    @test 1==1
end


text:
module ArrayInitialization

using DocStringExtensions
using Unitful
using AxisArrays
using ArgCheck
using ..ConstantArrays: ConstantArray
using ..UnitfulHartree
using ..Traits: components, ColinearSpin, SpinDegenerate, SpinCategory, ColinearSpin,
                ColinearSpinFirst, ColinearSpinLast, is_spin_polarized, concretize_type,
                SpinAware
using ..Dispatch
export wrap

macro lintpragma(s) end
@lintpragma("Ignore unused args")
@lintpragma("Ignore unused T")

const DD = Dispatch.Dimensions

_size(T::Type{<: DD.Scalars.All}, ::SpinDegenerate,
      ::SpinDegenerate, object::DD.AxisArrays.All) = size(object)
_size(T::Type{<: DD.Scalars.All}, ::ColinearSpin,
      ::SpinDegenerate, object::DD.AxisArrays.All) =
    (size(object)..., length(components(object)))

add_spin_axis(::SpinDegenerate, tup::Tuple, ::Any) = tup
add_spin_axis(::ColinearSpinLast, tup::Tuple, axis::Any) = tup..., axis
add_spin_axis(::ColinearSpin{N}, tup::Tuple, axis::Any) where N = begin
    @lintpragma("Ignore use of undeclared variable N")
    @lintpragma("Ignore use of undeclared variable axis")
    @lintpragma("Ignore use of undeclared variable tup")
    top, down = Base.IteratorsMD.split(tup, Val{N})

    Base.front(top)..., axis, Base.last(top), down...
end

replace_spin_axis(T::Type{<: DD.Scalars.All}, ::SpinDegenerate,
                  sizes::Tuple, ax::Tuple) = begin
    @assert length(sizes.parameters) == length(ax.parameters)
    @assert findfirst(is_spin_polarized, ax.parameters) == 0
    sizes, ax
end

"""
Position for the spin axis for a given spin category

$(SIGNATURES)

`n` is the number of dimensions in the array, and `s` is the location of the current spin
axis, or 0 if there are none.
"""
spin_axis_position(::Type{ColinearSpin{T}}, ::Integer, ::Integer) where T = T
spin_axis_position(::Type{ColinearSpinLast}, n::Integer, s::Integer) =
    !(1 ≤ s ≤ n) ? n + 1: n
spin_axis_position(S::SpinCategory, n::Integer, s::Integer) =
    spin_axis_position(typeof(S), n, s)
spin_axis_position(array::DD.AxisArrays.All) = findfirst(is_spin_polarized, axes(array))
spin_axis_position(T::Type{DD.AxisArrays.All}) =
    findfirst(is_spin_polarized, T.parameters[end].parameters)

function _move_and_re(name::Symbol, expr::Expr, orig::Integer, final::Integer, n::Integer)
    result = Expr(:tuple)
    for i in 1:n
        i > orig && push!(result.args, :($name[$i]))
        i == final && push!(result.args, expr)
        i < orig && push!(result.args, :($name[$i]))
    end
    final == n + 1 && push!(result.args, expr)
    result
end

@generated replace_spin_axis(T::Type{<: DD.Scalars.All},
                             S::ColinearSpin, sizes::Tuple, ax::Tuple) = begin
    @lintpragma("Ignore use of undeclared variable comps")
    @assert length(sizes.parameters) == length(ax.parameters)
    current = findfirst(is_spin_polarized, ax.parameters)
    final = spin_axis_position(S, length(sizes.parameters), current)
    @assert final ≠ 0
    d = _move_and_re(:sizes, :(length(comps)), current, final, length(sizes.parameters))
    a = _move_and_re(:ax, :(Axis{:spin}(comps)), current, final, length(sizes.parameters))
    quote
        comps = components(T, S)
        ($d, $a)
    end
end

@generated remove_spin_axis(sizes::Tuple, ax::Tuple) = begin
    @lintpragma("Ignore use of undeclared variable comps")
    @assert length(sizes.parameters) == length(ax.parameters)
    current = findfirst(is_spin_polarized, ax.parameters)
    d = _move_and_re(:sizes, :(length(comps)), current, 0, length(sizes.parameters))
    a = _move_and_re(:ax, :(Axis{:spin}(comps)), current, 0, length(sizes.parameters))
    :(($d, $a))
end

for extension in [:zeros, :ones, :rand]
    @eval begin
        """
        Creates an array for the given DFT quantity

        The spin axis is automatically added if required. Neither `dims` nor `ax` should
        include the spin dimension.
        """
        Base.$extension(T::Type{<:DD.Scalars.All}, C::SpinCategory,
                        dims::Tuple, ax::Tuple) = begin
            @argcheck length(ax) <= (length(dims) + 1) "Too many axes"
            comps = components(T, C)
            data = $extension(T.parameters[1], add_spin_axis(C, dims, length(comps)))
            # we can use a dummy array here since the underlying array (and indices) will be
            # the standard one
            defaults = AxisArrays.default_axes(ConstantArray(0, dims), ax)
            AxisArray(reinterpret(T, data), add_spin_axis(C, defaults, Axis{:spin}(comps)))
        end
        Base.$extension(T::Type{<:DD.Scalars.All}, C::SpinCategory, dims::Tuple) = begin
            comps = components(T, C)
            data = $extension(T, add_spin_axis(C, dims, length(comps)))
            defaults = AxisArrays.default_axes(ConstantArray(0, dims))
            AxisArray(data, add_spin_axis(C, defaults, Axis{:spin}(comps)))
        end
        Base.$extension(T::Type{<:DD.Scalars.All}, C::SpinCategory,
                        ax::Tuple{Axis, Vararg{Axis}}) = begin
            @argcheck(SpinCategory(ax) isa SpinDegenerate,
                      "Do not include Spin axis explicitly when specifying spin category")
            $extension(T, C, length.(ax), ax)
        end
        Base.$extension(T::Type{<:DD.Scalars.All}, C::SpinCategory, ax::Axis...) =
            $extension(T, C, ax)
        Base.$extension(T::Type{<:DD.Scalars.All}, ax::Axis...) = $extension(T, ax)
        Base.$extension(T::Type{<:DD.Scalars.All}, ax::Tuple{Axis, Vararg{Axis}}) = begin
            const C = SpinCategory(ax)
            if C isa ColinearSpin
                index = findfirst(a -> a isa Axis{:spin}, ax)
                @argcheck(length(components(T, C)) == length(ax[index]),
                          "Incorrect spin components for given input type")
            end
            AxisArray(reinterpret(T, $extension(T.parameters[1], length.(ax))), ax)
        end

        """
        Creates an array for the given DFT quantity

        The spin components, if any, are added as the dimension specified by the spin trait.
        Note that the spin trait is the only way polarization should enter the input
        arguments, e.g. the spin-axis should not be given explicitly, neither in the axes
        nor in the dimensions.
        """
        Base.$extension(T::Type{<:DD.Scalars.All}, P::SpinCategory,
                        args::Vararg{Union{Integer, Axis}}; kwargs...) = begin
            @lintpragma("Ignore use of undeclared variable x")
            dims = (
                (x for x in args if typeof(x) <: Integer)...,
                (length(v) for (_, v) in kwargs)...,
            )
            ax = (
                (x for x in args if typeof(x) <: Axis)...,
                (Axis{k}(v) for (k, v) in kwargs)...,
            )
            $extension(T, P, dims, ax)
        end
        Base.$extension(T::Type{<:DD.Scalars.All}, polarized::Bool,
                        args::Vararg{Union{Integer, Axis}}; kwargs...) =
            $extension(T, polarized ? ColinearSpin(): SpinDegenerate(), args...; kwargs...)
    end
end

for extension in [:zeros, :ones, :similar]
    private = Symbol(:_, extension)
    @eval begin
        $private(T::Type{<:DD.Scalars.All}, ::SpinDegenerate, ::SpinDegenerate,
            array::DD.AxisArrays.All) = $extension(array, concretize_type(T, array))

        $private(T::Type{<:DD.Scalars.All}, C::ColinearSpin,
                 ::SpinDegenerate, array::DD.AxisArrays.All) = begin
            comps = components(T, C)
            dims = add_spin_axis(C, size(array), length(comps))
            data = $extension(array.data, concretize_type(T, array), dims)
            defaults = AxisArrays.default_axes(ConstantArray(0, size(array)), axes(array))
            AxisArray(data, add_spin_axis(C, defaults, Axis{:spin}(comps)))
        end

        $private(T::Type{<:DD.Scalars.All}, C::ColinearSpin,
                 ::ColinearSpin, array::DD.AxisArrays.All) = begin
            dims = replace_spin_axis(T, C, size(array), axes(array))
            AxisArray($extension(array.data, concretize_type(T, array), dims[1]), dims[2])
        end

        $private(T::Type{<:DD.Scalars.All}, ::SpinDegenerate,
                    ::ColinearSpin, array::DD.AxisArrays.All) = begin
            dims = remove_spin_axis(size(array), axes(array))
            AxisArray($extension(array.data, concretize_type(T, array), dims[1]), dims[2])
        end

        Base.$extension(array::DD.AxisArrays.All, T::Type{<:DD.Scalars.All}, ::SpinAware) =
            $private(T, SpinCategory(array), SpinCategory(array), array)
        Base.$extension(array::DD.AxisArrays.All, ::SpinAware) =
            $private(eltype(array), SpinCategory(array), SpinCategory(array), array)
        Base.$extension(array::DD.AxisArrays.All, C::SpinCategory) =
            $private(eltype(array), C, SpinCategory(array), array)
        Base.$extension(array::DD.AxisArrays.All, T::Type{<:DD.Scalars.All},
                        wanted::SpinCategory) =
            $private(T, wanted, SpinCategory(array), array)

        Base.$extension(array::DD.DenseArrays.All, T::Type{<:DD.Scalars.All},
                        ::SpinDegenerate) = similar(array, T)
        Base.$extension(array::DD.DenseArrays.All, T::Type{<:DD.Scalars.All},
                        C::ColinearSpin{N}) where N = begin
            @lintpragma("Ignore use of undeclared variable array")
            @lintpragma("Ignore use of undeclared variable C")
            @lintpragma("Ignore use of undeclared variable T")
            @argcheck 1 ≤ N ≤ ndims(array)
            left, right = Base.IteratorsMD.split(size(array), Val{N})
            n = Base.last(left)
            @argcheck n == length(components(eltype(array), C))

            n′ = length(components(T, C))
            similar(array, T, Base.front(left)..., n′, right...)
        end
        Base.$extension(array::DD.DenseArrays.All, T::Type{<:DD.Scalars.All},
                        C::ColinearSpinLast) = begin
            @argcheck size(array, ndims(array)) == length(components(eltype(array), C))
            n′ = length(components(T, C))
            similar(array, T, Base.front(size(array))..., n′)
        end
    end
end

Base.reinterpret(T::Type{<: DD.Scalars.All}, ::SpinDegenerate, array::DenseArray) =
    AxisArray(reinterpret(concretize_type(T, array), array), axes(array))
Base.reinterpret(T::Type{<: DD.Scalars.All}, ::ColinearSpinFirst, array::DenseArray) =
    AxisArray(reinterpret(concretize_type(T, array), array),
              Axis{:spin}(components(T, ColinearSpinFirst())), Base.tail(axes(array))...)
Base.reinterpret(T::Type{<: DD.Scalars.All}, C::ColinearSpin, array::DenseArray) =
    AxisArray(reinterpret(concretize_type(T, array), array),
              Base.front(axes(array))..., Axis{:spin}(components(T, C)))

"""
    wrap([T::Type{<: Dispatch.Scalars.All}], [P::SpinCategory], array)

Wraps an existing array as a DFT array. The default spin-category is `SpinDegenerate`.
If the type and units are explicitly specified (first argument)), then the array must be
dimensionless. If they are not, then the array must have an element type `T` such that
`T <: Dispatch.Scalars.All`.
"""
wrap(T::Type{<: DD.Scalars.All}, P::SpinCategory, array::DenseArray) = begin
    @argcheck dimension(eltype(array)) == NoDims "Input array should be dimensionless"
    reinterpret(T, P, array)
end
wrap(::SpinDegenerate, array::DenseArray{<: DD.Scalars.All}) = AxisArray(array)
wrap(P::ColinearSpinFirst, array::DenseArray{<: DD.Scalars.All}) = begin
    c = components(eltype(array), P)
    @argcheck length(c) == size(array, 1) "Size of spin-axis does not match expectations"
    AxisArray(array, Axis{:spin}(c))
end
wrap(P::ColinearSpin, array::DenseArray{<: DD.Scalars.All}) = begin
    c = components(eltype(array), P)
    @argcheck(length(c) == size(array, ndims(array)),
              "Size of spin-axis does not match expectations")
    AxisArray(array, Base.front(axes(array))..., Axis{:spin}(c))
end
wrap(array::DenseArray{<: DD.Scalars.All}) = wrap(SpinDegenerate(), array)
wrap(T::Type{<: DD.Scalars.All}, array::DenseArray) =
    wrap(concretize_type(T, array), SpinDegenerate(), array)

""" Helper functions to convert between arrays with different units and memory layout """
_convert(T::Type{<:DD.Scalars.All}, C::ColinearSpin,
         C′::ColinearSpin, array::DD.AxisArrays.All) = begin
    T′ = concretize_type(T, array)
    (T′ == T && C == C′) ? array: copy!(similar(array, T, C), array)
end

"""
Converts axis to the requisite spin-axis location

$(SIGNATURES)

If the spin-axis does not move, then a reference to the original array is returned.
Otherwise, a new array is returned.
"""
Base.convert(::SpinAware, array::DD.AxisArrays.All) = array
Base.convert(T::Type{<: DD.Scalars.All}, ::SpinAware, array::DD.AxisArrays.All) =
    _convert(concretize_type(T, array), SpinCategory(array), SpinCategory(array), array)
Base.convert(C::ColinearSpin, array::DD.AxisArrays.All) =
    _convert(eltype(array), C, SpinCategory(array), array)
Base.convert(T::Type{<: DD.Scalars.All}, C::SpinCategory, array::DD.AxisArrays.All) =
    _convert(T, C, SpinCategory(array), array)
Base.convert(T::Type{<: DD.Scalars.All}, array::DD.AxisArrays.All) =
    _convert(T, SpinCategory(array), SpinCategory(array), array)

Unitful.uconvert(u::Unitful.Units, array::AxisArray) = begin
    data = similar(array.data, typeof(uconvert(u, oneunit(eltype(array)))))
    data .= array.data
    AxisArray(data, axes(array))
end

end

max_stars_repo_name: mdavezac/DFTShims

max_issues_repo_name: mdavezac/DFTShims

max_forks_repo_name: mdavezac/DFTShims


text:
using Distributions
using Random

function sample_infection_time(;omicron=false)
    """
    Sample infection times for num individuals based on the generation
    interval distribution, Gamma(shape_gen, scale_gen).
    """

    (shape_gen, scale_gen) = (2.75, 1.00)
    (shape_gen_omicron, scale_gen_omicron) = (1.58, 1.32)

    shape = (1-omicron)*shape_gen + omicron*shape_gen_omicron
    scale = (1-omicron)*scale_gen + omicron*scale_gen_omicron

    infection_time = ceil(Int, rand(Gamma(shape, scale)))

    return infection_time
end

function sample_onset_time(;omicron=false)
    """
    Sample incubation times for num individuals based on incubation period
    distribution, Gamma(shape_inc, scale_inc).
    """

    (shape_inc, scale_inc) = (5.807, 0.948)
    (shape_inc_omicron, scale_inc_omicron) = (3.33, 1.34)

    shape = (1-omicron)*shape_inc + omicron*shape_inc_omicron
    scale = (1-omicron)*scale_inc + omicron*scale_inc_omicron

    onset_time = ceil(Int, rand(Gamma(shape, scale)))

    return onset_time
end

function set_simulation_constants(state)
    """
    Contains the assumptions for simulation parameters. This includes all the dynamical constants:
        - k = heterogeneity parameter
        - p_symp = probability of symptoms
        - γ = relative infectiousness of asymptomatic individuals
        - p_symp_given_detect = probability of symptoms given detection
        - p_asymp_given_detect = probability of being asymptomatic given detection
        - consistency_multiplier = chosen such that sim_cases < consistency_multiplier*actual cases
            results in cases being injected into the simulation. This is used to account
            for superspreading events after periods of low incidence.
    These values are stored in sim_constants which is a dictionary indexed by the parameter name and
    ultimately stored on the stack in a SimulationParameters object.
    """

    # overdispersion parameter
    k = 0.15
    # assumptions surrouding the probability of symptomatic,
    # relative infectiousness γ and the ratio of Reff (α's)
    p_symp = 0.4

    p_detect_given_symp_dict = Dict{String, Float64}(
        "NSW" => 0.7,
        "QLD" => 0.7,
        "SA" => 0.7,
        "TAS" => 0.7,
        "WA" => 0.7,
        "ACT" => 0.7,
        "NT" => 0.7,
        "VIC" => 0.7,
    )

    p_detect_given_asymp_dict = Dict{String, Float64}(
        "NSW" => 0.467,
        "QLD" => 0.467,
        "SA" => 0.467,
        "TAS" => 0.467,
        "WA" => 0.467,
        "ACT" => 0.467,
        "NT" => 0.467,
        "VIC" => 0.467,
    )

    γ = 0.5     # relative infectiousness of asymptomatic
    # solve the system α_s*ps + α_a(1-ps) = 1 with α_a = γ*α_s
    p_detect_given_symp = p_detect_given_symp_dict[state]
    p_detect_given_asymp = p_detect_given_asymp_dict[state]
    # prob of detection
    p_detect = p_symp*p_detect_given_symp + (1-p_symp)*p_detect_given_asymp
    # prob symptomatic given detect
    p_symp_given_detect = p_detect_given_symp*p_symp/p_detect

    α_s = 1/(p_symp + γ*(1-p_symp))
    α_a = γ * α_s
    consistency_multiplier = 5.0

    # prob of detecting an international import
    qi = 0.98
    # prior parametes for the import model
    prior_alpha = 0.5
    prior_beta = 0.2
    # ema smoothing factor
    ϕ = 0.1

    simulation_constants = Constants(
        k,
        p_symp,
        p_detect_given_symp,
        p_detect_given_asymp,
        p_detect,
        p_symp_given_detect,
        γ,
        α_s,
        α_a,
        qi,
        prior_alpha,
        prior_beta,
        ϕ,
        consistency_multiplier,
    )

    # mapping between types
    individual_type_map = IndividualTypeMap(1, 2, 3)

    return (simulation_constants, individual_type_map)
end


text:
using Revise, CUDA, LinearAlgebra, Random

using CUDA.CUBLAS
using CUDA.CUBLAS: band, bandex

Random.seed!(11)
function mwe()

  local m = 20
  local n = 35
  local k = 13

  elty=Float32

  local alpha = rand(elty)
  local beta = rand(elty)

  local A = triu(rand(elty, m, m))
  local B = rand(elty,m,n)
  local C = zeros(elty,m,n)
  local dA = CuArray(A)
  local dB = CuArray(B)
  local dC = CuArray(C)
  local failed=false

  try
    C = alpha*(A\B)
    dC = copy(dB)
    CUBLAS.xt_trsm!('L','U','N','N',alpha,dA,dC)
    CUDA.synchronize()
    # move to host and compare
    h_C = Array(dC)
    @assert C ≈ h_C
  catch err
    @warn "xt_trsm! gpu failed!! error: $err"
    failed=true
  end

  try
    C  = alpha*(A\B)
    h_C = CUBLAS.xt_trsm('L','U','N','N',alpha,Array(dA),Array(dB))
    CUDA.synchronize()
    @assert C ≈ h_C
  catch err
    @warn "xt_trsm cpu failed!! error: $err"
    failed=true
  end

  return failed
end

for i ∈ 1:10^3
  if mwe()
    @info "Failed on iteration $i"
    break
  end
end
#=
try
    # generate matrices
    syrkx_A = rand(elty, n, k)
    syrkx_B = rand(elty, n, k)
    d_syrkx_A = CuArray(syrkx_A)
    d_syrkx_B = CuArray(syrkx_B)
    d_syrkx_C = CUBLAS.xt_syrkx('U','N',d_syrkx_A,d_syrkx_B)
    final_C = syrkx_A*transpose(syrkx_B)
    # move to host and compare
    h_C = Array(d_syrkx_C)
    @assert triu(final_C) ≈ triu(h_C)
catch err
    @warn "xt_syrkx gpu failed!! error: $err"
end=#

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# With your submarine's subterranean subsystems subsisting suboptimally, the
# only way you're getting out of this cave anytime soon is by finding a path
# yourself. Not just a path - the only way to know if you've found the best
# path is to find all of them.
#
# Fortunately, the sensors are still mostly working, and so you build a rough
# map of the remaining caves (your puzzle input). For example:
#
# start-A
# start-b
# A-c
# A-b
# b-d
# A-end
# b-end
#
# This is a list of how all of the caves are connected. You start in the cave
# named start, and your destination is the cave named end. An entry like b-d
# means that cave b is connected to cave d - that is, you can move between
# them.
#
# So, the above cave system looks roughly like this:
#
#     start
#     /   \
# c--A-----b--d
#     \   /
#      end
#
# Your goal is to find the number of distinct paths that start at start, end at
# end, and don't visit small caves more than once. There are two types of
# caves: big caves (written in uppercase, like A) and small caves (written in
# lowercase, like b). It would be a waste of time to visit any small cave more
# than once, but big caves are large enough that it might be worth visiting
# them multiple times. So, all paths you find should visit small caves at most
# once, and can visit big caves any number of times.
#
# Given these rules, there are 10 paths through this example cave system:
#
# start,A,b,A,c,A,end
# start,A,b,A,end
# start,A,b,end
# start,A,c,A,b,A,end
# start,A,c,A,b,end
# start,A,c,A,end
# start,A,end
# start,b,A,c,A,end
# start,b,A,end
# start,b,end
#
# (Each line in the above list corresponds to a single path; the caves visited
# by that path are listed in the order they are visited and separated by
# commas.)
#
# Note that in this cave system, cave d is never visited by any path: to do so,
# cave b would need to be visited twice (once on the way to cave d and a second
# time when returning from cave d), and since cave b is small, this is not
# allowed.
#
# Here is a slightly larger example:
#
# dc-end
# HN-start
# start-kj
# dc-start
# dc-HN
# LN-dc
# HN-end
# kj-sa
# kj-HN
# kj-dc
#
# The 19 paths through it are as follows:
#
# start,HN,dc,HN,end
# start,HN,dc,HN,kj,HN,end
# start,HN,dc,end
# start,HN,dc,kj,HN,end
# start,HN,end
# start,HN,kj,HN,dc,HN,end
# start,HN,kj,HN,dc,end
# start,HN,kj,HN,end
# start,HN,kj,dc,HN,end
# start,HN,kj,dc,end
# start,dc,HN,end
# start,dc,HN,kj,HN,end
# start,dc,end
# start,dc,kj,HN,end
# start,kj,HN,dc,HN,end
# start,kj,HN,dc,end
# start,kj,HN,end
# start,kj,dc,HN,end
# start,kj,dc,end
#
# Finally, this even larger example has 226 paths through it:
#
# fs-end
# he-DX
# fs-he
# start-DX
# pj-DX
# end-zg
# zg-sl
# zg-pj
# pj-he
# RW-he
# fs-DX
# pj-RW
# zg-RW
# start-pj
# he-WI
# zg-he
# pj-fs
# start-RW
#
# How many paths through this cave system are there that visit small caves at
# most once?

# load data into a dictionary
caves = Dict{String, Vector{String}}()
foreach(eachline("input.txt")) do l
  (k, v) = split(l, '-')
  if haskey(caves, k)
    push!(caves[k], v)
  else
    caves[k] = [v]
  end
  if haskey(caves, v)
    push!(caves[v], k)
  else
    caves[v] = [k]
  end
end

# A dictionary of small caves names to an arbitrary index. Used to map small
# caves to a BitVector of visited caves.
small_cave_idxs = Dict([(k,i) for (i,k) in enumerate(filter(c -> all(islowercase, c), keys(caves)))])

# Implements a BFS of caves - returns the number of ways to reach end
function visit_cave(cave::String, visited::BitVector)
  # have we reached end?
  if cave == "end"
    return 1
  end

  # if cave is a small cave, mark it visited
  if haskey(small_cave_idxs, cave)
    visited = copy(visited)
    visited[small_cave_idxs[cave]] = true
  end

  # find adjacent large caves, or small caves we haven't yet visited
  adjacent_caves = filter(caves[cave]) do adjacent_cave
    !haskey(small_cave_idxs, adjacent_cave) || !visited[small_cave_idxs[adjacent_cave]]
  end

  # recurse
  foldl((a,c) -> a + visit_cave(c, visited), adjacent_caves; init=0)
end

# Do eet!
cnt = visit_cave("start", falses(length(small_cave_idxs)))
println("Result: ", cnt)