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

# ╔═╡ 3a4957ec-8723-11eb-22a0-8b35322596e2
html"""

Section 2.2

 Sampling and Random Variables

Lecture Video

"""

# ╔═╡ 06d2666a-8723-11eb-1395-0febdf3dc2a4
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"),

			Pkg.PackageSpec(name="Colors"),

			Pkg.PackageSpec(name="StatsBase"),

			Pkg.PackageSpec(name="Distributions")

			])

	using Plots, PlutoUI, Colors, Images, StatsBase, Distributions
	using Statistics

end

# ╔═╡ 0a70bca4-8723-11eb-1bcf-e9abb9b1ab75
PlutoUI.TableOfContents(aside=true)

# ╔═╡ 472a41d2-8724-11eb-31b3-0b81612f0083
md"""
## Julia: Useful tidbits
"""

# ╔═╡ aeb99f72-8725-11eb-2efd-d3e44686be03
md"""
The package that each function comes from is shown in brackets, unless it comes from `Base`.
"""

# ╔═╡ 4f9bd326-8724-11eb-2c9b-db1ac9464f1e
md"""

- Julia `Base` library (no `using` required):

  - `if...else...end`
  - `Dict`: Julia's dictionary type
  - `÷` or `div`:  integer division (type `\div` + <tab>)

  - `sum(S)`: sum of elements in the collection `S`, e.g. an array
  - `count(S)`: count the number of true elements of a Boolean collection

  - `rand(S)`: random sampling from a collection `S`

- `Statistics.jl` (pre-loaded standard library; just needs `using`)

  - `mean(S)`: calculate the mean of a collection `S`
  - `std(S)`:  calculate the standard deviation of a collection `S`

- `StatsBase.jl`:
  - `countmap`

- `Plots.jl`:

  - `histogram(x)`: Plot a histogram of data vector `x` [Plots]

  - `bar(d)`: Plot a bar graph of categorical data [Plots]

- `Colors.jl`:
  - `distinguishable_colors(n)`: Make `n` distinguishable colours [Colors]

"""

# ╔═╡ db2d25de-86b1-11eb-0c78-d1ee52e019ca
md"""
## Random sampling with `rand`
"""

# ╔═╡ e33fe4c8-86b1-11eb-1031-cf45717a3dc9
md"""
The `rand` function in Julia is quite versatile: it tries to generate, or **sample**, a random object from the argument that you pass in:
"""

# ╔═╡ f49191a2-86b1-11eb-3eab-b392ba058415
rand(1:6)

# ╔═╡ 1abda6c4-86b2-11eb-2aa3-4d1148bb52b7
rand([2, 3, 5, 7, 11])

# ╔═╡ 30b12f28-86b2-11eb-087b-8d50ec429b89
rand("MIT")

# ╔═╡ 4ce946c6-86b2-11eb-1820-0728798665ab
rand('a':'z')

# ╔═╡ fae3d138-8743-11eb-1014-b3a2a9b49aba
typeof('a':'z')

# ╔═╡ 6cdea3ae-86b2-11eb-107a-17bea3f54bc9
rand()   # random number between 0 and 1

# ╔═╡ 1c769d58-8744-11eb-3bd3-ab11ea1503ed
rand

# ╔═╡ 297fdfa0-8744-11eb-1934-9fe31e8be534
methods(rand);

# ╔═╡ 776ec3f2-86b3-11eb-0216-9b71d07e99f3
md"""
We can take random objects from a collection of objects of *any* type, for example:
"""

# ╔═╡ 5fcf8d4e-8744-11eb-080e-cba749004b08
distinguishable_colors(10)

# ╔═╡ 4898106a-8744-11eb-128a-35fec741e6b8
typeof(distinguishable_colors(10))

# ╔═╡ 0926366a-86b2-11eb-0f6d-31ae6981598c
rand(distinguishable_colors(3))   # from Colors.jl package

# ╔═╡ 7c8d7b72-86b2-11eb-2dd5-4f77bc5fb8ff
md"""
### Several random objects
"""

# ╔═╡ 2090b7f2-86b3-11eb-2a99-ed98800e1d63
md"""
To sample several random objects from the same collection, we could write an array comprehension:
"""

# ╔═╡ a7dff55c-86b2-11eb-330f-3d6279347095
[rand(1:6) for i in 1:10]

# ╔═╡ 2db33022-86b3-11eb-17dd-13c534ac9892
md"""
But in fact, just adding another argument to `rand` does the trick:
"""

# ╔═╡ 0de6f23e-86b2-11eb-39ff-318bbc4ecbcf
rand(1:6, 10)

# ╔═╡ 36c3da4a-86b3-11eb-0b2f-fffdde06fcd2
md"""
In fact, you can also generate not only random vectors, but also random matrices:
"""

# ╔═╡ 940c2bf6-86b2-11eb-0a5e-011abdd6352b
rand(1:6, 10, 12)

# ╔═╡ 5a4e7fc4-86b3-11eb-3376-0941b79574aa
rand(distinguishable_colors(5), 10, 10)

# ╔═╡ c433104e-86b3-11eb-20bb-af608bb281cc
md"""
We can also use random images:
"""

# ╔═╡ 78dc94e2-8723-11eb-1ff2-bb7104b62033
penny_image = load(download("https://www.usacoinbook.com/us-coins/lincoln-memorial-cent.jpg"))

# ╔═╡ bb1465c4-8723-11eb-1abc-bdb5a7028cf2
begin
	head = penny_image[:, 1:end÷2]
	tail = penny_image[:, end÷2:end]
end;

# ╔═╡ e04f3828-8723-11eb-3452-09f821391ad0
rand( [head, tail], 5, 5)

# ╔═╡ b7793f7a-8726-11eb-11d8-cd928f1a3645
md"""
## Uniform sampling
"""

# ╔═╡ ba80cc78-8726-11eb-2f33-e364f19295d8
md"""
So far, each use of `rand` has done **uniform sampling**, i.e. each possible object as the *same* probability.

Let's *count* heads and tails using the `countmap` function from the `StatsBase.jl` package:
"""

# ╔═╡ 8fe715e4-8727-11eb-2e7f-15b723bb8d9d
tosses = rand( ["head", "tail"], 10000)

# ╔═╡ 9da4e6f4-8727-11eb-08cb-d55e3bbff0e4
toss_counts = countmap(tosses)

# ╔═╡ 9647840c-8745-11eb-33fc-898ae351ddfe
"tail" => 3

# ╔═╡ 9d782aa6-8745-11eb-034e-0f8b01be987b
typeof("tail" => 3)

# ╔═╡ a693582c-8745-11eb-261b-ef79e503420e
toss_counts["tail"]

# ╔═╡ c3255d7a-8727-11eb-0421-d99faf27c7b0
md"""

We see that `countmap` returns a **dictionary** (`Dict`), which maps **keys** to **values**; we will say more about dictionaries in another chapter.
"""

# ╔═╡ e2037806-8727-11eb-3c36-1500f1dd545b
prob_tail = toss_counts["tail"] / length(tosses)

# ╔═╡ 28ff621c-8728-11eb-2ec0-2579cb313315
md"""
As we increase the number of tosses, we "expect" the probability to get closer to $1/2$.
"""

# ╔═╡ 57125768-8728-11eb-3c3b-1dd37e1ac189
md"""
## Tossing a weighted coin
"""

# ╔═╡ 6a439c78-8728-11eb-1969-27a19d254704
md"""
How could we model a coin that is **weighted**, so that it is more likely to come up heads? We want to assign a probability $p = 0.7$ to heads, and $q = 0.3$ to tails.
"""

# ╔═╡ 9f8ac08c-8728-11eb-10ad-f93ca225ce38
md"""
One way would be to generate random integers between 1 and 10 and assign heads to a subset of the possible results with the desired probability, e.g. 1:7 get heads, and 8:10 get tails:
"""

# ╔═╡ 062b400a-8729-11eb-16c5-235cef648edb
function simple_weighted_coin()
	outcome = if rand(1:10) ≤ 7
		"heads"
	else      # could have elseif
		"tails"
	end

	return outcome
end

# ╔═╡ 7c099606-8746-11eb-37fe-c3befde06e9d
function simple_weighted_coin2()
	if rand(1:10) ≤ 7
		"heads"
	else      # could have elseif
		"tails"
	end
end

# ╔═╡ 97cb3dde-8746-11eb-0b00-690f20cb26dc
result = for i in 1:10

end

# ╔═╡ 9d8cdc8e-8746-11eb-2b9a-b30a52026f09
result == nothing

# ╔═╡ 81a30c9e-8746-11eb-38c8-9be4f6ba2e80
simple_weighted_coin2()

# ╔═╡ 2a81ba88-8729-11eb-3dcb-1db26e468066
md"""
(Note that `if` statements have a return value in Julia.)
"""

# ╔═╡ 5ea5838a-8729-11eb-1749-030533fb0656
md"""
How could we generalise this to an arbitrary probability $p ∈ [0, 1]$?

We can generate a uniform floating-point number between 0 and 1, and check if it is less than $p$. This is called a **Bernoulli trial**.
"""

# ╔═╡ c9f21046-8746-11eb-27c6-910807240fd1
rand()

# ╔═╡ 806e6aae-8729-11eb-19ea-33722c60edf0
rand() < 0.314159

# ╔═╡ 90642564-8729-11eb-3cd7-b3d41a1553b4
md"""
Note that comparisons also return a value in Julia. Here we have switched from heads/tails to true/false as the output.
"""

# ╔═╡ 9a426a20-8729-11eb-0c0f-31e1d4dc91bc
md"""
Let's make that into a function:
"""

# ╔═╡ a4870b14-8729-11eb-20ee-e531d4a7108d
bernoulli(p) = rand() < p

# ╔═╡ 081e3796-8747-11eb-32ec-dfd998605737
bernoulli(0.7)

# ╔═╡ 008a40d2-872a-11eb-224d-5b3331f29c99
md"""
p = $(@bind p Slider(0.0:0.01:1.0, show_value=true, default=0.7))
"""

# ╔═╡ baed5908-8729-11eb-00e0-9f749406c30c
countmap( [bernoulli(p) for i in 1:1000] )

# ╔═╡ ed94eae8-86b3-11eb-3f1b-15c7a54903f5
md"""
## Bar charts and histograms
"""

# ╔═╡ f20504de-86b3-11eb-3125-3140e0e060b0
md"""
Once we have generated several random objects, it is natural to want to count **how many times** each one occurred.
"""

# ╔═╡ 7f1a3766-86b4-11eb-353f-b13acaf1503e
md"""
Let's roll a (fair) die 1000 times:
"""

# ╔═╡ 371838f8-86b4-11eb-1633-8d282e42a085
md"""
An obvious way to find the counts would be to run through the data looking for 1s, then run through again looking for 2s, etc.:
"""

# ╔═╡ 94688c1a-8747-11eb-13a3-eb36f731674c
rolls .== 1

# ╔═╡ ad701cdc-8747-11eb-3804-63a0fc881547
count(rolls .== 1)

# ╔═╡ 2405eb68-86b4-11eb-31b0-dff8e355d88e
counts = [count(rolls .== i) for i in 1:6]

# ╔═╡ 9e9d3556-86b5-11eb-3dfb-916e625da235
md"""
Note that this is *not* the most efficient algorithm!
"""

# ╔═╡ 90844738-8738-11eb-0604-3d23662152d9
md"""
We can plot **categorical data** using a **bar chart**, `bar` in Plots.jl. This counts each discrete item.
"""

# ╔═╡ 02d03642-86b4-11eb-365a-63ff61ddd3b5
rolls = rand(1:6, 100000)   # try modifying 100 by adding more zeros

# ╔═╡ 2d71fa88-86b5-11eb-0e55-35566c2246d7
begin
	bar(counts, alpha=0.5, leg=false, size=(500, 300))
	hline!([length(rolls) / 6], ls=:dash)
	title!("number of die rolls = $(length(rolls))")
	ylims!(0, length(rolls) / 3)
end

# ╔═╡ cb8a9762-86b1-11eb-0484-6b6cc8b1b14c
md"""
## Probability densities

### Rolling multiple dice
"""

# ╔═╡ d0c9814e-86b1-11eb-2f29-1d041bccc649
roll_dice(n) = sum( rand(1:12, n) )

# ╔═╡ b81b1090-8735-11eb-3a52-2dca4d4ed472
experiment() = roll_dice(n)

# experiment() = sum([randn()^2 for i in 1:n])

# ╔═╡ 7a16b674-86b7-11eb-3aa5-83712cdc8580
trials = 10^6

# ╔═╡ e8e811de-86b6-11eb-1cbf-6d4aeaee510a
data = [experiment() for t in 1:trials]

# ╔═╡ 2bfa712a-8738-11eb-3248-6f9bb93154e8
md"""
### Converging shape
"""

# ╔═╡ e4abcbf4-86b8-11eb-167a-d97c61e07837
data

# ╔═╡ 6c133ab6-86b7-11eb-15f6-7780da5afc31
md"""
n = $(@bind n Slider(1:50, show_value=true))
"""

# ╔═╡ 514f6be0-86b8-11eb-30c9-d1020f783afe
histogram(data, alpha=0.5, legend=false, bins=200, c=:lightsalmon1, title="n = $n")
# c = RGB(0.1, 0.2, 0.3))

# ╔═╡ a15fc456-8738-11eb-25bd-b15c2b16d461
md"""
Here we have switched from a bar chart to a **histogram**, which counts the number of items falling into a given range or **bin**. When $n$ is small this tends to look like a bar chart, but it looks like a "smooth bar chart" as $n$ gets larger, due to the **aggregation**.
"""

# ╔═╡ dd753568-8736-11eb-1f20-1b81110ae807
md"""
Does the above histogram look like a bell to you?
"""

# ╔═╡ 8ab9001a-8737-11eb-1009-5717fbe83af7
begin
	bell = load(download("https://encrypted-tbn0.gstatic.com/images?q=tbn:ANd9GcRmxRAIQt_L-X99A_4FoP3vsC-l_WHlC3TtAw&usqp=CAU"))

	bell[1:end*9÷10, :]
end

# ╔═╡ f8b2dd20-8737-11eb-1593-43659c693109
md"""
### Normalising the $y$-axis
"""

# ╔═╡ 89bb366a-8737-11eb-10a6-e754ee817f9a
md"""
Notice that the *shape* of the curve seems to converge to a bell-shaped curve, but the axes do not. What can we do about this?

We need to **normalise** so that the total shaded *area* is 1. We can use the histogram option `norm=true` to do so.
"""

# ╔═╡ 14f0a090-8737-11eb-0ccf-391249267401
histogram(data, alpha=0.5, legend=false, bins=50, norm=true,
			c=:lightsalmon1, title="n = $n", ylims=(0, 0.05))

# ╔═╡ e305467e-8738-11eb-1213-eb11aaebe151
md"""
### Normalising the $x$ axis
"""

# ╔═╡ e8341288-8738-11eb-27ae-0795fa7e4a7e
md"""
As we changed $n$ above, the range of values on the $x$ axis also changed. We want to find a way to rescale $x$ at the same time so that both axes and shape stop changing as we increase $n$.

We need to make sure that the data is centred in the same place -- we will choose 0. And we need to make sure that the width is the same -- we will divide by the standard deviation.
"""

# ╔═╡ 07e8ae34-873b-11eb-1df2-175392ac4678
σ = std(data)

# ╔═╡ 77616afe-873a-11eb-3f11-1bc417f53138
mean(data), std(data)

# ╔═╡ 16023bea-8739-11eb-1e32-79f2d006b093
normalised_data = ( data .- mean(data) ) ./ std(data)

# ╔═╡ bfc52118-873b-11eb-0fbb-952c60bc7cc2
histogram(normalised_data, bins=-5 - (1/(2σ)):(1/σ):5, norm=true,
	alpha=0.5, leg=false, ylim=(0, 0.41), size=(500, 300))

# ╔═╡ aa6126f8-873b-11eb-3b4a-0f96fe07b7fb
plot!(x -> exp(-x^2 / 2) / √(2π), lw=2, c=:red, alpha=0.5)

# ╔═╡ 308547c6-873d-11eb-3a42-833f8bf496ae
md"""
Note that in the limit, the data becomes *continuous*, no longer discrete. The probability of any particular value is 0 ! We then talk about a **probability density function**, $f(x)$; the integral of this function over the interval $[a, b]$ gives the probability of being in that interval.
"""

# ╔═╡ cb2fb68e-8749-11eb-29ea-9729ac0c63b4
md"""
### Options for plotting functions
"""

# ╔═╡ e0a1863e-8735-11eb-1182-1b3c59b1e05a
md"""
Other options for the `histogram` function options:
- `legend=false` turns off the legend (key), i.e. the plot labels
- `bins=50` specifies the number of **bins**; can also be a vector of bin edges
- `linetype=:stephist`: use steps instead of bars

- `alpha`: a general plot option specifying transparency (0=invisible; 1=opaque)
- `c` or `color`: a general plot option for the colour
- `lw`: line width (default=1)

There are several different ways to specify colours. [Here](http://juliagraphics.github.io/Colors.jl/stable/namedcolors/) is a list of named colours, but you can also specify `RGB(0.1, 0.2, 0.3)`.
"""

# ╔═╡ 900af0c8-86ba-11eb-2270-71b1869b9a1a
md"""
Note that `linetype=:stephist` will give a stepped version of the histogram:
"""

# ╔═╡ be6e4c00-873c-11eb-1413-5326aba54216
md"""
## Sampling from other distributions
"""

# ╔═╡ 9a1136c2-873c-11eb-124f-c3939972ce4a
md"""
dof = $(@bind dof Slider(1:50, show_value=true))
"""

# ╔═╡ e01b6f70-873c-11eb-04a1-ad8e86578982
chisq_data = rand( Chisq(dof), 100000 )

# ╔═╡ b5251f76-873c-11eb-38cb-7db300c8fe3c
histogram( chisq_data, norm=true, bins=100, size=(500, 300), leg=false, alpha=0.5,
	xlims=(0, 10*√(dof)))

# ╔═╡ da62fd1c-873c-11eb-0758-e7cb48e964f1
histogram( [ sum( randn().^2 for _=1:dof )  for _ = 1:100000], norm=true,
	alpha=0.5, leg=false)

# ╔═╡ Cell order:
# ╟─3a4957ec-8723-11eb-22a0-8b35322596e2
# ╠═06d2666a-8723-11eb-1395-0febdf3dc2a4
# ╠═0a70bca4-8723-11eb-1bcf-e9abb9b1ab75
# ╟─472a41d2-8724-11eb-31b3-0b81612f0083
# ╟─aeb99f72-8725-11eb-2efd-d3e44686be03
# ╟─4f9bd326-8724-11eb-2c9b-db1ac9464f1e
# ╟─db2d25de-86b1-11eb-0c78-d1ee52e019ca
# ╟─e33fe4c8-86b1-11eb-1031-cf45717a3dc9
# ╠═f49191a2-86b1-11eb-3eab-b392ba058415
# ╠═1abda6c4-86b2-11eb-2aa3-4d1148bb52b7
# ╠═30b12f28-86b2-11eb-087b-8d50ec429b89
# ╠═4ce946c6-86b2-11eb-1820-0728798665ab
# ╠═fae3d138-8743-11eb-1014-b3a2a9b49aba
# ╠═6cdea3ae-86b2-11eb-107a-17bea3f54bc9
# ╠═1c769d58-8744-11eb-3bd3-ab11ea1503ed
# ╠═297fdfa0-8744-11eb-1934-9fe31e8be534
# ╟─776ec3f2-86b3-11eb-0216-9b71d07e99f3
# ╠═5fcf8d4e-8744-11eb-080e-cba749004b08
# ╠═4898106a-8744-11eb-128a-35fec741e6b8
# ╠═0926366a-86b2-11eb-0f6d-31ae6981598c
# ╟─7c8d7b72-86b2-11eb-2dd5-4f77bc5fb8ff
# ╟─2090b7f2-86b3-11eb-2a99-ed98800e1d63
# ╠═a7dff55c-86b2-11eb-330f-3d6279347095
# ╟─2db33022-86b3-11eb-17dd-13c534ac9892
# ╠═0de6f23e-86b2-11eb-39ff-318bbc4ecbcf
# ╟─36c3da4a-86b3-11eb-0b2f-fffdde06fcd2
# ╠═940c2bf6-86b2-11eb-0a5e-011abdd6352b
# ╠═5a4e7fc4-86b3-11eb-3376-0941b79574aa
# ╟─c433104e-86b3-11eb-20bb-af608bb281cc
# ╠═78dc94e2-8723-11eb-1ff2-bb7104b62033
# ╠═bb1465c4-8723-11eb-1abc-bdb5a7028cf2
# ╠═e04f3828-8723-11eb-3452-09f821391ad0
# ╟─b7793f7a-8726-11eb-11d8-cd928f1a3645
# ╟─ba80cc78-8726-11eb-2f33-e364f19295d8
# ╠═8fe715e4-8727-11eb-2e7f-15b723bb8d9d
# ╠═9da4e6f4-8727-11eb-08cb-d55e3bbff0e4
# ╠═9647840c-8745-11eb-33fc-898ae351ddfe
# ╠═9d782aa6-8745-11eb-034e-0f8b01be987b
# ╠═a693582c-8745-11eb-261b-ef79e503420e
# ╟─c3255d7a-8727-11eb-0421-d99faf27c7b0
# ╠═e2037806-8727-11eb-3c36-1500f1dd545b
# ╟─28ff621c-8728-11eb-2ec0-2579cb313315
# ╟─57125768-8728-11eb-3c3b-1dd37e1ac189
# ╟─6a439c78-8728-11eb-1969-27a19d254704
# ╟─9f8ac08c-8728-11eb-10ad-f93ca225ce38
# ╠═062b400a-8729-11eb-16c5-235cef648edb
# ╠═7c099606-8746-11eb-37fe-c3befde06e9d
# ╠═97cb3dde-8746-11eb-0b00-690f20cb26dc
# ╠═9d8cdc8e-8746-11eb-2b9a-b30a52026f09
# ╠═81a30c9e-8746-11eb-38c8-9be4f6ba2e80
# ╟─2a81ba88-8729-11eb-3dcb-1db26e468066
# ╟─5ea5838a-8729-11eb-1749-030533fb0656
# ╠═c9f21046-8746-11eb-27c6-910807240fd1
# ╠═806e6aae-8729-11eb-19ea-33722c60edf0
# ╟─90642564-8729-11eb-3cd7-b3d41a1553b4
# ╟─9a426a20-8729-11eb-0c0f-31e1d4dc91bc
# ╠═a4870b14-8729-11eb-20ee-e531d4a7108d
# ╠═081e3796-8747-11eb-32ec-dfd998605737
# ╟─008a40d2-872a-11eb-224d-5b3331f29c99
# ╠═baed5908-8729-11eb-00e0-9f749406c30c
# ╟─ed94eae8-86b3-11eb-3f1b-15c7a54903f5
# ╟─f20504de-86b3-11eb-3125-3140e0e060b0
# ╟─7f1a3766-86b4-11eb-353f-b13acaf1503e
# ╟─371838f8-86b4-11eb-1633-8d282e42a085
# ╠═94688c1a-8747-11eb-13a3-eb36f731674c
# ╠═ad701cdc-8747-11eb-3804-63a0fc881547
# ╠═2405eb68-86b4-11eb-31b0-dff8e355d88e
# ╟─9e9d3556-86b5-11eb-3dfb-916e625da235
# ╟─90844738-8738-11eb-0604-3d23662152d9
# ╠═02d03642-86b4-11eb-365a-63ff61ddd3b5
# ╠═2d71fa88-86b5-11eb-0e55-35566c2246d7
# ╟─cb8a9762-86b1-11eb-0484-6b6cc8b1b14c
# ╠═d0c9814e-86b1-11eb-2f29-1d041bccc649
# ╠═b81b1090-8735-11eb-3a52-2dca4d4ed472
# ╠═7a16b674-86b7-11eb-3aa5-83712cdc8580
# ╠═e8e811de-86b6-11eb-1cbf-6d4aeaee510a
# ╟─2bfa712a-8738-11eb-3248-6f9bb93154e8
# ╠═e4abcbf4-86b8-11eb-167a-d97c61e07837
# ╟─6c133ab6-86b7-11eb-15f6-7780da5afc31
# ╠═514f6be0-86b8-11eb-30c9-d1020f783afe
# ╟─a15fc456-8738-11eb-25bd-b15c2b16d461
# ╟─dd753568-8736-11eb-1f20-1b81110ae807
# ╠═8ab9001a-8737-11eb-1009-5717fbe83af7
# ╟─f8b2dd20-8737-11eb-1593-43659c693109
# ╟─89bb366a-8737-11eb-10a6-e754ee817f9a
# ╠═14f0a090-8737-11eb-0ccf-391249267401
# ╟─e305467e-8738-11eb-1213-eb11aaebe151
# ╟─e8341288-8738-11eb-27ae-0795fa7e4a7e
# ╠═07e8ae34-873b-11eb-1df2-175392ac4678
# ╠═77616afe-873a-11eb-3f11-1bc417f53138
# ╠═16023bea-8739-11eb-1e32-79f2d006b093
# ╠═bfc52118-873b-11eb-0fbb-952c60bc7cc2
# ╠═aa6126f8-873b-11eb-3b4a-0f96fe07b7fb
# ╟─308547c6-873d-11eb-3a42-833f8bf496ae
# ╟─cb2fb68e-8749-11eb-29ea-9729ac0c63b4
# ╟─e0a1863e-8735-11eb-1182-1b3c59b1e05a
# ╟─900af0c8-86ba-11eb-2270-71b1869b9a1a
# ╟─be6e4c00-873c-11eb-1413-5326aba54216
# ╟─9a1136c2-873c-11eb-124f-c3939972ce4a
# ╠═e01b6f70-873c-11eb-04a1-ad8e86578982
# ╠═b5251f76-873c-11eb-38cb-7db300c8fe3c
# ╠═da62fd1c-873c-11eb-0758-e7cb48e964f1


text:
using Documenter, TSVD

DocMeta.setdocmeta!(
    TSVD,
    :DocTestSetup,
    :(using TSVD, MatrixDepot; matrixdepot("LPnetlib/lp_osa_30"));
    recursive=true)

makedocs(sitename="TSVD Documentation")

deploydocs(
    repo = "github.com/JuliaLinearAlgebra/TSVD.jl.git",
)

ext: jl


text:
module SCP_VB

include("../Misc.jl")
using .Misc

using Distributions, SpecialFunctions

export standard_VB

function standard_VB(X::Array{ℜ,3}, X⁺::Array{ℜ,2}, R::Ƶ ; μ::ℜ=3.5, γ::ℜ=0.1, EPOCHS::Ƶ=1, ϵ::ℜ=1e-16) where {ℜ<:Real, Ƶ<:Int}
    I::Ƶ, J::Ƶ, K::Ƶ = size(X)

    a::ℜ, b::ℜ = I*J*K*γ, γ/μ
    α_R, α_KR, α_IKR, α_JKR = fill(a/R,R), fill(a/(K*R),K,R), fill(a/(I*K*R),I,K,R), fill(a/(J*K*R),J,K,R)

    log_λ = log(rand(Gamma(a,1.0/b)))
    log_θ_R = log.(rand(Dirichlet(α_R)))
    log_θ_KR = reshape([log(θ_kr) for r=1:R for θ_kr=rand(Dirichlet(α_KR[:,r] .+ 1.0/K ))],K,R)
    log_θ_IKR = reshape([log(θ_ikr) for k=1:K, r=1:R for θ_ikr=rand(Dirichlet(α_IKR[:,k,r] .+ 1.0/I))],I,K,R)
    log_θ_JKR = reshape([log(θ_jkr) for k=1:K, r=1:R for θ_jkr=rand(Dirichlet(α_JKR[:,k,r] .+ 1.0/J))],J,K,R)

    s, S_R, S_KR, S_IKR, S_JKR, S⁺ = zeros(ℜ,K,R), zeros(ℜ,R), zeros(ℜ,K,R), zeros(ℜ,I,K,R), zeros(ℜ,J,K,R), sum(X⁺)
    π, log_π, log_ρ = Array{ℜ}(undef,K,R), Array{ℜ}(undef,K,R), Array{ℜ}(undef,K,R)
    log_ρ_ij::ℜ, log_ρ_ijk::Array{ℜ} = 0.0, Array{ℜ}(undef,K)

    ELBO::Array{ℜ} = zeros(ℜ,EPOCHS)
    X_full = copy(X)

    for eph=1:EPOCHS
        S_R .= 0.0
        S_KR .= 0.0
        S_IKR .= 0.0
        S_JKR .= 0.0

        for j=1:J, i=1:I #order of traversal is important

            log_ρ .= log_λ .+ log_θ_R' .+ log_θ_KR .+ log_θ_IKR[i,:,:] .+ log_θ_JKR[j,:,:]

            if isnan(X[i,j,1])
                log_ρ_ij = logsumexp(log_ρ)
                log_π .= log_ρ .- log_ρ_ij
                π .= exp.(log_π)
                s = X⁺[i,j] .* π
                X_full[i,j,:] .= sum(s,dims=2)[:,1]
                ELBO[eph] += X⁺[i,j] * log_ρ_ij - lgamma(X⁺[i,j] + 1.0)
            else
                log_ρ_ijk = logsumexp(log_ρ;dims=2)
                log_π .= log_ρ .- log_ρ_ijk
                π .= exp.(log_π)
                s = X[i,j,:] .* π
                ELBO[eph] += sum(X[i,j,:] .* log_ρ_ijk .- lgamma.(X[i,j,:] .+ 1.0))
            end

            S_R .+= sum(s,dims=1)[1,:]
            S_KR .+= s
            S_JKR[j,:,:] .+= s
            S_IKR[i,:,:] .+= s

        end

        ELBO[eph] += a*log(b) - (S⁺+a)*log(b+1.0) + sum(lgamma.(α_KR)) - sum(lgamma.(α_KR .+ S_KR))
        ELBO[eph] += sum(lgamma.(α_IKR .+ S_IKR)) - sum(lgamma.(α_IKR))
        ELBO[eph] += sum(lgamma.(α_JKR .+ S_JKR)) - sum(lgamma.(α_JKR))
        ELBO[eph] -= S⁺ * log_λ + sum(S_R .* log_θ_R) + sum(S_KR .* log_θ_KR)
        ELBO[eph] -= sum(S_JKR .* log_θ_JKR) + sum(S_IKR .* log_θ_IKR)

        log_λ = digamma(S⁺+a) - log(b+1.0)
        log_θ_R .= digamma.(S_R.+α_R) .- digamma(S⁺+a)
        log_θ_KR .= digamma.(S_KR.+α_KR) .- digamma.(S_R.+α_R)'
        log_θ_IKR .= digamma.(S_IKR.+α_IKR) .- reshape(digamma.(S_KR.+α_KR),1,K,R)
        log_θ_JKR .= digamma.(S_JKR.+α_JKR) .- reshape(digamma.(S_KR.+α_KR),1,K,R)

    end
    return ELBO, X_full, log_λ, (log_θ_R, log_θ_KR, log_θ_IKR, log_θ_JKR)
end

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
mutable struct DummyGraph <: AbstractGraph{Int} end
mutable struct DummyDiGraph <: AbstractGraph{Int} end
mutable struct DummyEdge <: AbstractEdge{Int} end

@testset "Interface" begin
    dummygraph = DummyGraph()
    dummydigraph = DummyDiGraph()
    dummyedge = DummyEdge()

    @test_throws LightGraphs.NotImplementedError is_directed(DummyGraph)
    @test_throws LightGraphs.NotImplementedError zero(DummyGraph)

    for edgefun in [src, dst, Pair, Tuple, reverse]
        @test_throws LightGraphs.NotImplementedError edgefun(dummyedge)
    end

    for edgefun2edges in [==]
        @test_throws LightGraphs.NotImplementedError edgefun2edges(dummyedge, dummyedge)
     end

    for graphfunbasic in [
        nv, ne, vertices, edges, is_directed,
        edgetype, eltype
    ]
        @test_throws LightGraphs.NotImplementedError graphfunbasic(dummygraph)
    end

    for graphfun1int in [
        has_vertex, inneighbors, outneighbors
    ]
        @test_throws LightGraphs.NotImplementedError graphfun1int(dummygraph, 1)
    end
    for graphfunedge in [
        has_edge,
      ]
        @test_throws LightGraphs.NotImplementedError graphfunedge(dummygraph, dummyedge)
        @test_throws LightGraphs.NotImplementedError graphfunedge(dummygraph, 1, 2)
    end

    # Implementation error
    impl_error = LightGraphs.NotImplementedError(edges)
    @test impl_error isa LightGraphs.NotImplementedError{typeof(edges)}
    io = IOBuffer()
    Base.showerror(io, impl_error)
    @test String(take!(io)) == "method $edges not implemented."

end # testset


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id: condenser.mso 555 2008-07-18 19:01:13Z rafael $
#*--------------------------------------------------------------------
type condenserSubcooled
	condenserSubcooled()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Components"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Vapour inlet stream",
				:PosX=>0.16,
				:PosY=>0,
				:Symbol=>"_{in}^{Vapour}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Liquid outlet stream",
				:PosX=>0.53,
				:PosY=>1,
				:Symbol=>"_{out}^{Liquid}"
			)),
			temperature(Dict{Symbol,Any}(
				:Brief=>"Condensate temperature (subcooled)",
				:Symbol=>"T_{sub}"
			)),
			temp_delta(Dict{Symbol,Any}(
				:Brief=>"Subcooling Degree",
				:Symbol=>"\\Delta T_{sub}"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Calculated condenser duty for desired subcooling",
				:Protected=>true,
				:Symbol=>"Q_{cond}"
			)),
			[
				:(InletVapour.F = OutletLiquid.F),
				:(InletVapour.z = OutletLiquid.z),
				:(OutletLiquid.P = InletVapour.P),
				:(OutletLiquid.T = T_sub),
				:(SubcoolingDegree = InletVapour.T - T_sub),
				:(OutletLiquid.h = PP.LiquidEnthalpy(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z)),
				:(CondenserDuty = OutletLiquid.F*OutletLiquid.h - InletVapour.F*InletVapour.h),
				:(OutletLiquid.v = 0),
			],
			[
				"Molar Flow Balance","Molar Composition Balance","Pressure Drop","Subcooled Temperature","Degree of subcooling","Liquid enthalpy","Condenser Duty","Vapourisation Fraction",
			],
			[:PP,:NComp,],
			[:InletVapour,:OutletLiquid,:T_sub,:SubcoolingDegree,:CondenserDuty,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	InletVapour::stream
	OutletLiquid::stream
	T_sub::temperature
	SubcoolingDegree::temp_delta
	CondenserDuty::power
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export condenserSubcooled
function setEquationFlow(in::condenserSubcooled)
	addEquation(1)
	addEquation(2)
	#"Energy Balance"
	#InletVapour.F*InletVapour.h  + InletQ = OutletLiquid.F*OutletLiquid.h;
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
end
function atributes(in::condenserSubcooled,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/CondenserSteady"
	fields[:Brief]="Model of a  Steady State total condenser with specified outlet temperature conditions."
	fields[:Info]="A simple model of a Steady State total condenser with specified temperature (or subcooling degree), with a real
calculation of the output stream enthalpy. The subcooling degree is considered to be the difference between the
inlet vapour and the outlet liquid temperatures.

== ASSUMPTIONS ==
* perfect mixing of both phases;
* saturated vapour at the Inlet;
* no thermodynamics equilibrium;
* no pressure drop in the condenser.

== SPECIFY ==
* the InletVapour stream;
* the subcooled temperature OR the the degree of subcooling.

"
	drive!(fields,_)
	return fields
end
condenserSubcooled(_::Dict{Symbol,Any})=begin
	newModel=condenserSubcooled()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
__precompile__()

module IntervalUnions

# Imports

import Base: ==, ∈, ⊆, ∩, ∪, isless, show, string, setdiff, symdiff, empty

# Exports
export
	OrderedPair, Interval, IntervalUnion,

	# Query
	==, ∈, ⊆, isless, left, right, empty, limits, disjoint,

	# Display
	show, string,

	# Size
	number_of_components, cardinal, inverse_cardinal,

	# Operations
	∩, ∪, complement, jaccard, overlap_coefficient, dice_coefficient,
	setdiff, symdiff, compact, compactness, restricted_complement,
	superset, super_complement, complement_cardinal, sample

# Includes

include("base.jl")

end


text:
using DelimitedFiles

# doing a hello world, because why not
println("hello world")

# reading the input
data = readdlm("input_day1",Int)

#would give the tuples whith numbers that summ gives 2020
#since the for is going through the whole list, it will give
#2 tuples following that condition
#I'm ok with that for now :)

numbers = [(i,j) for i in data for j in data if i+j==2020]

print("part one!")
println(numbers)
println(numbers[1][1] * numbers[1][2])

numbers = [(i,j,k) for i in data for j in data for k in data if i+j+k==2020]
println("part2")
println(numbers)
println(numbers[1][1]*numbers[1][2]*numbers[1][3])


text:
using Test
using LoopVectorization
using LinearAlgebra

@testset "LoopVectorization.jl" begin

    function logsumexp!(r::AbstractArray{T}, x::AbstractArray{T}) where {T}
        n = length(x)
        length(r) == n || throw(DimensionMismatch())
        isempty(x) && return -T(Inf)
        1 == stride(r,1) == stride(x,1) || throw(error("Arrays not strided"))

        u = maximum(x)                                       # max value used to re-center
        abs(u) == Inf && return any(isnan, x) ? T(NaN) : u   # check for non-finite values
        s = zero(T)
        @avx for i = 1:n
            tmp = exp(x[i] - u)
            r[i] = tmp
            s += tmp
        end

        invs = inv(s)
        r .*= invs

        return log1p(s-1) + u
    end

    x = collect(1:1_000) ./ 10;
    r = similar(x);

    @test logsumexp!(r, x) ≈ 102.35216846104409

    @testset "GEMM" begin
        using LoopVectorization, Test
        U, T = LoopVectorization.VectorizationBase.REGISTER_COUNT == 16 ? (3, 4) : (4, 4)
        AmulBq1 = :(for m ∈ 1:size(A,1), n ∈ 1:size(B,2)
                    C[m,n] = zeroB
               for k ∈ 1:size(A,2)
                   C[m,n] += A[m,k] * B[k,n]
               end
               end)
        lsAmulB1 = LoopVectorization.LoopSet(AmulBq1);
        @test LoopVectorization.choose_order(lsAmulB1) == (Symbol[:n,:m,:k], :m, U, T)
        AmulBq2 = :(for m ∈ 1:M, n ∈ 1:N
               C[m,n] = zero(eltype(B))
               for k ∈ 1:K
                    C[m,n] += A[m,k] * B[k,n]
               end
                end)
        lsAmulB2 = LoopVectorization.LoopSet(AmulBq2);
        @test LoopVectorization.choose_order(lsAmulB2) == (Symbol[:n,:m,:k], :m, U, T)
        AmulBq3 = :(for m ∈ 1:size(A,1), n ∈ 1:size(B,2)
                ΔCₘₙ = zero(eltype(C))
                for k ∈ 1:size(A,2)
                    ΔCₘₙ += A[m,k] * B[k,n]
                end
                C[m,n] += ΔCₘₙ
           end)
        lsAmulB3 = LoopVectorization.LoopSet(AmulBq3);
        @test LoopVectorization.choose_order(lsAmulB3) == (Symbol[:n,:m,:k], :m, U, T)

        function AmulB!(C, A, B)
            C .= 0
            for k ∈ 1:size(A,2), j ∈ 1:size(B,2)
                @simd ivdep for i ∈ 1:size(A,1)
                    @inbounds C[i,j] += A[i,k] * B[k,j]
                end
            end
        end
        function AmulBavx1!(C, A, B)
            @avx for m ∈ 1:size(A,1), n ∈ 1:size(B,2)
                Cₘₙ = zero(eltype(C))
                for k ∈ 1:size(A,2)
                    Cₘₙ += A[m,k] * B[k,n]
                end
                C[m,n] = Cₘₙ
            end
        end
        function AmulBavx2!(C, A, B)
            z = zero(eltype(C))
            @avx for m ∈ 1:size(A,1), n ∈ 1:size(B,2)
                C[m,n] = z
                for k ∈ 1:size(A,2)
                    C[m,n] += A[m,k] * B[k,n]
                end
            end
        end
        function AmulBavx3!(C, A, B)
            @avx for m ∈ 1:size(A,1), n ∈ 1:size(B,2)
                C[m,n] = zero(eltype(C))
                for k ∈ 1:size(A,2)
                    C[m,n] += A[m,k] * B[k,n]
                end
            end
        end
        function AmuladdBavx!(C, A, B, factor = 1)
            @avx for m ∈ 1:size(A,1), n ∈ 1:size(B,2)
                ΔCₘₙ = zero(eltype(C))
                for k ∈ 1:size(A,2)
                    ΔCₘₙ += A[m,k] * B[k,n]
                end
                C[m,n] += ΔCₘₙ * factor
            end
        end

        # function AtmulB!(C, A, B)
        #     for j ∈ 1:size(C,2), i ∈ 1:size(C,1)
        #         Cᵢⱼ = zero(eltype(C))
        #         @simd ivdep for k ∈ 1:size(A,1)
        #             @inbounds Cᵢⱼ += A[k,i] * B[k,j]
        #         end
        #         C[i,j] = Cᵢⱼ
        #     end
        # end
        AtmulBq = :(for n ∈ 1:size(C,2), m ∈ 1:size(C,1)
                Cₘₙ = zero(eltype(C))
                for k ∈ 1:size(A,1)
                    Cₘₙ += A[k,m] * B[k,n]
                end
                C[m,n] = Cₘₙ
            end)
        lsAtmulB = LoopVectorization.LoopSet(AtmulBq);
        # LoopVectorization.choose_order(lsAtmulB)
        @test LoopVectorization.choose_order(lsAtmulB) == (Symbol[:m,:n,:k], :k, U, T)

        function AtmulBavx!(C, A, B)
            @avx for n ∈ 1:size(C,2), m ∈ 1:size(C,1)
                Cₘₙ = zero(eltype(C))
                for k ∈ 1:size(A,1)
                    Cₘₙ += A[k,m] * B[k,n]
                end
                C[m,n] = Cₘₙ
            end
        end
        function rank2AmulB!(C, Aₘ, Aₖ, B)
            @inbounds for m ∈ 1:size(C,1), n ∈ 1:size(C,2)
                Cₘₙ = zero(eltype(C))
                @fastmath for k ∈ 1:size(B,1)
                    Cₘₙ += (Aₘ[m,1]*Aₖ[1,k]+Aₘ[m,2]*Aₖ[2,k]) * B[k,n]
                end
                C[m,n] = Cₘₙ
            end
        end
        r2ambq = :(for m ∈ 1:size(C,1), n ∈ 1:size(C,2)
                Cₘₙ = zero(eltype(C))
                for k ∈ 1:size(B,1)
                    Cₘₙ += (Aₘ[m,1]*Aₖ[1,k]+Aₘ[m,2]*Aₖ[2,k]) * B[k,n]
                end
                C[m,n] = Cₘₙ
                   end)
        lsr2amb = LoopVectorization.LoopSet(r2ambq)
        function rank2AmulBavx!(C, Aₘ, Aₖ, B)
            @avx for m ∈ 1:size(C,1), n ∈ 1:size(C,2)
                Cₘₙ = zero(eltype(C))
                for k ∈ 1:size(B,1)
                    Cₘₙ += (Aₘ[m,1]*Aₖ[1,k]+Aₘ[m,2]*Aₖ[2,k]) * B[k,n]
                end
                C[m,n] = Cₘₙ
            end
        end

        function mulCAtB_2x2block!(C, A, B)
            M, N = size(C); K = size(B,1)
            @assert size(C, 1) == size(A, 2)
            @assert size(C, 2) == size(B, 2)
            @assert size(A, 1) == size(B, 1)
            T = eltype(C)
            if mod(M, 2) == 0 && mod(N, 2) == 0
	        for m ∈ 1:2:M
	    	    m1 = m + 1
	    	    for n ∈ 1:2:N
	    		n1 = n + 1
		    	C11, C21, C12, C22 = zero(T), zero(T), zero(T), zero(T)
		    	@avx for k ∈ 1:K
		    	    C11 += A[k,m] * B[k,n]
		    	    C21 += A[k,m1] * B[k,n]
		    	    C12 += A[k,m] * B[k,n1]
		    	    C22 += A[k,m1] * B[k,n1]
		    	end
		    	C[m,n] = C11
		    	C[m1,n] = C21
		    	C[m,n1] = C12
		    	C[m1,n1] = C22
		    end
	        end
	    else
		@inbounds for n ∈ 1:N, m ∈ 1:M
	    	    Cmn = 0.0
	    	    @inbounds for k ∈ 1:K
	    		Cmn += A[k,m] * B[k,n]
	    	    end
	    	    C[m,n] = Cmn
	        end
	    end
            return C
        end

        for T ∈ (Float32, Float64, Int32, Int64)
            M, K, N = 72, 75, 68;
            TC = sizeof(T) == 4 ? Float32 : Float64
            R = T <: Integer ? (T(1):T(1000)) : T
            C = Matrix{TC}(undef, M, N);
            A = rand(R, M, K); B = rand(R, K, N);
            C2 = similar(C);
            AmulB!(C2, A, B)
            AmulBavx1!(C, A, B)
            @test C ≈ C2
            fill!(C, 999.99); AmulBavx2!(C, A, B)
            @test C ≈ C2
            fill!(C, 999.99); AmulBavx3!(C, A, B)
            @test C ≈ C2
            fill!(C, 0.0); AmuladdBavx!(C, A, B)
            @test C ≈ C2
            AmuladdBavx!(C, A, B)
            @test C ≈ 2C2
            AmuladdBavx!(C, A, B, -1)
            @test C ≈ C2
            At = copy(A');
            fill!(C, 9999.999); AtmulBavx!(C, At, B)
            @test C ≈ C2
            fill!(C, 9999.999); mulCAtB_2x2block!(C, At, B);
            @test C ≈ C2
            Aₘ= rand(R, M, 2); Aₖ = rand(R, 2, K);
            rank2AmulBavx!(C, Aₘ, Aₖ, B)
            rank2AmulB!(C2, Aₘ, Aₖ, B)
            @test C ≈ C2
        end
    end

    @testset "dot" begin
        dotq = :(for i ∈ eachindex(a,b)
                 s += a[i]*b[i]
                 end)
        lsdot = LoopVectorization.LoopSet(dotq);
        @test LoopVectorization.choose_order(lsdot) == (Symbol[:i], :i, 4, -1)

        function mydot(a, b)
            s = zero(eltype(a))
            @inbounds @simd for i ∈ eachindex(a,b)
                s += a[i]*b[i]
            end
            s
        end
        function mydotavx(a, b)
            s = zero(eltype(a))
            @avx for i ∈ eachindex(a,b)
                s += a[i]*b[i]
            end
            s
        end

        selfdotq = :(for i ∈ eachindex(a)
                     s += a[i]*a[i]
                     end)
        lsselfdot = LoopVectorization.LoopSet(selfdotq);
        @test LoopVectorization.choose_order(lsselfdot) == (Symbol[:i], :i, 8, -1)

        function myselfdot(a)
            s = zero(eltype(a))
            @inbounds @simd for i ∈ eachindex(a)
                s += a[i]*a[i]
            end
            s
        end
        function myselfdotavx(a)
            s = zero(eltype(a))
            @avx for i ∈ eachindex(a)
                s += a[i]*a[i]
            end
            s
        end

        # a = rand(400);
        for T ∈ (Float32, Float64)
            a = rand(T, 100); b = rand(T, 100);
            @test mydotavx(a,b) ≈ mydot(a,b)
            @test myselfdotavx(a) ≈ myselfdot(a)
        end

    end

    @testset "Special Functions" begin
        vexpq = :(for i ∈ eachindex(a)
                  b[i] = exp(a[i])
                  end)
        lsvexp = LoopVectorization.LoopSet(vexpq);
        @test LoopVectorization.choose_order(lsvexp) == (Symbol[:i], :i, 1, -1)

        function myvexp!(b, a)
            @inbounds for i ∈ eachindex(a)
                b[i] = exp(a[i])
            end
        end
        function myvexpavx!(b, a)
            @avx for i ∈ eachindex(a)
                b[i] = exp(a[i])
            end
        end

        vexpsq = :(for i ∈ eachindex(a)
                   s += exp(a[i])
                   end)
        lsvexps = LoopVectorization.LoopSet(vexpsq);
        @test LoopVectorization.choose_order(lsvexps) == (Symbol[:i], :i, 1, -1)

        function myvexp(a)
            s = zero(eltype(a))
            @inbounds for i ∈ eachindex(a)
                s += exp(a[i])
            end
            s
        end
        function myvexpavx(a)
            s = zero(eltype(a))
            @avx for i ∈ eachindex(a)
                s += exp(a[i])
            end
            s
        end

        for T ∈ (Float32, Float64)
            a = randn(T, 127);
            b1 = similar(a);
            b2 = similar(a);

            myvexp!(b1, a)
            myvexpavx!(b2, a)
            @test b1 ≈ b2
            @test myvexp(a) ≈ myvexpavx(a)
            @test b1 ≈ @avx exp.(a)
        end
    end

    @testset "GEMV" begin
        gemvq = :(for i ∈ eachindex(y)
                  yᵢ = 0.0
                  for j ∈ eachindex(x)
                  yᵢ += A[i,j] * x[j]
                  end
                  y[i] = yᵢ
                  end)
        lsgemv = LoopVectorization.LoopSet(gemvq);
        @test LoopVectorization.choose_order(lsgemv) == (Symbol[:i, :j], :i, 8, -1)

        function mygemv!(y, A, x)
            @inbounds for i ∈ eachindex(y)
                yᵢ = zero(eltype(y))
                @simd for j ∈ eachindex(x)
                    yᵢ += A[i,j] * x[j]
                end
                y[i] = yᵢ
            end
        end
        function mygemvavx!(y, A, x)
            @avx for i ∈ eachindex(y)
                yᵢ = zero(eltype(y))
                for j ∈ eachindex(x)
                    yᵢ += A[i,j] * x[j]
                end
                y[i] = yᵢ
            end
        end
        M, K = 51, 49
        for T ∈ (Float32, Float64)
            A = randn(T, M, K);
            x = randn(T, K);
            y1 = Vector{T}(undef, M); y2 = similar(y1);
            mygemv!(y1, A, x)
            mygemvavx!(y2, A, x)

            @test y1 ≈ y2
        end
    end

@testset "Miscellaneous" begin

    dot3q = :(for m ∈ 1:M, n ∈ 1:N
            s += x[m] * A[m,n] * y[n]
              end)
    lsdot3 = LoopVectorization.LoopSet(dot3q);
    LoopVectorization.choose_order(lsdot3)

    @static if VERSION < v"1.4"
        dot3(x, A, y) = dot(x, A * y)
    else
        dot3(x, A, y) = dot(x, A, y)
    end
    function dot3avx(x, A, y)
        M, N = size(A)
        s = zero(promote_type(eltype(x), eltype(A), eltype(y)))
        @avx for m ∈ 1:M, n ∈ 1:N
            s += x[m] * A[m,n] * y[n]
        end
        s
    end

    subcolq = :(for i ∈ 1:size(A,2), j ∈ eachindex(x)
                B[j,i] = A[j,i] - x[j]
                end)
    lssubcol = LoopVectorization.LoopSet(subcolq);
    @test LoopVectorization.choose_order(lssubcol) == (Symbol[:j,:i], :j, 4, -1)
    ## @avx is SLOWER!!!!
    ## need to fix!
    function mysubcol!(B, A, x)
        @inbounds for i ∈ 1:size(A,2)
            @simd for j ∈ eachindex(x)
                B[j,i] = A[j,i] - x[j]
            end
        end
    end
    function mysubcolavx!(B, A, x)
        @avx for i ∈ 1:size(A,2), j ∈ eachindex(x)
            B[j,i] = A[j,i] - x[j]
        end
    end

    colsumq = :(for i ∈ 1:size(A,2), j ∈ eachindex(x)
                x[j] += A[j,i] - 0.25
                end)
    lscolsum = LoopVectorization.LoopSet(colsumq);
    @test LoopVectorization.choose_order(lscolsum) == (Symbol[:j,:i], :j, 8, -1)

    # my colsum is wrong (by 0.25), but slightly more interesting
    function mycolsum!(x, A)
        @. x = 0
        @inbounds for i ∈ 1:size(A,2)
            @simd for j ∈ eachindex(x)
                x[j] += A[j,i] - 0.25
            end
        end
    end
    function mycolsumavx!(x, A)
        @avx for j ∈ eachindex(x)
            xⱼ = zero(eltype(x))
            for i ∈ 1:size(A,2)
                xⱼ += A[j,i] - 0.25
            end
            x[j] = xⱼ
        end
    end

    varq = :(for j ∈ eachindex(s²), i ∈ 1:size(A,2)
             δ = A[j,i] - x̄[j]
             s²[j] += δ*δ
             end)
    lsvar = LoopVectorization.LoopSet(varq);
    LoopVectorization.choose_order(lsvar)
    @test LoopVectorization.choose_order(lsvar) == (Symbol[:j,:i], :j, 5, -1)

    function myvar!(s², A, x̄)
        @. s² = 0
        @inbounds for i ∈ 1:size(A,2)
            @simd for j ∈ eachindex(s²)
                δ = A[j,i] - x̄[j]
                s²[j] += δ*δ
            end
        end
    end
    function myvaravx!(s², A, x̄)
        @avx for j ∈ eachindex(s²)
            s²ⱼ = zero(eltype(s²))
            x̄ⱼ = x̄[j]
            for i ∈ 1:size(A,2)
                δ = A[j,i] - x̄ⱼ
                s²ⱼ += δ*δ
            end
            s²[j] = s²ⱼ
        end
    end

    for T ∈ (Float32, Float64)
        A = randn(T, 199, 498);
        x = randn(T, size(A,1));
        B1 = similar(A); B2 = similar(A);

        mysubcol!(B1, A, x)
        mysubcolavx!(B2, A, x)

        @test B1 ≈ B2
        x1 = similar(x); x2 = similar(x);
        mycolsum!(x1, A)
        mycolsumavx!(x2, A)
        @test x1 ≈ x2

        x̄ = x1 ./ size(A,2);
        myvar!(x1, A, x̄)
        myvaravx!(x2, A, x̄)
        @test x1 ≈ x2

        M, N = 47, 73;
        x = rand(T, M); A = rand(T, M, N); y = rand(T, N);
        @test dot3avx(x, A, y) ≈ dot3(x, A, y)

    end
end

@testset "broadcast" begin
    M, N = 37, 47
    # M = 77;
    for T ∈ (Float32, Float64)

        a = rand(T,100,100,100);
        b = rand(T,100,100,1);
        bl = LowDimArray{(true,true,false)}(b);
        br = reshape(b, (100,100));
        c1 = a .+ b;
        c2 = @avx a .+ bl;
        @test c1 ≈ c2
        fill!(c2, 99999.9);
        @avx c2 .= a .+ br;
        @test c1 ≈ c2
        br = reshape(b, (100,1,100));
        bl = LowDimArray{(true,false,true)}(br);
        @. c1 = a + br;
        fill!(c2, 99999.9);
        @avx @. c2 = a + bl;
        @test c1 ≈ c2
        br = reshape(b, (1,100,100));
        bl = LowDimArray{(false,true,true)}(br);
        @. c1 = a + br;
        fill!(c2, 99999.9);
        @avx @. c2 = a + bl;
        @test c1 ≈ c2

        a = rand(T, M); B = rand(T, M, N); c = rand(T, N); c′ = c';
        d1 =      @. a + B * c′;
        d2 = @avx @. a + B * c′;

        @test d1 ≈ d2
        @.      d1 = a + B * c′;
        @avx @. d2 = a + B * c′;
        @test d1 ≈ d2

        d3 = a .+ B * c;
        # no method matching _similar_for(::UnitRange{Int64}, ::Type{Any}, ::Product)
        d4 = @avx a .+ B *ˡ c;
        @test d3 ≈ d4

        fill!(d3, -1000.0);
        fill!(d4, 91000.0);

        d3 .= a .+ B * c;
        @avx d4 .= a .+ B *ˡ c;
        @test d3 ≈ d4

        fill!(d4, 91000.0);
        @avx @. d4 = a + B *ˡ c;
        @test d3 ≈ d4

        M, K, N = 77, 83, 57;
        A = rand(T,M,K); B = rand(T,K,N); C = rand(T,M,N);

        D1 = C .+ A * B;
        D2 = @avx C .+ A *ˡ B;
        @test D1 ≈ D2

        D3 = exp.(B');
        D4 = @avx exp.(B');
        @test D3 ≈ D4

        fill!(D3, -1e3); fill!(D4, 9e9);
        Bt = Transpose(B);
        @. D3 = exp(Bt);
        @avx @. D4 = exp(Bt);
        @test D3 ≈ D4

        D1 = similar(B); D2 = similar(B);
        D1t = Transpose(D1);
        D2t = Transpose(D2);
        @. D1t = exp(Bt);
        @avx @. D2t = exp(Bt);
        @test D1t ≈ D2t

        fill!(D1, -1e3);
        fill!(D2, 9e9);
        @. D1' = exp(Bt);
        lset = @avx @. D2' = exp(Bt);

        @test D1 ≈ D2

        a = rand(137);
        b1 = @avx @. 3*a + sin(a) + sqrt(a);
        b2 =      @. 3*a + sin(a) + sqrt(a);
        @test b1 ≈ b2
        three = 3; fill!(b1, -9999.999);
        @avx @. b1 = three*a + sin(a) + sqrt(a);
        @test b1 ≈ b2

        C = rand(100,10,10);
        D1 = C .^ 0.3;
        D2 = @avx C .^ 0.3;
        @test D1 ≈ D2

    end
end

@testset "map" begin
    foo(x, y) = exp(x) - sin(y)
    N = 37
    for T ∈ (Float32,Float64)
        a = rand(T, N); b = rand(T, N)
        c1 = map(foo, a, b)
        c2 = vmap(foo, a, b)
        @test c1 ≈ c2
    end
end

end


text:
# ---
# notebook: true
# ---
x = 1 + 2


text:
########################################################
#           ONE LAYER NEURONAL NETWORK TYPE
#   and evaluation with forward propagation
########################################################
# Constructores
export RandomWeightsNN
export FromMatrixNN
# Tipo
export AbstractOneLayerNeuralNetwork

"""
    AbstractOneLayerNeuralNetwork
The basic elements that define a one layer neural network
Must have two matrix:
W1: Matrix n x (d+1)
W2: Matris s x n
"""
abstract type AbstractOneLayerNeuralNetwork end

"""
    RandomWeightsNN
Return a random initialized Neuronal Network
"""
mutable struct RandomWeightsNN  <: AbstractOneLayerNeuralNetwork
    entry_dimesion :: Int
    number_of_hide_units :: Int
    output_dimension :: Int
    W1  # pesos de la entrada a la capa oculta A S (sesgo última columna)
    W2  # pesos de la capa oculta a la salida

    function RandomWeightsNN(entry_dimesion,
                                                number_of_hide_units,
                                                output_dimension)

        W1 = rand(Float64, number_of_hide_units, entry_dimesion+1)
        W2 = rand(Float64, output_dimension, number_of_hide_units)
        return new(
            entry_dimesion,
            number_of_hide_units,
            output_dimension,
            W1,
            W2
        )

    end
end

"""
    RandomWeightsNN
Return a  Neuronal Network inizialized  by three matrix
"""
mutable struct FromMatrixNN <: AbstractOneLayerNeuralNetwork
    W1 :: Matrix # pesos de la entrada a la capa oculta
    W2 :: Matrix# pesos de la capa oculta a la salida
    function FromMatrixNN(S,A,B)
        # Comprobación de que los tipos son correctos
        if !( typeof(S) <: Vector && typeof(A) <: Matrix && typeof(B) <: Matrix )
            throw(ArgumentError("El tipo de los argumentos no es el correcto\n
            Debería de ser:\n
             typeof(S) <: Vector && typeof(A) <: Matrix && typeof(B) <: Matrix \n
            pero se ha encontrado: \n
                typeof(S)= $(typeof(S))  typeof(A) $(typeof(A))  typeof(B) $(typeof(B))
            "))
        end
        # Comprobaciones de que los tamaños son coherentes
        (n_a,d_a) = size(A)
        l_s = length(S)
        (s_b, n_b) = size(B)
        ## Coherencia A y S
        if n_a != l_s
            throw(ArgumentError("El número de filas de A (que es $(n_a))debe de ser igual que
            la longitud de S (que es $(l_s))
            Los tamaños encontrados son:
                size(S)=$(size(S)).
                size(A)=$(size(A))
            "))
        end
        # Coherencia A y B
        if n_a != n_b
            throw(ArgumentError("El número de fila de A (que es $(n_a)) no es coherente con el número de columnas de B (que es $(n_b)).
            Ambos debería de ser iguales."))
        end
        return new(hcat(A,S), B)
    end
end

"""
    Base.show(io::IO, h <:AbstractOneLayerNeuralNetwork)
Implementamos el algoritmo de visualización de nuestras matrices
"""
function Base.show(io::IO, h ::AbstractOneLayerNeuralNetwork)
    display(Text("La matrix de pesos de las neuronas, W1, es:\n"))
    display(h.W1)
    display(Text("\nLa matrix de pesos de la salida, W2, es:\n"))
    display(h.W2)
end


text:
device(v::T) where T <: CuArray = gpu
device(v::T) where T <: AbstractArray = cpu
device(v::SubArray{T,N,P,I,L}) where {T, N, P <: CuArray, I, L} = gpu
device(v::SubArray{T,N,P,I,L}) where {T, N, P <: AbstractArray, I, L} = cpu
function device(c)
    p = Flux.params(c)
    length(p) > 0 && p[1] isa CuArray ? gpu : cpu
end

# Call F with input x but ensure they are both on the device of F
gpucall(F, x::CuArray) = F(x)
gpucall(F, x::SubArray{T,N,P,I,L}) where {T, N, P <: CuArray, I, L} = F(x)

gpucall(F, x::Array) = cpu(F(gpu(x)))
gpucall(F, x::SubArray{T,N,P,I,L}) where {T, N, P <: AbstractArray, I, L} = cpu(F(gpu(collect(x))))

cpucall(F, x::Array) = F(x)
cpucall(F, x::SubArray{T,N,P,I,L}) where {T, N, P <: AbstractArray, I, L} = F(x)

cpucall(F, x::CuArray) = gpu(F(cpu(x)))
cpucall(F, x::SubArray{T,N,P,I,L}) where {T, N, P <: CuArray, I, L} = gpu(F(cpu(x)))

mdcall(F, x, device) = device == gpu ? gpucall(F,x) : cpucall(F, x)

@inline function bslice(v, i)
    nd = ndims(v)
    if nd == 2
        return view(v,:,i)
    elseif nd == 3
        return view(v, :, :, i)
    elseif nd == 4
        return view(v, :, :, :, i)
    else
        return view(v, ntuple(x->:, nd-1)..., i)
    end
end


text:
using WaveSimulator
using BenchmarkTools

const SUITE = BenchmarkGroup()

SUITE["backend"] = BenchmarkGroup(["init", "allocation", "kernel"])
for (resource, tags) in (
        (CPU1(), ["CPU", "CPU1"]),
        (CPUThreads((100,100,1)), ["CPU", "CPUThreads"]),
        (CUDALibs(), ["GPU"]))
    GRP = SUITE["backend"][string(typeof(resource).name.name)] = BenchmarkGroup(tags)
    wave = UniformWave{3}(fmax=2e3)
    sim = Simulator(wave, resource=resource, duration=0.01)
    domain = BoxDomain(6,8,4, gamma=0.05)
    backend = WaveSimulator.backend_init(sim.resource, domain, sim)
    GRP["state_init_zeros"] = @benchmarkable(
        WaveSimulator.state_init($backend, $domain, $sim),
        teardown = gc(),
        seconds = 30,
        samples = 10
    )
    GRP["state_init_gaus"] = @benchmarkable(
        WaveSimulator.state_init(f0, $backend, $domain, $sim),
        setup = begin
            f0 = WaveSimulator.gauss($domain)
        end,
        teardown = gc(),
        seconds = 30,
        samples = 10
    )
    GRP["update!"] = @benchmarkable(
        WaveSimulator.backend_update!(state, backend, sim),
        setup = begin
            sim = Simulator($wave, resource=$(resource), duration=0.01)
            domain = BoxDomain(6,8,4, gamma=0.05)
            f0 = WaveSimulator.gauss(domain)
            backend = WaveSimulator.backend_init(sim.resource, domain, sim)
            state   = WaveSimulator.state_init(f0, backend, domain, sim)
        end,
        teardown = begin
            backend = nothing
            state = nothing
            gc()
        end,
        seconds = 30,
        samples = 50
    )
end

SUITE["integration"] = BenchmarkGroup(["simulate!"])
for (resource, tags) in (
        (CPU1(), ["CPU", "CPU1"]),
        (CPUThreads((100,100,1)), ["CPU", "CPUThreads"]),
        (CUDALibs(), ["GPU"]))
    GRP = SUITE["integration"][string(typeof(resource).name.name)] = BenchmarkGroup(tags)
    GRP["simulate!"] = @benchmarkable(
        simulate!(state, backend, sim),
        setup = begin
            wave = UniformWave{3}(fmax=2e3)
            sim = Simulator(wave, resource=$(resource), duration=0.01)
            domain = BoxDomain(6,8,4, gamma=0.05)
            f0 = WaveSimulator.gauss(domain)
            backend = WaveSimulator.backend_init(sim.resource, domain, sim)
            state   = WaveSimulator.state_init(f0, backend, domain, sim)
        end,
        teardown = begin
            backend = nothing
            state = nothing
            gc()
        end,
        seconds = 30,
        samples = 50
    )
end


text:
include("toml_consistency.jl")
include("param_boxes.jl")
include("types_from_file.jl")

include("old/planet.jl")
include("old/subgrid_scale.jl")
include("old/edmf.jl")
include("old/microphysics.jl")
include("old/atmos_subgrid_scale.jl")
include("old/surface_fluxes.jl")
include("old/override_defaults.jl")

max_stars_repo_name: climate-machine/CLIMAParameters

max_issues_repo_name: climate-machine/CLIMAParameters

max_forks_repo_name: climate-machine/CLIMAParameters


text:
# do nothing for other types
"""
    codegen_ast(def)

Generate Julia AST object `Expr` from a given syntax type.

# Example

One can generate the Julia AST object from a `JLKwStruct` syntax
type.

```julia
julia> def = @expr JLKwStruct struct Foo{N, T}
                  x::T = 1
              end
#= kw =# struct Foo{N, T}
    #= REPL[19]:2 =#
    x::T = 1
end

julia> codegen_ast(def)|>rm_lineinfo
quote
    struct Foo{N, T}
        x::T
    end
    begin
        function Foo{N, T}(; x = 1) where {N, T}
            Foo{N, T}(x)
        end
        function Foo{N}(; x::T = 1) where {N, T}
            Foo{N, T}(x)
        end
    end
end
```
"""
codegen_ast(ex) = ex

function codegen_ast(def::JLFor)
    lhead = Expr(:block)
    for (var, itr) in zip(def.vars, def.iterators)
        push!(lhead.args, :($var = $itr))
    end
    return Expr(:for, lhead, codegen_ast(def.kernel))
end

function codegen_ast(def::JLIfElse)
    isempty(def.conds) && return def.otherwise
    stmt = ex = Expr(:if)
    for (k, (cond, action)) in enumerate(def)
        push!(stmt.args, cond)
        push!(stmt.args, Expr(:block, codegen_ast(action)))

        if k !== length(def)
            push!(stmt.args, Expr(:elseif))
            stmt = stmt.args[end]
        end
    end
    def.otherwise === nothing || push!(stmt.args, codegen_ast(def.otherwise))
    return ex
end

function codegen_ast(fn::JLFunction)
    # handle anonymous syntax: function (x; kw=value) end
    if fn.head === :function && fn.name === nothing && fn.kwargs !== nothing &&
        isone(length(fn.args)) && isone(length(fn.kwargs))
        kw = fn.kwargs[1].args[1]
        va = fn.kwargs[1].args[2]
        call = Expr(:block, fn.args[1], :($kw = $va))
    else
        if fn.name === nothing
            call = Expr(:tuple)
        else
            call = Expr(:call, fn.name)
        end

        if fn.kwargs !== nothing
            push!(call.args, Expr(:parameters, fn.kwargs...))
        end
        append!(call.args, fn.args)
    end

    if fn.rettype !== nothing
        call = Expr(:(::), call, fn.rettype)
    end

    if fn.whereparams !== nothing
        call = Expr(:where, call, fn.whereparams...)
    end

    fn_def = Expr(fn.head, call, maybe_wrap_block(codegen_ast(fn.body)))
    return codegen_ast_docstring(fn, fn_def)
end

function codegen_ast(def::JLStruct)
    return codegen_ast_struct(def)
end

function codegen_ast(def::JLKwStruct)
    quote
        $(codegen_ast_struct(def))
        $(codegen_ast_kwfn(def))
    end
end

function maybe_wrap_block(ex::Expr)
    ex.head === :block && return ex
    return Expr(:block, ex)
end

"""
    codegen_ast_kwfn(def[, name = nothing])

Generate the keyword function from a Julia struct definition.

# Example

```julia
julia> def = @expr JLKwStruct struct Foo{N, T}
                  x::T = 1
              end
#= kw =# struct Foo{N, T}
    #= REPL[19]:2 =#
    x::T = 1
end

julia> codegen_ast_kwfn(def)|>prettify
quote
    function Foo{N, T}(; x = 1) where {N, T}
        Foo{N, T}(x)
    end
    function Foo{N}(; x::T = 1) where {N, T}
        Foo{N, T}(x)
    end
end

julia> def = @expr JLKwStruct struct Foo
                  x::Int = 1
              end
#= kw =# struct Foo
    #= REPL[23]:2 =#
    x::Int = 1
end

julia> codegen_ast_kwfn(def)|>prettify
quote
    function Foo(; x = 1)
        Foo(x)
    end
    nothing
end
```
"""
function codegen_ast_kwfn(def, name = nothing)
    quote
        $(codegen_ast_kwfn_plain(def, name))
        $(codegen_ast_kwfn_infer(def, name))
    end
end

"""
    codegen_ast_kwfn_plain(def[, name = nothing])

Generate the plain keyword function that does not infer type variables.
So that one can use the type conversions defined by constructors.
"""
function codegen_ast_kwfn_plain(def, name = nothing)
    isempty(def.fields) && return
    struct_name = struct_name_plain(def)

    if name === nothing # constructor method
        name = struct_name
        args = []
        whereparams = isempty(def.typevars) ? nothing : name_only.(def.typevars)
    else
        @gensym T
        args = [:(::Type{$T}), ]
        whereparams = [name_only.(def.typevars)..., :($T <: $struct_name)]
    end

    # do not generate kwfn if it's defined by the user
    has_kwfn_constructor(def, name) && return

    kwfn_def = JLFunction(;
        name = name,
        args = args,
        # NOTE:
        # do not use type annotations so that
        # we can use type conversion defined
        # by the constructors.
        kwargs = codegen_ast_fields(def.fields; just_name=true),
        whereparams = whereparams,
        body = Expr(:call, struct_name, [field.name for field in def.fields]...)
    )

    return codegen_ast(kwfn_def)
end

"""
    codegen_ast_kwfn_infer(def, name = nothing)

Generate the keyword function that infers the type.
"""
function codegen_ast_kwfn_infer(def, name = nothing)
    # no typevars to infer, use the plain one
    isempty(def.typevars) && return
    struct_name = struct_name_without_inferable(def)
    requires = uninferrable_typevars(def)
    # all uninferrable, use plain constructor
    length(requires) == length(def.typevars) && return

    if name === nothing # constructor method
        name = struct_name
        args = []
        whereparams = isempty(requires) ? nothing : requires
    else
        @gensym T
        ub = isempty(requires) ? def.name : Expr(:curly, def.name, requires...)
        args = [:(::Type{$T}), ]
        whereparams = [requires..., :($T <: $ub)]
    end

    # do not generate kwfn if it's defined by the user
    has_kwfn_constructor(def, name) && return

    kwfn_def = JLFunction(;
        name = name,
        args = args,
        # NOTE:
        # enable type annotations to infer typevars
        kwargs = codegen_ast_fields(def.fields; just_name=true),
        whereparams = whereparams,
        body = Expr(:call, struct_name,
            [field.name for field in def.fields]...)
    )
    return codegen_ast(kwfn_def)
end

"""
    codegen_ast_fields(fields; just_name::Bool=true)

Generate a list of Julia AST object for each field, only generate
a list of field names by default, option `just_name` can be turned
off to call [`codegen_ast`](@ref) on each field object.
"""
function codegen_ast_fields(fields; just_name::Bool=true)
    map(fields) do field
        name = just_name ? field.name : codegen_ast(field)
        support_default(field) || return name

        if field.default === no_default
            name
        else
            Expr(:kw, name, field.default)
        end
    end
end

"""
    struct_name_plain(def)

Plain constructor name. See also [`struct_name_without_inferable`](@ref).

# Example

```julia
julia> def = @expr JLKwStruct struct Foo{N, Inferable}
    x::Inferable = 1
end

julia> struct_name_plain(def)
:(Foo{N, Inferable})
```
"""
function struct_name_plain(def)
    isempty(def.typevars) && return def.name
    return Expr(:curly, def.name, name_only.(def.typevars)...)
end

"""
    struct_name_without_inferable(def; leading_inferable::Bool=true)

Constructor name that assume some of the type variables is inferred.
See also [`struct_name_plain`](@ref). The kwarg `leading_inferable`
can be used to configure whether to preserve the leading inferable
type variables, the default is `true` to be consistent with the
default julia constructors.

# Example

```julia
julia> def = @expr JLKwStruct struct Foo{N, Inferable}
    x::Inferable = 1
end

julia> struct_name_without_inferable(def)
:(Foo{N})

julia> def = @expr JLKwStruct struct Foo{Inferable, NotInferable}
    x::Inferable
end

julia> struct_name_without_inferable(def; leading_inferable=true)
:(Foo{Inferable, NotInferable})

julia> struct_name_without_inferable(def; leading_inferable=false)
:(Foo{NotInferable})
```
"""
function struct_name_without_inferable(def; leading_inferable::Bool=true)
    isempty(def.typevars) && return def.name
    required_typevars = uninferrable_typevars(def; leading_inferable=leading_inferable)
    isempty(required_typevars) && return def.name
    return Expr(:curly, def.name, required_typevars...)
end

function codegen_ast_docstring(def, body)
    def.doc === nothing && return body
    Expr(:macrocall, GlobalRef(Core, Symbol("@doc")), def.line, def.doc, body)
end

"""
    codegen_ast_struct_head(def)

Generate the struct head.

# Example

```julia
julia> using Expronicon

julia> def = JLStruct(:(struct Foo{T} end))
struct Foo{T}
end

julia> codegen_ast_struct_head(def)
:(Foo{T})

julia> def = JLStruct(:(struct Foo{T} <: AbstractArray end))
struct Foo{T} <: AbstractArray
end

julia> codegen_ast_struct_head(def)
:(Foo{T} <: AbstractArray)
```
"""
function codegen_ast_struct_head(def)
    head = def.name::Symbol
    if !isempty(def.typevars)
        head = Expr(:curly, head, def.typevars...)
    end

    if def.supertype !== nothing
        head = Expr(:(<:), head, def.supertype)
    end
    return head
end

"""
    codegen_ast_struct_body(def)

Generate the struct body.

# Example

```julia
julia> def = JLStruct(:(struct Foo
           x::Int

           Foo(x::Int) = new(x)
       end))
struct Foo
    x::Int
end

julia> codegen_ast_struct_body(def)
quote
    #= REPL[15]:2 =#
    x::Int
    Foo(x::Int) = begin
            #= REPL[15]:4 =#
            new(x)
        end
end
```
"""
function codegen_ast_struct_body(def)
    body = Expr(:block)
    for field in def.fields
        field.line === nothing || push!(body.args, field.line)
        field.doc === nothing || push!(body.args, field.doc)
        push!(body.args, codegen_ast(field))
    end

    for constructor in def.constructors
        push!(body.args, codegen_ast(constructor))
    end

    body = flatten_blocks(body)
    def.misc === nothing || append!(body.args, def.misc)
    return body
end

"""
    codegen_ast_struct(def)

Generate pure Julia struct `Expr` from struct definition. This is equivalent
to `codegen_ast` for `JLStruct`. See also [`codegen_ast`](@ref).

# Example

```julia-repl
julia> def = JLKwStruct(:(struct Foo
           x::Int=1

           Foo(x::Int) = new(x)
       end))
struct Foo
    x::Int = 1
end

julia> codegen_ast_struct(def)
:(struct Foo
      #= REPL[21]:2 =#
      x::Int
      Foo(x::Int) = begin
              #= REPL[21]:4 =#
              new(x)
          end
  end)
```
"""
function codegen_ast_struct(def)
    head = codegen_ast_struct_head(def)
    body = codegen_ast_struct_body(def)
    ex = Expr(:struct, def.ismutable, head, body)
    return codegen_ast_docstring(def, ex)
end

function codegen_ast(def::Union{JLField, JLKwField})
    return if def.type === Any
        def.name
    else
        :($(def.name)::$(def.type))
    end
end

# X functions, the x<name> functions from Zygote, IRTools etc.

"""
    xtuple(xs...)

Create a `Tuple` expression.
"""
xtuple(xs...) = Expr(:tuple, xs...)

"""
    xnamedtuple(;kw...)

Create a `NamedTuple` expression.
"""
function xnamedtuple(;kw...)
    ex = Expr(:tuple)
    for (k, v) in kw
        push!(ex.args, :($k = $v))
    end
    return ex
end

"""
    xcall(name, args...; kw...)

Create a function call to `name`.
"""
function xcall(name, args...; kw...)
    isempty(kw) && return Expr(:call, name, args...)
    p = Expr(:parameters)
    for (k, v) in kw
        push!(p.args, Expr(:kw, k, v))
    end
    Expr(:call, name, p, args...)
end

"""
    xcall(m::Module, name::Symbol, args...; kw...)

Create a function call to `GlobalRef(m, name)`.

!!! tip

    due to [Revise/#616](https://github.com/timholy/Revise.jl/issues/616),
    to make your macro work with Revise, we use the dot expression
    `Expr(:., <module>, QuoteNode(<name>))` instead of `GlobalRef` here.
"""
function xcall(m::Module, name::Symbol, args...; kw...)
    # NOTE: not use GlobalRef due to Revise#616
    xcall(Expr(:., m, QuoteNode(name)), args...; kw...)
end

"""
    xgetindex(collection, key...)

Create a function call expression to `Base.getindex`.
"""
xgetindex(coll, key...) = xcall(Base, :getindex, coll, key...)

"""
    xpush(collection, items...)

Create a function call expression to `Base.push!`.
"""
function xpush(collection, items...)
    xcall(Base, :(push!), collection, items...)
end

"""
    xfirst(collection)

Create a function call expression to `Base.first`.
"""
xfirst(collection) = xcall(Base, :first, collection)

"""
    xlast(collection)

Create a function call expression to `Base.last`.
"""
xlast(collection) = xcall(Base, :last, collection)

"""
    xprint(xs...)

Create a function call expression to `Base.print`.
"""
xprint(xs...) = xcall(Base, :print, xs...)

"""
    xprintln(xs...)

Create a function call expression to `Base.println`.
"""
xprintln(xs...) = xcall(Base, :println, xs...)

"""
    xmap(f, xs...)

Create a function call expression to `Base.map`.
"""
xmap(f, xs...) = xcall(Base, :map, f, xs...)

"""
    xmapreduce(f, op, xs...)

Create a function call expression to `Base.mapreduce`.
"""
xmapreduce(f, op, xs...) = xcall(Base, :mapreduce, f, op, xs...)

"""
    xiterate(it)

Create a function call expression to `Base.iterate`.
"""
xiterate(it) = xcall(Base, :iterate, it)

"""
    xiterate(it, st)

Create a function call expression to `Base.iterate`.
"""
xiterate(it, st) = xcall(Base, :iterate, it, st)


text:
"""
Inverse Fourier transform in zonal direction.
"""
function fourier_inverse(   input::Array{Complex{NF},2},
                            G::GeoSpectral{NF},
                            scale::Bool=false) where {NF<:AbstractFloat}

    @unpack nlon, nlat, nlon_half, cosg⁻¹ = G.geometry
    @unpack mx = G.spectral

    # preallocate, TODO turn into plan_irfft for performance
    output = zeros(NF, nlon, nlat)
    nlonNF = convert(NF,nlon)

    for j in 1:nlat
        # Do inverse FFT then multiply by number of longitudes
        # add the truncated wavenumbers with zeros
        output[:,j] = nlonNF*irfft(vcat(input[:,j], zeros(Complex{NF}, nlon_half+1-mx)), nlon)
    end

    # Scale by cosine(lat) if needed
    if scale
        for j in 1:nlat
            output[:,j] *= cosg⁻¹[j]
        end
    end

    return output
end

"""
Fourier transform in the zonal direction.
"""
function fourier(   input::Array{NF,2},
                    G::GeoSpectral{NF}) where {NF<:AbstractFloat}

    @unpack nlon, nlat = G.geometry
    @unpack mx = G.spectral

    # preallocate output
    #TODO pass on output array as argument, turn fourier into fourier!
    output = zeros(Complex{NF}, mx, nlat)

    one_over_nlon = convert(NF,1/nlon)

    # Copy grid-point data into working array
    for j in 1:nlat
        # Do direct FFT then divide by number of longitudes
        #TODO use plan_rfft and use a type-flexible FFT instead of FFTW
        output[:,j] = one_over_nlon*rfft(input[:,j])[1:mx]  # truncate to mx ?
    end

    return output
end


text:
## utility functions ##

function reorder(values::AbstractVector{T}) where T<:Real
    newOrder = sortperm(values)
    oldOrder = sortperm(newOrder)
    return newOrder, oldOrder
end

function sort_if_needed(x; kws...)
    if issorted(x; kws...)
        return x
    else
        return sort(x; kws...)
    end
end

function sort_if_needed!(x; kws...)
    if !issorted(x; kws...)
        sort!(x; kws...)
    end
end

function unsort(x; kws...)
    y = copy(x)
    while issorted(y; kws...)
        shuffle!(y)
    end
    return y
end

function unorder(x; kws...)
    ord = collect(1:length(x))
    while issorted(x[ord]; kws...)
        shuffle!(ord)
    end
    return ord
end

function valid_pvalues(x::AbstractVector{T}) where T<:AbstractFloat
    if !isin(x)
        throw(DomainError("p-values must all be in [0, 1]"))
    end
end

function isin(x::Real, lower::Real = 0., upper::Real = 1.)
    x >= lower && x <= upper
end

function isin(x::AbstractVector{T}, lower::Real = 0., upper::Real = 1.) where T<:Real
    ex = extrema(x)
    ex[1] >= lower && ex[2] <= upper
end

function isotonic_regression_reference(y::AbstractVector{T}, w::AbstractVector{T}) where T<:AbstractFloat
    # TODO ignore zero weights
    y = copy(y)
    w = copy(w)
    cnts = fill(1, length(y))
    i = 2
    while i <= length(y)
        if y[i] < y[i-1]
            y[i-1] = (w[i]*y[i]+w[i-1]*y[i-1])/(w[i]+w[i-1])
            w[i-1] = w[i]+w[i-1]
            cnts[i-1] += cnts[i]
            deleteat!(y, i)
            deleteat!(w, i)
            deleteat!(cnts, i)
            i = max(i-2, 1)
        end
        i += 1
    end
    yisotonic = vcat([y[idx] .* fill(1.0, cnt) for (idx, cnt) in enumerate(cnts)]...)
    return yisotonic
end

function isotonic_regression_reference(y::AbstractVector{T}) where T<:AbstractFloat
    isotonic_regression_reference(y, fill(1.0, size(y)))
end

function isotonic_regression(y::AbstractVector{T}, weights::AbstractVector{T}) where T<:AbstractFloat
    n = length(y)
    if n <= 1
        return y
    end
    if n != length(weights)
        throw(DimensionMismatch("Lengths of values and weights mismatch"))
    end
    @inbounds begin
        n -= 1
        while true
            i = 1
            is_pooled = false
            while i <= n
                k = i
                while k <= n && y[k] >= y[k+1]
                    k += 1
                end
                if y[i] != y[k]
                    numerator = 0.0
                    denominator = 0.0
                    for j in i:k
                        numerator += y[j] * weights[j]
                        denominator += weights[j]
                    end
                    m = numerator / denominator
                    for j in i:k
                        y[j] = m
                    end
                    is_pooled = true
                end
                i = k + 1
            end
            if !is_pooled
               break
            end
        end
    end
    return y
end

function isotonic_regression(y::AbstractVector{T}) where T<:AbstractFloat
    isotonic_regression(y, fill(1.0, size(y)))
end

function grenander(pv::AbstractVector{T}) where T<:AbstractFloat
    pv_sorted = sort(pv)
    # ecdf that handles duplicated values
    pv_sorted_unique, counts = rle(pv_sorted)
    ecdf_value = cumsum(counts)
    ecdf_value = ecdf_value ./ ecdf_value[end]

    Δx = diff(pv_sorted_unique)
    Δy = diff(ecdf_value)

    f = Δy ./ Δx
    f = -isotonic_regression(-f, Δx)
    F  = ecdf_value[1] .+ vcat(0, cumsum(f .* Δx))
    f = push!(f, f[end])

    return pv_sorted_unique, f, F
end


text:
# Copyright © 2021 Alexander L. Hayes
# Apache 2.0 License

"""Request copies of relational datasets.
"""

using HTTP
using ZipFile

include("base.jl")
include("types.jl")

"""
Available datasets from the
[`srlearn/datasets`](https://github.com/srlearn/datasets) repository.
"""
DATASETS = [
    "toy_cancer",
    "toy_father",
    "citeseer",
    "cora",
    "uwcse",
    "webkb",
    "financial_nlp_small",
    "nell_sports",
    "icml",
    "boston_housing",
    "drug_interactions",
    "toy_machines",
]

"""Default download version.

If a "version" parameter is not passed to `load`, a dataset of this version
is downloaded by default.
"""
LATEST_VERSION = "v0.0.5"

"""
    load(name::String, version::Union{String, Nothing} = nothing; fold::Int64 = 1)

Load the training and test folds for a dataset.
"""
function load(name::String, version::Union{String, Nothing} = nothing; fold::Int64 = 1)
    data_location = fetch(name, version)
    return deserialize_zipfile(data_location, name, fold=fold)
end

function deserialize_zipfile(data_location::String, name::String; fold::Int64 = 1)

    reader = ZipFile.Reader(data_location)

    folds = _n_folds(reader)

    train_pos, train_neg, train_facts = nothing, nothing, nothing
    test_pos, test_neg, test_facts = nothing, nothing, nothing

    # TODO(hayesall): The ZipFile.jl doesn't make it quite as easy to seek a file by name.
    #   There's a package for building iterators out of tar archives, but it appears that
    #   the only way to do a "file name lookup" is by building an intermediate dictionary
    #   and jumping to the index of the name you want.
    file_num = Dict{String, Int64}()
    for (i, fname) in enumerate(reader.files)
        file_num[fname.name] = i
    end

    if folds == 0
        train_pos = readlines(reader.files[file_num["$(name)/train/train_pos.txt"]])
        train_neg = readlines(reader.files[file_num["$(name)/train/train_neg.txt"]])
        train_facts = readlines(reader.files[file_num["$(name)/train/train_facts.txt"]])

        test_pos = readlines(reader.files[file_num["$(name)/test/test_pos.txt"]])
        test_neg = readlines(reader.files[file_num["$(name)/test/test_neg.txt"]])
        test_facts = readlines(reader.files[file_num["$(name)/test/test_facts.txt"]])

    elseif fold > folds
        throw(error("Fold ($(fold)) does not exist in $(data_location)"))

    else

        train_pos = readlines(reader.files[file_num["$(name)/fold$(fold)/train/train_pos.txt"]])
        train_neg = readlines(reader.files[file_num["$(name)/fold$(fold)/train/train_neg.txt"]])
        train_facts = readlines(reader.files[file_num["$(name)/fold$(fold)/train/train_facts.txt"]])

        test_pos = readlines(reader.files[file_num["$(name)/fold$(fold)/test/test_pos.txt"]])
        test_neg = readlines(reader.files[file_num["$(name)/fold$(fold)/test/test_neg.txt"]])
        test_facts = readlines(reader.files[file_num["$(name)/fold$(fold)/test/test_facts.txt"]])

    end

    close(reader)

    return RelationalDataset((
        train_pos,
        train_neg,
        train_facts,
    )), RelationalDataset((
        test_pos,
        test_neg,
        test_facts,
    ))

end

"""
    fetch(name::String, version::Union{String, Nothing} = nothing)

Download a dataset with `name` and `version`, write the zipfile to the cache
directory, and return a string representing the path to the file on disk.
"""
function fetch(name::String, version::Union{String, Nothing} = nothing)::String

    data_file = _make_file_path(name, version)
    if Base.Filesystem.isfile(data_file)
        return data_file
    end

    download_url = _make_data_url(name, version)
    request = HTTP.request("GET", download_url)

    io = open(data_file, "w")
    write(io, request.body)
    close(io)

    return data_file
end

function _make_file_path(name::String, version::Union{String, Nothing} = nothing)::String
    @assert name in DATASETS
    if version == nothing
        return Base.Filesystem.joinpath(
            get_data_home(),
            "$(name)_$(LATEST_VERSION).zip",
        )
    end
    return Base.Filesystem.joinpath(
        get_data_home(),
        "$(name)_$(version).zip",
    )
end

function _make_data_url(name::String, version::Union{String, Nothing} = nothing)::String
    @assert name in DATASETS
    if version == nothing
        version = LATEST_VERSION
    end

    # TODO(hayesall): Enforce no v0.0.2 or v0.0.3
    return "https://github.com/srlearn/datasets/releases/download/$(version)/$(name)_$(version).zip"
end

"""Does this zipfile contain CV folds?"""
function _has_folds(zip::ZipFile.Reader)::Bool
    return any([occursin("fold", file.name) for file in zip.files])
end

"""How many folds does this contain?"""
function _n_folds(zip::ZipFile.Reader)::Int64

    if !_has_folds(zip)
        return 0
    end

    numbers = Int64[]
    for path in zip.files
        if occursin("fold", path.name)
            append!(numbers, parse(Int64, split(split(path.name, "fold")[2], "/")[1]))
        end
    end

    return maximum(numbers)
end


text:
##Differentiation and integration

Base.sum(f::Fun{Laurent{DD,RR}}) where {DD<:PeriodicSegment,RR} = coefficient(f,1).*arclength(domain(f))
Base.sum(f::Fun{Laurent{DD,RR}}) where {DD<:Circle,RR} = coefficient(f,3 - domain(f).orientation).*complexlength(domain(f))

Base.sum(f::Fun{Fourier{DD,RR}}) where {DD<:PeriodicSegment,RR} = coefficient(f,1).*arclength(domain(f))
Base.sum(f::Fun{Fourier{DD,RR}}) where {DD<:Circle,RR} =
    (im*coefficient(f,2) + coefficient(f,3))/2*complexlength(domain(f))

function linesum(f::Fun{Laurent{DD,RR}}) where {DD,RR}
    d=domain(f)
    if isa(d,Circle)
        sum(setcanonicaldomain(f))*d.radius
    else
        sum(f)
    end
end

linesum(f::Fun{Fourier{DD,RR}}) where {DD<:Circle,RR} = sum(setcanonicaldomain(f))*domain(f).radius
linesum(f::Fun{Fourier{DD,RR}}) where {DD<:PeriodicSegment,RR} = sum(f) #TODO: Complex periodic interval

differentiate(f::Fun{Taylor{DD,RR}}) where {DD<:PeriodicSegment,RR} =
    Fun(f.space,im*tocanonicalD(f,0)*taylor_diff(f.coefficients))
differentiate(f::Fun{Hardy{false,DD,RR}}) where {DD<:PeriodicSegment,RR} =
    Fun(f.space,im*tocanonicalD(f,0)*hardyfalse_diff(f.coefficients))
differentiate(f::Fun{Laurent{DD,RR}}) where {DD<:PeriodicSegment,RR} =
    Fun(f.space,im*tocanonicalD(f,0)*laurentdiff(f.coefficients))

differentiate(f::Fun{CosSpace{DD,RR}}) where {DD<:PeriodicSegment,RR} =
    Fun(SinSpace(domain(f)),tocanonicalD(f,0)*cosspacediff(f.coefficients))
differentiate(f::Fun{SinSpace{DD,RR}}) where {DD<:PeriodicSegment,RR} =
    Fun(CosSpace(domain(f)),tocanonicalD(f,0)*sinspacediff(f.coefficients))
differentiate(f::Fun{Fourier{DD,RR}}) where {DD<:PeriodicSegment,RR} =
    Fun(f.space,tocanonicalD(f,0)*fourierdiff(f.coefficients))

differentiate(f::Fun{Laurent{DD,RR}}) where {DD,RR} = Derivative(space(f))*f
differentiate(f::Fun{Fourier{DD,RR}}) where {DD,RR} = Derivative(space(f))*f

function integrate(f::Fun{Hardy{false,D,R}}) where {D,R}
    if isa(domain(f),Circle) # drop -1 term if zero and try again
        @assert ncoefficients(f)==0 || abs(f.coefficients[1])<100eps()
        integrate(Fun(f,space(f)|(2:∞)))
    else  # Probably periodic itnerval
        Integral(space(f))*f
    end
end

function integrate(f::Fun{Taylor{D,R}}) where {D,R}
    if isa(domain(f),Circle)
        Integral(space(f))*f
    else  # Probably periodic itnerval  drop constant term if zero
        @assert ncoefficients(f)==0 || abs(f.coefficients[1])<100eps()
        Fun(integrate(Fun(f,space(f)|(2:∞))),space(f))
    end
end

Base.sum(f::Fun{CosSpace{DD,RR}}) where {DD<:PeriodicSegment,RR} =
    f.coefficients[1]*complexlength(domain(f))

linesum(f::Fun{CosSpace{DD,RR}}) where {DD<:PeriodicSegment,RR} =
    f.coefficients[1]*arclength(domain(f))

function integrate(f::Fun{CS}) where CS<:CosSpace
    if isa(domain(f),Circle)
        error("Integrate not implemented for CosSpace on Circle")
    else  # Probably periodic itnerval, drop constant term if zero
        tol=1E-14 #TODO: smart tolerance.  Here relative is a bit tricky
                  # since this is called by Fourier integrate
        if abs(f.coefficients[1])<tol
            integrate(Fun(f,space(f)|(2:∞)))
        else
            d=domain(f)
            @assert isa(d,PeriodicSegment)
            x=Fun(identity, Interval(d))
            (f.coefficients[1]*x)⊕integrate(Fun(f,space(f)|(2:∞)))
        end
    end
end

function integrate(f::Fun{SS}) where SS<:SinSpace
    if isa(domain(f),Circle) # drop term containing z^(-1)
        integrate(Fun(f,space(f)|(2:∞)))
    else  # Probably periodic itnerval\
        Integral(space(f))*f
    end
end

#TODO: This is a hack to make sure Fourier maps to Fourier
# we don't have banded differentiate from CosSpace/SinSpace on a circle
for OP in (:differentiate,:integrate)
    @eval $OP(f::Fun{Fourier{D,R},T}) where {T,D<:Circle,R} = $OP(Fun(f,Laurent))
end

integrate(f::Fun{Fourier{D,R},T}) where {T,D<:PeriodicSegment,R} =
    integrate(component(f,2))⊕integrate(component(f,1))

# O(min(m,n)) Laurent line integral

function linebilinearform(f::Fun{Laurent{D,R},T},g::Fun{Laurent{D,R},T}) where {T,D<:Circle,R}
    @assert domain(f) == domain(g)
    u,v,mn = f.coefficients,g.coefficients,min(ncoefficients(f),ncoefficients(g))
    if mn > 1
        ret = u[1]*v[1]
        for i=2:2:mn-1
            ret += u[i]*v[i+1] + u[i+1]*v[i]
        end
        return arclength(domain(f))*ret
    elseif mn > 0
        return arclength(domain(f))*u[1]*v[1]
    else
        return zero(T)
    end
end

function bilinearform(f::Fun{Laurent{D,R},T},g::Fun{Laurent{D,R},T}) where {T,D<:Circle,R}
    @assert domain(f) == domain(g)
    u,v,mn = f.coefficients,g.coefficients,min(ncoefficients(f),ncoefficients(g))
    if mn > 2
        ret = u[1]*v[2] + u[2]*v[1]
        for i=3:2:mn-1
            ret += u[i]*v[i+1] + u[i+1]*v[i]
        end
        return complexlength(domain(f))*ret
    else
        return zero(T)
    end
end


text:
export QdldlKKTSolver, CholmodKKTSolver
# -------------------------------------
# abstract type defs
# -------------------------------------
abstract type AbstractKKTSolver end
abstract type AbstractKKTShape end

# NB: all concrete examples of this type should
# implement refactor! and solve! methods and should
# implement a constructor taking (P,A,sigma,rho) as
# arguments

# -------------------------------------
# some internal utility functions
# -------------------------------------

"Check the dimensions of the problem data used to assemble the KKT matrix `K`."
function _kktutils_check_dims(P::SparseMatrixCSC{Tv, Ti}, A::AbstractMatrix{Tv}, sigma::Tv, rho::Union{Tv, AbstractVector{Tv}}) where {Tv <: AbstractFloat, Ti <: Integer}

    n = size(P, 1)
    m = size(A, 1)

    size(A,2) == n || throw(DimensionMismatch())

    length(rho)   == m || length(rho)   == 1 || throw(DimensionMismatch())
    length(sigma) == n || length(sigma) == 1 || throw(DimensionMismatch())

    return m, n
end

"Count the number of nnz in each column of the lower triangle of the KKT matrix `K`."
function _count_lower_triangle!(Kcolnz::Vector{Ti}, P::SparseMatrixCSC{Tv, Ti}, A::SparseMatrixCSC{Tv, Ti}, n::Ti) where {Tv <: AbstractFloat, Ti <: Integer}
    # count strict lower triangular values of P
    @inbounds for cidx = 1:n
        for j = (P.colptr[cidx]):(P.colptr[cidx+1]-1)
            ridx = P.rowval[j]
            if (ridx > cidx)
                Kcolnz[cidx] += 1
            end
        end
    end
    # add number of column entries of A
    @inbounds for cidx = 1:n
        Kcolnz[cidx] += A.colptr[cidx+1] - A.colptr[cidx]
    end

    #every element on the diagonal is also nonzero
    Kcolnz .+= 1
end

"Count the number of nnz in each column of the upper triangle of the KKT matrix `K`."
function _count_upper_triangle!(Kcolnz::Vector{Ti}, P::SparseMatrixCSC{Tv, Ti}, A::SparseMatrixCSC{Tv, Ti}, n::Ti) where {Tv <: AbstractFloat, Ti <: Integer}
    # count nonzeros in the strict upper triangle of (P)
    @inbounds for cidx = 1:n
        for j = (P.colptr[cidx]):(P.colptr[cidx+1]-1)
            ridx = P.rowval[j]
            if (ridx < cidx)
                Kcolnz[cidx] += 1
            else
                break
            end
        end
    end

    #count the nonzeros in columns in A'
    @inbounds for idx in A.rowval
        Kcolnz[idx + n] += 1
    end

    #every element on the diagonal is also nonzero
    Kcolnz .+= 1
end

"Determine and fill in the rowvals `Ki` and nzvals `Kx` of the lower triangle of `K`."
function _fill_lower_triangle!(Ki::Vector{Ti}, Kp::Vector{Ti}, Kx::Vector{Tv}, P::SparseMatrixCSC{Tv, Ti}, A::SparseMatrixCSC{Tv, Ti}, R::Union{Tv, AbstractVector{Tv}}, S::Union{Tv, AbstractVector{Tv}}, KNextInCol::Vector{Ti}, m::Ti, n::Ti) where {Tv <: AbstractFloat, Ti <: Integer}
    # tril(K) = | tril(P) + σI |         |
    #           |   A          | - 1/ρ I |

    #fill in the Sigma part
    @inbounds for cidx = 1:n
            Kx[KNextInCol[cidx]]      = S[cidx]
            Ki[KNextInCol[cidx]]      = cidx
            KNextInCol[cidx] += 1
    end

    #fill in tril(P)
    @inbounds for cidx = 1:n
        for j = (P.colptr[cidx]):(P.colptr[cidx+1]-1)
            ridx = P.rowval[j]
            if ridx > cidx
                Kidx = KNextInCol[cidx]
                Ki[Kidx] = ridx
                Kx[Kidx] = P.nzval[j]
                KNextInCol[cidx] += 1
            elseif ridx == cidx # add diagonal entries of P onto the sigma value
                Kidx = KNextInCol[cidx] - 1
                Kx[Kidx] += P.nzval[j]
            end
        end
    end

    #fill in A in the lower LHC
    @inbounds for cidx = 1:n
        for j = (A.colptr[cidx]):(A.colptr[cidx+1]-1)
            ridx = A.rowval[j]
            ridxpn = ridx + n
            Kidx = KNextInCol[cidx]
            Ki[Kidx] = ridxpn
            Kx[Kidx] = A.nzval[j]
            KNextInCol[cidx] += 1
        end
    end

    #fill in the rho part
    @inbounds for idx = 1:m
        kidx = KNextInCol[idx + n]
        Kx[kidx] = R[idx]
        Ki[kidx] = idx + n
    end
end

"Determine and fill in the rowvals `Ki` and nzvals `Kx` of the upper triangle of `K`."
function _fill_upper_triangle!(Ki::Vector{Ti}, Kp::Vector{Ti}, Kx::Vector{Tv}, P::SparseMatrixCSC{Tv, Ti}, A::SparseMatrixCSC{Tv, Ti}, R::Union{Tv, AbstractVector{Tv}}, S::Union{Tv, AbstractVector{Tv}}, KNextInCol::Vector{Ti}, m::Ti, n::Ti) where {Tv <: AbstractFloat, Ti <: Integer}
    # triu(K) = | triu(P) + σI |    A'   |
    #           |              | - 1/ρ I |

    #fill in triu(P)
    @inbounds for cidx = 1:n
        for j = (P.colptr[cidx]):(P.colptr[cidx+1]-1)
            ridx = P.rowval[j]
            if ridx <= cidx
                Kidx = KNextInCol[cidx]
                Ki[Kidx] = ridx
                Kx[Kidx] = P.nzval[j]
                KNextInCol[cidx] += 1
            else
                break
            end
        end
    end

    #fill in A' in the upper RHS
    @inbounds for cidx = 1:n
        for j = (A.colptr[cidx]):(A.colptr[cidx+1]-1)
            ridx = A.rowval[j]
            ridxpn = ridx + n
            Kidx = KNextInCol[ridxpn]
            KNextInCol[ridxpn] += 1
            Ki[Kidx] = cidx
            Kx[Kidx] = A.nzval[j]
        end
    end

    #fill in the Sigma part
    @inbounds for cidx = 1:n
        if KNextInCol[cidx] == Kp[cidx+1]       #P had diagonal term at cidx)
            Kx[KNextInCol[cidx] - 1] += S[cidx]
        else
            Kx[KNextInCol[cidx]]      = S[cidx]
            Ki[KNextInCol[cidx]]      = cidx
        end
    end
    #fill in the rho part
    @inbounds for idx = 1:m
        kidx = KNextInCol[idx+n]
        Kx[kidx] = R[idx]
        Ki[kidx] = idx+n
    end
end

# handle the different forms of inputs for sigma and rho: here both are vectors (*)
function assemble_kkt_triangle(P::SparseMatrixCSC{Tv, Ti}, A::SparseMatrixCSC{Tv, Ti}, sigma::AbstractVector{Tv}, rho::AbstractVector{Tv}, shape::Symbol = :U) where {Tv <: AbstractFloat, Ti <: Integer}
        S = sigma
        R  = Array{Tv}(undef, size(A, 1))
     @. R  = -one(Tv) / rho
     # pass on to main function
    return _assemble_kkt_triangle(P, A, S, R, shape)
end
# only sigma is a scalar
function assemble_kkt_triangle(P::SparseMatrixCSC{Tv, Ti}, A::SparseMatrixCSC{Tv, Ti}, sigma::Tv, rho::AbstractVector{Tv}, shape::Symbol = :U) where {Tv <: AbstractFloat, Ti <: Integer}
    S = Vector{Tv}(undef, size(P, 1))
    fill!(S, sigma)
    # put again through (*)
    return assemble_kkt_triangle(P, A, S, rho, shape)
end

# only rho is a scalar
function assemble_kkt_triangle(P::SparseMatrixCSC{Tv, Ti}, A::SparseMatrixCSC{Tv, Ti}, sigma::AbstractVector{Tv}, rho::Tv, shape::Symbol = :U) where {Tv <: AbstractFloat, Ti <: Integer}
    R  = Array{Tv}(undef, size(A, 1))
    fill!(R, -one(Tv) ./ rho)
    # pass on to main function
    return _assemble_kkt_triangle(P, A, sigma, R, shape)
end

# rho and sigma are scalars
function assemble_kkt_triangle(P::SparseMatrixCSC{Tv, Ti}, A::SparseMatrixCSC{Tv, Ti}, sigma::Tv, rho::Tv, shape::Symbol = :U) where {Tv <: AbstractFloat, Ti <: Integer}
    S = Vector{Tv}(undef, size(P, 1))
    fill!(S, sigma)

    R  = Array{Tv}(undef, size(A, 1))
    fill!(R, -one(Tv) ./ rho)
    # pass on to main function
    return _assemble_kkt_triangle(P, A, S, R, shape)
end
"Given `P`, `A`, `sigma` and `rho` return the upper / lower triangle, defined by `shape`,  of the KKT condition matrix `K`."
function _assemble_kkt_triangle(P::SparseMatrixCSC{Tv, Ti}, A::SparseMatrixCSC{Tv, Ti}, S::AbstractVector{Tv}, R::AbstractVector{Tv}, shape::Symbol = :U) where {Tv <: AbstractFloat, Ti <: Integer}

    #   K = | P + σI |  A'     |
    #       |   A    | - 1/ρ I |
    #         left      right
    n = size(P, 1)
    m = size(A, 1)

    # make the column counter of K
    Kcolnz = zeros(Ti, m + n)

    #work out how many nonzeros are in K.
    if shape == :U
        _count_upper_triangle!(Kcolnz, P, A, n)
    elseif shape == :L
        _count_lower_triangle!(Kcolnz, P, A, n)
    end

    #make the column pointer and nonzero count
    Kp   = Vector{Ti}(undef, m + n + 1)
    Kp[1] = 1
    @inbounds for i = 1:(m+n)
        Kp[i+1] = Kp[i] + Kcolnz[i]
    end

    KNextInCol = copy(Kp);        #next value in column goes here
    nnzK = Kp[end] - 1
    Ki   = Vector{Ti}(undef, nnzK)
    Kx   = Vector{Tv}(undef, nnzK)

    # fill in the rowvals and nzvals
    if shape == :U
        _fill_upper_triangle!(Ki, Kp, Kx, P, A, R, S, KNextInCol, m, n)
    elseif shape == :L
        _fill_lower_triangle!(Ki, Kp, Kx, P, A, R, S, KNextInCol, m, n)
    end

    #assemble the matrix
    K = SparseMatrixCSC(m + n, m + n, Kp, Ki, Kx)

    return K
end

"Given `P`, `A`, `sigma` and `rho` return the full KKT condition matrix `K`."
function _assemble_kkt_full(P::SparseMatrixCSC{Tv, Ti}, A::AbstractMatrix{Tv}, sigma::Tv, rho::Union{Tv, AbstractVector{Tv}}) where {Tv <: AbstractFloat, Ti <: Integer}

    n = size(P, 1)
    m = size(A, 1)
    S = length(sigma) == 1 ? (sigma[1]) * I : Diagonal(sigma)
    rhoinv  = Vector{Tv}(undef, m)
    rhoinv .= (-one(Tv)./rho)
    D       = SparseMatrixCSC(Diagonal(rhoinv))

    #compute the full KKT matrix
    K = [P + S SparseMatrixCSC(A'); A D]

    return K
end

"Update the diagonal `rho` entries in the lower right hand corner of the KKT matrix `K`."
function update_kkt_matrix!(K::SparseMatrixCSC{Tv, Ti}, n::Ti, m::Ti, rho::Tv) where {Tv <: AbstractFloat, Ti <: Integer}
    @inbounds @simd for i = (n + 1):(n + m)
        K[i, i] = -one(Tv) / rho
    end
end

function update_kkt_matrix!(K::SparseMatrixCSC{Tv, Ti}, n::Ti, m::Ti, rho_vec::AbstractVector{Tv}) where {Tv <: AbstractFloat, Ti <: Integer}
    @inbounds @simd for i = (n + 1):(n + m)
        K[i, i] = -one(Tv) / rho_vec[i - n]
    end
end

# -------------------------------------
# QDLDL solver
# -------------------------------------

struct QdldlKKTSolver{Tv, Ti} <: AbstractKKTSolver

    m::Ti
    n::Ti
    ldlfact::QDLDL.QDLDLFactorisation{Tv, Ti}

    function QdldlKKTSolver(P::SparseMatrixCSC{Tv, Ti}, A::SparseMatrixCSC{Tv, Ti}, sigma::Tv, rho::Union{Tv, AbstractVector{Tv}}) where {Tv <: AbstractFloat, Ti <: Integer}
        m, n = _kktutils_check_dims(P, A, sigma, rho)
        K    = assemble_kkt_triangle(P, A, sigma, rho, :U)
        ldlfact = qdldl(K)

        #check for exactly n positive eigenvalues
        positive_inertia(ldlfact) == n || error("Objective function is not convex.")

        new{Tv, Ti}(m, n, ldlfact)
    end
end

function solve!(s::QdldlKKTSolver, lhs, rhs)
    @. lhs = rhs
    QDLDL.solve!(s.ldlfact, lhs)
end

function update_rho!(s::QdldlKKTSolver{Tv, Ti}, rho::Union{Tv, AbstractVector{Tv}}) where {Tv <: AbstractFloat, Ti <: Integer}
    QDLDL.update_diagonal!(s.ldlfact, (s.n+1):(s.n+s.m), (-one(Ti) ./ rho))
end

# -------------------------------------
# Julia Native solver (CHOLMOD based)
# -------------------------------------
mutable struct CholmodKKTSolver{Tv, Ti} <: AbstractKKTSolver

    fact::SuiteSparse.CHOLMOD.Factor{Float64} # SuiteSparse's cholmod only supports doubles
    K::SparseMatrixCSC{Tv, Ti}
    m::Ti
    n::Ti

    function CholmodKKTSolver(P::SparseMatrixCSC{Tv, Ti}, A::SparseMatrixCSC{Tv, Ti}, sigma::Tv, rho) where {Tv <: AbstractFloat, Ti <: Integer}

         m,n  = _kktutils_check_dims(P, A, sigma, rho)
         K    = _assemble_kkt_full(P, A, sigma, rho)
         fact = ldlt(K)

        return new{Tv, Ti}(fact, K, m, n)

    end
end

solve!(s::CholmodKKTSolver, lhs, rhs) = lhs .= s.fact \ rhs

function update_rho!(s::CholmodKKTSolver{Tv, Ti}, rho::Union{Tv, AbstractVector{Tv}}) where {Tv <: AbstractFloat, Ti <: Integer}
    update_kkt_matrix!(s.K, s.n, s.m, rho)
    #complete restart
    s.fact = ldlt(s.K)
end

free_memory!(s::AbstractKKTSolver) = nothing


text:
using ComparativeAutograder

@testset "Test runtestsuite" begin
    @testset "test suite for println" begin
        suite = TestSuite([
            TestCase("foo"),
        ])
        output = runtestsuite(println, suite)
        result = output.results[1]
        @test result.result === nothing
        @test result.exception === nothing
        @test result.stdout == "foo\n"
        @test result.stderr == ""
        @test output.elapsed_time > 0.0
    end

    @testset "Test runtestsuite on function that throws" begin
        myfunc = (x::Bool) -> x ? error("Throwing!") : "Not throwing."
        suite = TestSuite([
            TestCase(false),
            TestCase(true),
        ])
        output = runtestsuite(myfunc, suite)

        @test output.results[1].result == "Not throwing."
        @test output.results[2].exception !== nothing
    end
end


text:
using DiffEqSensitivity, SafeTestsets
using Test

const GROUP = get(ENV, "GROUP", "All")
const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR")
const is_TRAVIS = haskey(ENV,"TRAVIS")

@time begin
if GROUP == "All" || GROUP == "Core" || GROUP == "Downstream"
    @time @safetestset "Forward Sensitivity" begin include("forward.jl") end
    @time @safetestset "Adjoint Sensitivity" begin include("adjoint.jl") end
    @time @safetestset "Second Order Sensitivity" begin include("second_order.jl") end
    @time @safetestset "Concrete Solve Derivatives" begin include("concrete_solve_derivatives.jl") end
end

if GROUP == "All" || GROUP == "GSA"
    @time @safetestset "Morris Method" begin include("morris_method.jl") end
    @time @safetestset "Sobol Method" begin include("sobol_method.jl") end
    @time @safetestset "DGSM Method" begin include("DGSM.jl") end
    @time @safetestset "eFAST Method" begin include("eFAST_method.jl") end
    @time @safetestset "RegressionGSA Method" begin include("regression_sensitivity.jl") end
end

if GROUP == "DiffEqFlux"
    using Pkg
    if is_TRAVIS
      using Pkg
      Pkg.add("DiffEqFlux")
    end
    @time Pkg.test("DiffEqFlux")
end
end