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

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

# ╔═╡ e0c0dc94-277e-11eb-379e-83d064a93413
begin
	using Plots
	using PlutoUI
end

# ╔═╡ 00877a4a-277c-11eb-3ec0-e71e4094b404
md"""
# Evolution in time and space: Advection and diffusion in 1D
"""

# ╔═╡ 1b25b916-277c-11eb-0274-4b4fb946258d
md"""
So far we have been looking at dynamics in time, for example how does the temperature of the Earth change over time. But the Earth does not have a single, uniform temperature; rather, at a particular moment in time, different places on Earth are at different temperatures, and those different temperatures change over time due to several mechanisms.

In this notebook we will look at two fundamental mechanisms: **advection** and **diffusion**. Let's think about the temperature in the ocean. Since the ocean is a fluid that is in motion, a warm "parcel" of water can flow (move) to a new location due to the physical motion of the water itself; this is **advection**.  And even if the water doesn't move, temperature or a higher concentration of some substance dissolved in the fluid can spread out due to molecular mechanisms; this is **diffusion**.
"""

# ╔═╡ 956f5104-277d-11eb-291d-1faef485a5aa
md"""
In this notebook we will restrict ourselves to 1 space dimension.
So we will think about the temperature $T$, for example, being a function

$$T(t, x)$$

of

- time: $t$
- space: $x$

We want to calculate a value of the temperature $T$ at each point $(t, x)$.

"""

# ╔═╡ 36328b70-277d-11eb-02c7-2f854c1466cc
md"""
## Temperature profiles and discretization
"""

# ╔═╡ 42190984-277d-11eb-1ac2-7d84516c3269
md"""
An ordinary differential equation needs an initial value for each variable. Similarly, we will need an initial *function* $T_0(x)$ that gives us the temperature at each position $x$. Let's suppose that the position is restricted to the interval $[0, L_x]$.

As usual, to represent the continuous function $T_0(x)$ on a computer, we will need to **discretise** it in some way, i.e. *approximate* the continuous function by a *finite* set of numbers in the computer.

The simplest (but far from the only!) discretisation method is to **sample** the function at discrete **grid points** (or **nodes**) $x_i$, for $i = 1, \dots, N_x$. For simplicity we will take these equally spaced, with spacing $x_{i+1} - x_i =: \delta x := L_x / N_x$.
"""

# ╔═╡ d2bed768-277e-11eb-32cf-41f1fedec3cb
md"""
For example, let's consider an initial temperature profile given by a sine curve:
"""

# ╔═╡ e6493da0-277e-11eb-22ff-29752652b576
T₀(x) = sin(2π * x)

# ╔═╡ 0d6271c0-286c-11eb-1c9c-3ba039b49d24
md"""
and define the grid points as follows:
"""

# ╔═╡ f17f7734-277e-11eb-25cf-5f2ba2db5aa3
begin
	Nₓ = 20
	Lₓ = 1.0

	δx = Lₓ / Nₓ

	xs = δx/2:δx:Lₓ
end

# ╔═╡ fa327c08-286b-11eb-0032-2384998a42a8
xs

# ╔═╡ 0db43be2-284c-11eb-2740-4379437fd70c
md"""
It turns out to be a good idea to take the grid points at the *centre* of each interval, so that we have $N_x$ intervals and $N_x$ grid points, starting at $x_1 = \delta x/2$ and finishing at $x_N = L_x - \delta x / 2$.
"""

# ╔═╡ 468a0590-2780-11eb-045c-d1f468fc4e50
md"""
We call such a function of $x$ at a given time a **temperature profile**. Let's draw it as a function and as a heatmap:
"""

# ╔═╡ af30a0d0-2781-11eb-0274-ab423205facb
md"""
We will denote by $T^0_i$ the initial temperature at grid point number $i$.
"""

# ╔═╡ 646bc32e-284c-11eb-2ce8-5f64b1a49534
md"""
A useful way to think about $T^n_i$ is as the (spatial) average of $T(t_n, x)$ over the interval of positions between neighbouring grid points, so $T_i$ is the average over the interval between $x_i - \frac{\delta x}{2}$ and $x_i + \frac{\delta x}{2}$. We can thus think of the following **piecewise constant** approximation to the original continuous function:
"""

# ╔═╡ 79ce4b10-284c-11eb-2258-2155f850171d
let
	δx = xs[2] - xs[1]

	p = plot(0:0.001:Lₓ, T₀, label="T₀", lw=1, ls=:dash)
	scatter!(xs, T₀.(xs), label="sampled")
	scatter!(xs, zero.(xs), label="x nodes", alpha=0.5, ms=3, lw=2)

	for i in 1:length(xs)
		plot!([ (xs[i] - δx/2, T₀(xs[i])), (xs[i] + δx/2, T₀(xs[i])) ], c=:green, lw=4, lab=false)

		plot!([ (xs[i] - δx/2, 0), (xs[i] - δx/2, T₀(xs[i])), (xs[i] + δx/2, T₀(xs[i])), (xs[i] + δx/2, 0)], c=:green, lw=1, lab=false, ls=:dash, alpha=0.3
		)
	end

	xlabel!("x")
	ylabel!("T₀(x)")
end

# ╔═╡ 2494daaa-2780-11eb-3084-2317924048ea
md"""
## Advection
"""

# ╔═╡ 29444ffe-2780-11eb-0875-095302b5d486
md"""
Now let's think of this as the temperature in a fluid, and let's suppose that the fluid is moving to the right with a constant, uniform speed $U$. (**Uniform** here means that the speed is the same in all parts of the fluid.) Then the temperature profile should also *move with the fluid*. We call a quantity, such as the temperature, that is carried along with the fluid a **tracer**.

If we fix our attention at a single point in space, say the grid point $x_i$, the temperature there will vary over time due to the motion of the fluid. How it varies in time depends on the values at neighbouring grid points, since they determine how much heat will be *transported into* the current cell.
"""

# ╔═╡ 7256778a-2785-11eb-0369-f3b43d5dd203
md"""
Let's think about a time-stepping algorithm, and let $\delta t$ be the duration of a time step. We'll denote by $T^n_i$ the approximate value of $T$ to calculate at position $x_i$ and at the $n$th time step $t_n$, i.e. an approximation of $T(t_n, x_i)$, where $t_n = n \, \delta t$.

For example, $T^{n+1}_i \simeq T(t_n + \delta t, x_i)$ and $T^{n}_{i+1} \simeq T(t_n, x_i + \delta x)$.

[Note that the superscript $n$ in these algorithms does not denote a power; it's just a label for the time step. We could write $T_i^{(n)}$ instead, but that is annoying to both write and read, so we suppress the parentheses.]
"""

# ╔═╡ 44433a34-2782-11eb-0079-837c9306c5bd
md"""
Suppose the fluid is moving to the right with speed $U$. During a time step of duration $\delta t$, the temperature $T^n_i$ at cell $i$ changes due to the amount of heat that leaves cell $i$, and the amount that enters that cell.

Most of the fluid that starts within cell $i$ remains within that cell during the time step (if the time step is short enough). Only fluid that is within a distance $U \, \delta t$ from the boundary of the cell will cross over that cell. So a *proportion* $(U \, \delta t) / \delta x$ of the amount in cell $i$ crosses the boundary.

Hence, roughly an amount $T^n_i (U \delta t) / \delta x$ will leave cell number $i$ and cross into cell $i+1$ (the cell to the right). Similarly, an amount $T^n_{i-1} (U \delta t) / \delta x$ will *enter* cell $i$ from the neighbouring cell $i-1$ on the left. Hence we have

$$T^{n+1}_i = T^{n}_i + (T^n_{i-1} - T^n_{i})\, U \, \delta t / \delta x.$$
"""

# ╔═╡ 72c0ab0c-2781-11eb-1f59-9b22a52b0be0
md"""
Thus

$$\frac{T^{n+1}_i - T^{n}_i}{\delta t} =  \frac{T^n_{i-1} - T^n_{i}}{\delta x}\, U.$$
"""

# ╔═╡ e5761990-278b-11eb-134e-7954b577b1ac
md"""
Taking the continuum limit when $\delta t \to 0$ and $\delta x \to 0$, we recognise the definition of **partial derivatives** with respect to time and space variables from multivariable calculus. (Note the different indices that change on the two sides of the equation.)

Denoting these partial derivatives using $\partial$, we arrive at the **advection equation**:

$$\frac{\partial T(t, x)}{\partial t} = -U \frac{\partial T(t, x)}{\partial x},$$

or for short

$$\frac{\partial T}{\partial t} = -U \frac{\partial T}{\partial x}.$$

Since $T$ is a function of both $x$ and $t$, and this equation involves partial derivatives with respect to both of the independent variables, this is a **partial differential equation** (PDE). It describes how the function $T(t, x)$ changes continuously as a function both of time and space.

Although there are some analytical methods to solve PDEs, often it's necessary to use numerical methods. Here we'll look at simple numerical methods to solve such equations.
"""

# ╔═╡ 2033364e-278c-11eb-2936-17598ce14a41
md"""
## Numerics for the advection equation
"""

# ╔═╡ e9a37908-278c-11eb-278e-9bd155f0cae6
md"""
Let's return to the version of the equation in which the value at the *following* time step is isolated:

$$T^{n+1}_i = T^{n}_i - \left( U \frac{\delta t}{\delta x} \right) (T^n_{i} - T^n_{i-1}).$$

In the last term on the right-hand side, we see that we require combinations of values of $T$ at the *same* time step from *different* places, with certain coefficients.
"""

# ╔═╡ a10266c4-284e-11eb-3433-adee2e0e2d52
md"""
There are many approaches to implementing this numerically. The simplest is to directly transcribe the equation for the $i$th entry of the vector.

Calling `T` the current vector, i.e. $\mathbf{T}^n := (T^n_i)_{i=1, \ldots, N_x}$, and `T′` the new vector at the next time step, we have the following basic expression:

	T′[i] = T[i] + δt * U * (T[i-1] - T[i]) / δx

But now we realise a problem: What should we do when $i=1$? This will try to access the index 0 of the vector `T`, which does not exist!

This illustrates the necessity of choosing **boundary conditions** that specify what happens at the edge of the domain.

For simplicity we will choose to use **periodic boundary conditions**. This is a convenient mathematical fiction that allows us to treat all cells as being on the same footing, by wrapping the system around a torus, so that cells $i=1$ and $i=N_x$ are neighbours.
"""

# ╔═╡ e542a8da-284e-11eb-3297-6bbbf052284b
md"""
We can then write this as follows, where we separate out the case $i=1$:
"""

# ╔═╡ b15f4f44-284b-11eb-37c5-ab0153f7fe92
function advection(T, δt, δx, U)
	N = length(T)
	T′ = similar(T)  # create new vector of the same length

	for i in 2:N
		T′[i] = T[i] - δt * U * (T[i] - T[i-1]) / δx
	end

	# periodic boundary:
	T′[1] = T[1] - δt * U * (T[1] - T[N]) / δx

	return T′
end

# ╔═╡ 85540de0-286e-11eb-0d00-03435a8c7c83
md"""
We use a function `evolve` that carries out a generic time evolution. This is basically the Euler method again, except now we apply it to a vector of values at different points in space at each step.
"""

# ╔═╡ af79e360-286e-11eb-2a4d-3d6d7564088c
δt = 0.001;

# ╔═╡ dce9e53a-28f4-11eb-070b-17e10779a38b
U = 0.2;

# ╔═╡ 8c05e3cc-2858-11eb-1e1c-9781c30738c3
md"""
Unfortunately this does *not* behave as we expect: instead of preserving the shape of the profile over time, it is decaying. This is due to the way we are approximating.

A better way to discretize the spatial derivative is using the following **centered difference**:

$$\frac{\partial T(t_n, x_i)}{\partial x} \simeq \frac{T^n_{i+1} - T^n_{i-1}}{2 \delta x}$$

"""

# ╔═╡ a29fecac-285a-11eb-14b0-9313f8994fbb
function advection2(T, δt, δx, U)
	N = length(T)
	T′ = similar(T)  # create new vector of the same length

	for i in 2:N-1
		T′[i] = T[i] - δt * U * (T[i+1] - T[i-1]) / (2δx)
	end

	# periodic boundary:
	T′[1] = T[1] - δt * U * (T[2] - T[N]) / (2δx)
	T′[N] = T[N] - δt * U * (T[1] - T[N-1]) / (2δx)

	return T′
end

# ╔═╡ c59388ea-286e-11eb-0f21-eb18e5ba516f
md"""
## Diffusion
"""

# ╔═╡ 3c944998-2888-11eb-087d-492b9d0ee32e
md"""
Another key physical process is **diffusion**. This models how temperature or mass spreads out from hot or high concentration regions towards regions where it is cold or where there is a low concentration.

The PDE describing this is the **heat equation**, or **diffusion equation**:
"""

# ╔═╡ 6b7cea44-2888-11eb-0208-990860d6a152
md"""
$$\frac{\partial T}{\partial t} = D \frac{\partial^2 T}{\partial x^2}.$$
"""

# ╔═╡ 68db3372-2888-11eb-1b03-b5ebca4c2bd5
md"""
To obtain a numerical method to solve this equation, we again need to discretise this, in particular the second derivative. One possible discretisation is

$$\frac{\partial^2 T}{\partial x^2}(t_n, x_i) \simeq \frac{T^n_{i+1} - 2 T^n_i + T^n_{i-1}}{\delta x^2}.$$
"""

# ╔═╡ d6131ad0-2889-11eb-3085-15d17e33ee7a
md"""
This may again be transcribed directly into code:
"""

# ╔═╡ 630314bc-2868-11eb-1b93-b7b08a4b2887
function diffusion(T, δt, δx, D)
	N = length(T)
	T′ = similar(T)  # create new vector of the same length

	for i in 2:N-1
		T′[i] = T[i] + δt * D * (T[i+1] -2T[i] + T[i-1]) / (δx^2)
	end

	# periodic boundary:
	T′[1] = T[1] + δt * D * (T[2] - 2T[1] + T[N]) / (δx^2)
	T′[N] = T[N] + δt * D * (T[1] - 2T[N] + T[N-1]) / (δx^2)

	return T′
end

# ╔═╡ e63cfa84-2889-11eb-1ea2-51726645ddd9
md"""
## Advection--diffusion
"""

# ╔═╡ eee3008e-2889-11eb-088a-73aff304e736
md"""
Finally we can combine both mechanisms, to describe a tracer that is both being advected at a constant speed and diffusing. This basically utilises the composition of the advection and diffusion functions:
"""

# ╔═╡ ffd2a838-2889-11eb-1a7c-b35992543b8a
function advection_diffusion(T, δt, δx, (U, D))
	temp = advection2(T, δt, δx, U)
	return diffusion(temp, δt, δx, D)
end

# ╔═╡ 575a5f3c-2780-11eb-2119-27a4114ceac5
md"""
## Function library
"""

# ╔═╡ 5a3eec86-2780-11eb-0341-39a5c343fc52
function temperature_heatmap(x, T)

	p = heatmap(x, [0.], collect(T'),
			   clims=(-1., 1.), cbar=false, xticks=nothing, yticks=nothing)

	return p
end

# ╔═╡ 6de1859c-277f-11eb-1ead-8b4794832d59
begin
	p1 = plot(0:0.001:Lₓ, T₀, label="T₀", lw=3)
	scatter!(xs, T₀.(xs), label="sampled")
	scatter!(xs, zero.(xs), label="x nodes", alpha=0.5, ms=3)

	xlabel!("x")
	ylabel!("T₀")

	for x in xs
		plot!([ (x, 0), (x, T₀(x)) ], ls=:dash, c=:black, label="", alpha=0.5)
	end

	hline!([0], ls=:dash, lab=false)

	p2 = temperature_heatmap(xs, T₀.(xs))

	plot(p1, p2, layout = grid(2, 1, heights=[0.9, 0.1]))

end

# ╔═╡ 9187350a-2851-11eb-05f0-d3a6eef190fe
function evolve(method, xs, δt, U, t_final=10.0, f₀=T₀)

	T = f₀.(xs)
	δx = xs[2] - xs[1]

	t = 0.0
	ts = [t]

	results = [T]

	while t < t_final
		T′ = method(T, δt, δx, U)  # new
		push!(results, T′)

		t += δt
		push!(ts, t)

		T = copy(T′)

	end

	return ts, results
end

# ╔═╡ 02a893e4-2852-11eb-358a-371459191da7
ts, evolution = evolve(advection, xs, δt, U)

# ╔═╡ e6ae447e-2851-11eb-3fe1-096459167f2b
@bind n Slider(1:length(evolution), show_value=true)

# ╔═╡ 014e2530-2852-11eb-103f-1d647cb999b0
let
	p1 = plot(xs, evolution[n], m=:o, xlim=(0, 1), ylim=(-1.1, 1.1), title="t = $(round(ts[n], digits=2))", leg=false)

	p2 = temperature_heatmap(xs, evolution[n])

	plot(p1, p2, layout = grid(2, 1, heights=[0.9, 0.1]))
end

# ╔═╡ e42ec13e-285a-11eb-3cc0-7dc41ed5495b
ts2, evolution2 = evolve(advection2, xs, δt, 0.1)

# ╔═╡ f60a8b5e-285a-11eb-0d35-8daf23cf92ae
@bind n2 Slider(1:length(evolution2), show_value=true)

# ╔═╡ f1b5d130-285a-11eb-001c-67035925f43d
let
	p1 = plot(xs, evolution2[n2], m=:o, xlim=(0, 1), ylim=(-1.1, 1.1), title="t = $(round(ts2[n2], digits=2))", leg=false)

	p2 = temperature_heatmap(xs, evolution2[n2])

	plot(p1, p2, layout = grid(2, 1, heights=[0.9, 0.1]))
end

# ╔═╡ 09bc3c40-288a-11eb-0339-59f0b70e03a3
@bind n3 Slider(1:length(evolution2), show_value=true)

# ╔═╡ 121255d2-288a-11eb-1fa5-9db68af8c232
ts3, evolution3 = evolve(diffusion, xs, δt, 0.01)

# ╔═╡ 175d9902-288a-11eb-3700-390ccd1caa5b
let
	p1 = plot(xs, evolution3[n3], m=:o, xlim=(0, 1), ylim=(-1.1, 1.1), title="t = $(round(ts3[n3], digits=2))", leg=false)
	p2 = temperature_heatmap(xs, evolution3[n3])

	plot(p1, p2, layout = grid(2, 1, heights=[0.9, 0.1]), clim=(-1, 1))
end

# ╔═╡ 73efeed4-28fa-11eb-10c0-8fcfd53b812f
ts4, evolution4 = evolve(diffusion, xs, 0.0001, 0.01, 10.0, x -> (0.45 < x < 0.55) ? 1.0 : 0.0)

# ╔═╡ e74d920a-28fa-11eb-3c91-9133a01effc5
@bind n4 Slider(1:length(evolution4), show_value=true)

# ╔═╡ dc7b6328-28fa-11eb-38b7-71a6f8d0a751
let
	p1 = plot(xs, evolution4[n4], m=:o, xlim=(0, 1), ylim=(-0.1, 1.1), title="t = $(round(ts3[n3], digits=2))", leg=false)
	p2 = temperature_heatmap(xs, evolution4[n4])

	plot(p1, p2, layout = grid(2, 1, heights=[0.9, 0.1]), clim=(-1, 1))
end

# ╔═╡ f6fa3770-288d-11eb-32de-f95e03705791
ts5, evolution5 = evolve(advection_diffusion, xs, δt, (1.0, 0.01))

# ╔═╡ 6eb00a02-288d-11eb-354b-b56cf5a8380e
@bind n5 Slider(1:length(evolution5), show_value=true)

# ╔═╡ 65126bfc-288d-11eb-2bfc-493588365164
let
	p1 = plot(xs, evolution5[n5], m=:o, xlim=(0, 1), ylim=(-1.1, 1.1), title="t = $(round(ts3[n3], digits=2))", leg=false)
	p2 = temperature_heatmap(xs, evolution5[n5])

	plot(p1, p2, layout = grid(2, 1, heights=[0.9, 0.1]), clim=(-1, 1))
end

# ╔═╡ Cell order:
# ╠═e0c0dc94-277e-11eb-379e-83d064a93413
# ╟─00877a4a-277c-11eb-3ec0-e71e4094b404
# ╟─1b25b916-277c-11eb-0274-4b4fb946258d
# ╟─956f5104-277d-11eb-291d-1faef485a5aa
# ╟─36328b70-277d-11eb-02c7-2f854c1466cc
# ╟─42190984-277d-11eb-1ac2-7d84516c3269
# ╟─d2bed768-277e-11eb-32cf-41f1fedec3cb
# ╠═e6493da0-277e-11eb-22ff-29752652b576
# ╟─0d6271c0-286c-11eb-1c9c-3ba039b49d24
# ╠═f17f7734-277e-11eb-25cf-5f2ba2db5aa3
# ╠═fa327c08-286b-11eb-0032-2384998a42a8
# ╟─0db43be2-284c-11eb-2740-4379437fd70c
# ╟─468a0590-2780-11eb-045c-d1f468fc4e50
# ╟─6de1859c-277f-11eb-1ead-8b4794832d59
# ╟─af30a0d0-2781-11eb-0274-ab423205facb
# ╟─646bc32e-284c-11eb-2ce8-5f64b1a49534
# ╟─79ce4b10-284c-11eb-2258-2155f850171d
# ╟─2494daaa-2780-11eb-3084-2317924048ea
# ╟─29444ffe-2780-11eb-0875-095302b5d486
# ╟─7256778a-2785-11eb-0369-f3b43d5dd203
# ╟─44433a34-2782-11eb-0079-837c9306c5bd
# ╟─72c0ab0c-2781-11eb-1f59-9b22a52b0be0
# ╟─e5761990-278b-11eb-134e-7954b577b1ac
# ╟─2033364e-278c-11eb-2936-17598ce14a41
# ╟─e9a37908-278c-11eb-278e-9bd155f0cae6
# ╟─a10266c4-284e-11eb-3433-adee2e0e2d52
# ╟─e542a8da-284e-11eb-3297-6bbbf052284b
# ╠═b15f4f44-284b-11eb-37c5-ab0153f7fe92
# ╟─85540de0-286e-11eb-0d00-03435a8c7c83
# ╠═af79e360-286e-11eb-2a4d-3d6d7564088c
# ╠═dce9e53a-28f4-11eb-070b-17e10779a38b
# ╠═02a893e4-2852-11eb-358a-371459191da7
# ╠═e6ae447e-2851-11eb-3fe1-096459167f2b
# ╟─014e2530-2852-11eb-103f-1d647cb999b0
# ╟─8c05e3cc-2858-11eb-1e1c-9781c30738c3
# ╠═a29fecac-285a-11eb-14b0-9313f8994fbb
# ╠═e42ec13e-285a-11eb-3cc0-7dc41ed5495b
# ╠═f60a8b5e-285a-11eb-0d35-8daf23cf92ae
# ╟─f1b5d130-285a-11eb-001c-67035925f43d
# ╟─c59388ea-286e-11eb-0f21-eb18e5ba516f
# ╟─3c944998-2888-11eb-087d-492b9d0ee32e
# ╟─6b7cea44-2888-11eb-0208-990860d6a152
# ╟─68db3372-2888-11eb-1b03-b5ebca4c2bd5
# ╟─d6131ad0-2889-11eb-3085-15d17e33ee7a
# ╠═630314bc-2868-11eb-1b93-b7b08a4b2887
# ╠═121255d2-288a-11eb-1fa5-9db68af8c232
# ╠═09bc3c40-288a-11eb-0339-59f0b70e03a3
# ╟─175d9902-288a-11eb-3700-390ccd1caa5b
# ╠═73efeed4-28fa-11eb-10c0-8fcfd53b812f
# ╠═e74d920a-28fa-11eb-3c91-9133a01effc5
# ╟─dc7b6328-28fa-11eb-38b7-71a6f8d0a751
# ╟─e63cfa84-2889-11eb-1ea2-51726645ddd9
# ╟─eee3008e-2889-11eb-088a-73aff304e736
# ╠═ffd2a838-2889-11eb-1a7c-b35992543b8a
# ╠═f6fa3770-288d-11eb-32de-f95e03705791
# ╠═6eb00a02-288d-11eb-354b-b56cf5a8380e
# ╠═65126bfc-288d-11eb-2bfc-493588365164
# ╟─575a5f3c-2780-11eb-2119-27a4114ceac5
# ╟─5a3eec86-2780-11eb-0341-39a5c343fc52
# ╠═9187350a-2851-11eb-05f0-d3a6eef190fe

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using LinearAlgebra
using ValueShapes, Distributions, ArraysOfArrays
using ForwardDiff, Zygote, DistributionsAD
using InverseFunctions, ChangesOfVariables

@testset "test_distribution_transform" begin
    function test_back_and_forth(trg_d, src_d)
        @testset "transform $(typeof(trg_d).name) <-> $(typeof(src_d).name)" begin
            src_v = rand(src_d)
            trg_v = BAT.apply_dist_trafo(trg_d, src_d, src_v)
            src_v_reco = BAT.apply_dist_trafo(src_d, trg_d, trg_v)

            @test src_v ≈ src_v_reco

            let vs_trg = varshape(trg_d), vs_src = varshape(src_d)
                f = unshaped_x -> inverse(vs_trg)(BAT.apply_dist_trafo(trg_d, src_d, vs_src(unshaped_x)))
                ref_ladj = logpdf(src_d, src_v) - logpdf(trg_d, trg_v)
                @test ref_ladj ≈ logabsdet(ForwardDiff.jacobian(f, inverse(vs_src)(src_v)))[1]
            end
        end
    end

    function get_trgxs(trg_d, src_d, X)
        return (x -> BAT.apply_dist_trafo(trg_d, src_d, x)).(nestedview(X))
    end

    function get_trgxs(trg_d, src_d::Distribution{Univariate}, X)
        return (x -> BAT.apply_dist_trafo(trg_d, src_d, x)).(X)
    end

    function test_dist_trafo_moments(trg_d, src_d)
        @testset "check moments of trafo $(typeof(trg_d).name) <- $(typeof(src_d).name)" begin
            X = flatview(rand(src_d, 10^5))
            trgxs = get_trgxs(trg_d, src_d, X)
            unshaped_trgxs = broadcast(unshaped, trgxs, Ref(varshape(trg_d)))
            @test isapprox(mean(unshaped_trgxs), mean(unshaped(trg_d)), atol = 0.1)
            @test isapprox(cov(unshaped_trgxs), cov(unshaped(trg_d)), rtol = 0.1)
        end
    end

    stduvuni = BAT.StandardUvUniform()
    stduvnorm = BAT.StandardUvNormal()

    uniform1 = Uniform(-5.0, -0.01)
    uniform2 = Uniform(0.01, 5.0)

    normal1 = Normal(-10, 1)
    normal2 = Normal(10, 5)

    stdmvnorm1 = BAT.StandardMvNormal(1)
    stdmvnorm2 = BAT.StandardMvNormal(2)

    stdmvuni2 = BAT.StandardMvUniform(2)

    standnorm2_reshaped = ReshapedDist(stdmvnorm2, varshape(stdmvnorm2))

    mvnorm = MvNormal([0.3, -2.9], [1.7 0.5; 0.5 2.3])
    beta = Beta(3,1)
    gamma = Gamma(0.1,0.7)
    dirich = Dirichlet([0.1,4])

    ntdist = NamedTupleDist(
        a = uniform1,
        b = mvnorm,
        c = [4.2, 3.7],
        x = beta,
        y = gamma
    )

    test_back_and_forth(stduvuni, stduvuni)
    test_back_and_forth(stduvnorm, stduvnorm)
    test_back_and_forth(stduvuni, stduvnorm)
    test_back_and_forth(stduvnorm, stduvuni)

    test_back_and_forth(stdmvuni2, stdmvuni2)
    test_back_and_forth(stdmvnorm2, stdmvnorm2)
    test_back_and_forth(stdmvuni2, stdmvnorm2)
    test_back_and_forth(stdmvnorm2, stdmvuni2)

    test_back_and_forth(beta, stduvnorm)
    test_back_and_forth(gamma, stduvnorm)

    test_dist_trafo_moments(normal2, normal1)
    test_dist_trafo_moments(uniform2, uniform1)

    test_dist_trafo_moments(beta, gamma)

    test_dist_trafo_moments(beta, stduvnorm)
    test_dist_trafo_moments(gamma, stduvnorm)

    test_dist_trafo_moments(mvnorm, stdmvnorm2)
    test_dist_trafo_moments(dirich, stdmvnorm1)

    test_dist_trafo_moments(mvnorm, stdmvuni2)
    test_dist_trafo_moments(stdmvuni2, mvnorm)

    test_dist_trafo_moments(stdmvnorm2, stdmvuni2)

    test_dist_trafo_moments(mvnorm, standnorm2_reshaped)
    test_dist_trafo_moments(standnorm2_reshaped, mvnorm)
    test_dist_trafo_moments(stdmvnorm2, standnorm2_reshaped)
    test_dist_trafo_moments(standnorm2_reshaped, standnorm2_reshaped)

    test_back_and_forth(ntdist, BAT.StandardMvNormal(5))
    test_back_and_forth(ntdist, BAT.StandardMvUniform(5))

    let
        mvuni = product_distribution([Uniform(), Uniform()])

        x = rand()
        @test_throws ArgumentError BAT.apply_dist_trafo(stduvnorm, mvnorm, x)
        @test_throws ArgumentError BAT.apply_dist_trafo(stduvnorm, stdmvnorm1, x)
        @test_throws ArgumentError BAT.apply_dist_trafo(stduvnorm, stdmvnorm2, x)

        x = rand(2)
        @test_throws ArgumentError BAT.apply_dist_trafo(mvuni, mvnorm, x)
        @test_throws ArgumentError BAT.apply_dist_trafo(mvnorm, mvuni, x)
        @test_throws ArgumentError BAT.apply_dist_trafo(stduvnorm, mvnorm, x)
        @test_throws ArgumentError BAT.apply_dist_trafo(stduvnorm, stdmvnorm1, x)
        @test_throws ArgumentError BAT.apply_dist_trafo(stduvnorm, stdmvnorm2, x)
    end

    let
        primary_dist = NamedTupleDist(x = Normal(2), c = 5)
        f = x -> NamedTupleDist(y = Normal(x.x, 3), z = MvNormal([1.3 0.5; 0.5 2.2]))
        trg_d = @inferred(HierarchicalDistribution(f, primary_dist))
        src_d = BAT.StandardMvNormal(totalndof(varshape(trg_d)))
        test_back_and_forth(trg_d, src_d)
        test_dist_trafo_moments(trg_d, src_d)
    end

    #=
    using Cuba
    function integrate_over_unit(density::AbstractDensity)
        vs = varshape(density)
        f_cuba(source_x, y) = y[1] = exp(logdensityof(density)(vs(source_x)))
        Cuba.vegas(f_cuba, 1, 1).integral[1]
    end
    =#

    @testset "Custom cdf and quantile for dual numbers" begin
        Dual = ForwardDiff.Dual

        @test BAT._trafo_cdf(Normal(Dual(0, 1, 0, 0), Dual(1, 0, 1, 0)), Dual(0.5, 0, 0, 1)) == cdf(Normal(Dual(0, 1, 0, 0), Dual(1, 0, 1, 0)), Dual(0.5, 0, 0, 1))
        @test BAT._trafo_cdf(Normal(0, 1), Dual(0.5, 1)) == cdf(Normal(0, 1), Dual(0.5, 1))

        @test BAT._trafo_quantile(Normal(0, 1), Dual(0.5, 1)) == quantile(Normal(0, 1), Dual(0.5, 1))
        @test BAT._trafo_quantile(Normal(Dual(0, 1, 0, 0), Dual(1, 0, 1, 0)), Dual(0.5, 0, 0, 1)) == quantile(Normal(Dual(0, 1, 0, 0), Dual(1, 0, 1, 0)), Dual(0.5, 0, 0, 1))
    end

    for VT in (NamedTuple, ShapedAsNT)
        src_dist = unshaped(NamedTupleDist(VT, a = Weibull(), b = MvNormal([1.3 0.6; 0.6 2.4])))
        f = BAT.DistributionTransform(Normal, src_dist)
        x = rand(src_dist)
        InverseFunctions.test_inverse(f, x)
        ChangesOfVariables.test_with_logabsdet_jacobian(f, x, ForwardDiff.jacobian)
    end

    @testset "trafo broadcasting" begin
        dist = NamedTupleDist(a = Weibull(), b = Exponential())
        smpls = bat_sample(dist, IIDSampling(nsamples = 100)).result
        trafo = BAT.DistributionTransform(Normal, dist)
        @inferred(broadcast(trafo, smpls)) isa DensitySampleVector
        smpls_tr = trafo.(smpls)
        smpls_tr_cmp = [trafo(s) for s in smpls]
        @test smpls_tr == smpls_tr_cmp
	    @test @inferred(varshape(trafo)) == varshape(trafo.source_dist)
	    @test @inferred(trafo(varshape(trafo))) == varshape(trafo.target_dist)
    end

    @testset "trafo composition" begin
        dist1 = @inferred(NamedTupleDist(a = Normal(), b = Uniform(), c = Cauchy()))
        dist2 = @inferred(NamedTupleDist(a = Exponential(), b = Weibull(), c = Beta()))
        normal1 = Normal()
        normal2 = Normal(2)

        trafo = @inferred(BAT.DistributionTransform(dist1, dist2))
        inv_trafo = @inferred(inverse(trafo))

        composed_trafo = @inferred(∘(trafo, inv_trafo))
        @test composed_trafo.source_dist == composed_trafo.target_dist == dist1
        @test composed_trafo ∘ trafo == trafo
        @test_throws ArgumentError  trafo ∘ composed_trafo

        trafo = @inferred(BAT.DistributionTransform(normal1, normal2))
        @test_throws ArgumentError trafo ∘ trafo
    end

    @testset "full density transform" begin
        likelihood = NamedTupleDist(a = Normal(), b = Exponential())
        prior = NamedTupleDist(a = Normal(), b = Gamma())
        posterior_density = PosteriorDensity(likelihood, prior)

        posterior_density_trafod = @inferred(bat_transform(PriorToUniform(), posterior_density, FullDensityTransform()))

        @test posterior_density_trafod.result.orig.likelihood.dist == likelihood
        @test posterior_density_trafod.result.orig.prior.dist == prior

        @test posterior_density_trafod.result.trafo.target_dist isa BAT.StandardMvUniform

        lower_bounds = Float32.([-10, -10, -10])
        upper_bounds = Float32.([10, 10, 10])
        rect_bounds = @inferred(BAT.HyperRectBounds(lower_bounds, upper_bounds))
            mvn = @inferred(product_distribution(Normal.(randn(3))))
        dist_density = @inferred(BAT.DistributionDensity(mvn, rect_bounds))

        dist_density_trafod = @inferred(bat_transform(PriorToUniform(), dist_density, FullDensityTransform()))

        @test dist_density_trafod.trafo.target_dist isa BAT.StandardMvUniform
        @test dist_density_trafod.result.orig == dist_density
        @test dist_density_trafod.trafo.source_dist == dist_density_trafod.result.trafo.source_dist == mvn

        @test dist_density_trafod.trafo(varshape(dist_density_trafod.trafo)) == @inferred(varshape(dist_density))

        dist_density_trafod = @inferred(bat_transform(PriorToGaussian(), dist_density, FullDensityTransform()))

        @test dist_density_trafod.trafo.target_dist isa BAT.StandardMvNormal
        @test dist_density_trafod.result.orig == dist_density
        @test dist_density_trafod.trafo.source_dist == dist_density_trafod.result.trafo.source_dist == mvn

        @test dist_density_trafod.trafo(varshape(dist_density_trafod.trafo)) == @inferred(varshape(dist_density))
    end

    @testset "trafo autodiff pullbacks" begin
        # ToDo: Test for type stability and fix where necessary.

        xs = rand(5)
        @test Zygote.jacobian(BAT._pushfront, xs, 42)[1] ≈ ForwardDiff.jacobian(xs -> BAT._pushfront(xs, 1), xs)
        @test Zygote.jacobian(BAT._pushfront, xs, 42)[2] ≈ vec(ForwardDiff.jacobian(x -> BAT._pushfront(xs, x[1]), [42]))
        @test Zygote.jacobian(BAT._pushback, xs, 42)[1] ≈ ForwardDiff.jacobian(xs -> BAT._pushback(xs, 1), xs)
        @test Zygote.jacobian(BAT._pushback, xs, 42)[2] ≈ vec(ForwardDiff.jacobian(x -> BAT._pushback(xs, x[1]), [42]))
        @test Zygote.jacobian(BAT._rev_cumsum, xs)[1] ≈ ForwardDiff.jacobian(BAT._rev_cumsum, xs)
        @test Zygote.jacobian(BAT._exp_cumsum_log, xs)[1] ≈ ForwardDiff.jacobian(BAT._exp_cumsum_log, xs) ≈ ForwardDiff.jacobian(cumprod, xs)

        src_v = [0.6, 0.7, 0.8, 0.9]
        f = inverse(BAT.DistributionTransform(Uniform, DistributionsAD.TuringDirichlet([3.0, 4.0, 5.0, 6.0, 7.0])))
        @test isapprox(ForwardDiff.jacobian(f, src_v), Zygote.jacobian(f, src_v)[1], rtol = 10^-4)
        f = inverse(BAT.DistributionTransform(Uniform, Dirichlet([3.0, 4.0, 5.0, 6.0, 7.0])))
        @test isapprox(ForwardDiff.jacobian(f, src_v), Zygote.jacobian(f, src_v)[1], rtol = 10^-4)
        f = inverse(BAT.DistributionTransform(Normal, Dirichlet([3.0, 4.0, 5.0, 6.0, 7.0])))
        @test isapprox(ForwardDiff.jacobian(f, src_v), Zygote.jacobian(f, src_v)[1], rtol = 10^-4)
    end
end

@testset "bat_transform_defaults" begin
    mvn = @inferred(product_distribution([Normal(-1), Normal(), Normal(1)]))
    uniform_prior = @inferred(product_distribution([Uniform(-3, 1), Uniform(-2, 2), Uniform(-1, 3)]))

    posterior_uniform_prior = @inferred(PosteriorDensity(mvn, uniform_prior))
    posterior_gaussian_prior = @inferred(PosteriorDensity(mvn, mvn))

    @test @inferred(bat_transform(PriorToGaussian(), posterior_uniform_prior)).result.prior.dist == @inferred(BAT.StandardMvNormal(3))
    @test @inferred(bat_transform(PriorToUniform(), posterior_gaussian_prior)).result.prior.dist == @inferred(BAT.StandardMvUniform(3))
    @test @inferred(bat_transform(NoDensityTransform(), posterior_uniform_prior)).result.prior.dist == uniform_prior
    pd = @inferred(product_distribution([Uniform() for i in 1:3]))
    density = @inferred(BAT.DistributionDensity(pd))
    @test @inferred(bat_transform(NoDensityTransform(), density)).result.dist == density.dist

    # ToDo: Improve comparison for bounds so `.dist` is not required here:
    @inferred(bat_transform(PriorToUniform(), convert(AbstractDensity, BAT.StandardUvUniform()))).result.dist == convert(AbstractDensity, BAT.StandardUvUniform()).dist
    @inferred(bat_transform(PriorToUniform(), convert(AbstractDensity, BAT.StandardMvUniform(4)))).result.dist == convert(AbstractDensity, BAT.StandardMvUniform(4)).dist
    @inferred(bat_transform(PriorToGaussian(), convert(AbstractDensity, BAT.StandardUvNormal()))).result.dist == convert(AbstractDensity, BAT.StandardUvNormal()).dist
    @inferred(bat_transform(PriorToGaussian(), convert(AbstractDensity, BAT.StandardMvNormal(4)))).result.dist == convert(AbstractDensity, BAT.StandardMvNormal(4)).dist
end

ext: jl


text:
# simple program to test MCL clustering

using Base.Test
using Clustering

srand(34568)

# test for nonsquare matrices
@test_throws DimensionMismatch mcl(zeros(Float64, 4, 3))

nodes = [:bat, :bit, :cat, :fit, :hat, :hit]
edges = Tuple{Symbol, Symbol, Float64}[(:cat, :hat, 0.2), (:hat, :bat, 0.16),
                                       (:bat, :cat, 1.0), (:bat, :bit, 0.125),
                                       (:bit, :fit, 0.25), (:fit, :hit, 0.5),
                                       (:hit, :bit, 0.16)]
adj_matrix = zeros(Float64, length(nodes), length(nodes))
for edge in edges
    n1 = findfirst(nodes, edge[1])
    n2 = findfirst(nodes, edge[2])
    adj_matrix[n1, n2] = adj_matrix[n2, n1] = edge[3]
end
@assert issymmetric(adj_matrix)

res = mcl(adj_matrix, display=:none, inflation=1.8)
@test isa(res, MCLResult)
k = length(res.counts)
# println("k = $k")
@test k == 2
@test all(res.assignments .<= k)
@test length(res.assignments) == length(nodes)
@test length(res.counts) == k
for c = 1:k
    @test countnz(res.assignments .== c) == res.counts[c]
end
@test res.nunassigned == 0
@test res.assignments == [1, 2, 1, 2, 1, 2]

# test with integer inflation value
res = mcl(adj_matrix, display=:none, inflation=2)
@test isa(res, MCLResult)
@test length(res.assignments) == length(nodes)
@test res.nunassigned == 0

# test non-integral expansion (show not raise an exception)
res = mcl(adj_matrix, display=:none, inflation=1.5, expansion=1.5, save_final_matrix=true)
@test isa(res, MCLResult)
@test length(res.assignments) == length(nodes)
@test size(res.mcl_adj) == size(adj_matrix) # test that the matrix is returned

# test allow_singles
res = mcl(diagm([1.0, 1.0]), display=:none, allow_singles=true)
@test length(res.counts) == 2
@test res.assignments == [1, 2]
@test res.counts == [1, 1]
@test res.nunassigned == 0

res = mcl(diagm([1.0, 1.0]), display=:none, allow_singles=false)
@test length(res.counts) == 0
@test res.assignments == [0, 0]
@test res.nunassigned == 2

# use sparse matrix
res = mcl(sparse(adj_matrix), display=:none, expansion=2)
@test isa(res, MCLResult)
@test length(res.assignments) == length(nodes)
@test res.nunassigned == 0
@test eltype(res.mcl_adj) === Float64

@test_throws ArgumentError mcl(sparse(adj_matrix), display=:none, expansion=2.1)

# use Float32 input
res = mcl(convert(Matrix{Float32},adj_matrix), display=:none, expansion=2)
@test isa(res, MCLResult)
@test eltype(res.mcl_adj) === Float32


text:
module LotkaVolterra

using Iterators
using Match
using ODE

import Base.show

export AbstractLotkaVolterra,
       communitymatrix,intrinsicrate,
       issector,sectors,
       LogisticGrowth,GeneralizedLotkaVolterra,
       rate,specificrate,jacobian,
       somefixedpoint,isfixedpoint,fixedpoint,fixedpoints,
       somefeasible,isfeasible,feasibility,
       somestable,isstable,stability,
       someviable,isviable,viability,
       odeint

include("utils.jl")
include("abstract.jl")
include("logistic.jl")
include("general.jl")

end


text:
let
    # before testing the model, check whether missing/null data is rejected
    tdata = [rand(1:10, 5)..., missing]
    @test_throws ErrorException("Data contains `missing` values and a fit cannot be performed") LsqFit.check_data_health(tdata, tdata)
    tdata = [rand(1:10, 5)..., Inf]
    @test_throws ErrorException("Data contains `Inf` or `NaN` values and a fit cannot be performed") LsqFit.check_data_health(tdata, tdata)
    tdata = [rand(1:10, 5)..., NaN]
    @test_throws ErrorException("Data contains `Inf` or `NaN` values and a fit cannot be performed") LsqFit.check_data_health(tdata, tdata)

    # fitting noisy data to an exponential model
    # TODO: Change to `.-x` when 0.5 support is dropped
    model(x, p) = p[1] .* exp.(-x .* p[2])

    # some example data
    Random.seed!(12345)
    xdata = range(0, stop=10, length=20)
    ydata = model(xdata, [1.0, 2.0]) + 0.01*randn(length(xdata))
    p0 = [0.5, 0.5]

    for ad in (:finite, :forward, :forwarddiff)
        fit = curve_fit(model, xdata, ydata, p0; autodiff = ad)
        @assert norm(fit.param - [1.0, 2.0]) < 0.05
        @test fit.converged

        # can also get error estimates on the fit parameters
        errors = margin_error(fit, 0.1)
        @assert norm(errors - [0.017, 0.075]) < 0.01
    end
    # if your model is differentiable, it can be faster and/or more accurate
    # to supply your own jacobian instead of using the finite difference
    function jacobian_model(x,p)
        J = Array{Float64}(undef, length(x), length(p))
        J[:,1] = exp.(-x.*p[2])     #dmodel/dp[1]
        J[:,2] = -x.*p[1].*J[:,1]           #dmodel/dp[2]
        J
    end
    jacobian_fit = curve_fit(model, jacobian_model, xdata, ydata, p0)
    @assert norm(jacobian_fit.param - [1.0, 2.0]) < 0.05
    @test jacobian_fit.converged

    # some example data
    yvars = 1e-6*rand(length(xdata))
    ydata = model(xdata, [1.0, 2.0]) + sqrt.(yvars) .* randn(length(xdata))

    fit = curve_fit(model, xdata, ydata, 1 ./ yvars, [0.5, 0.5])
    println("norm(fit.param - [1.0, 2.0]) < 0.05 ? ", norm(fit.param - [1.0, 2.0]))
    @assert norm(fit.param - [1.0, 2.0]) < 0.05
    @test fit.converged

    # can also get error estimates on the fit parameters
    errors = margin_error(fit, 0.1)
    println("norm(errors - [0.017, 0.075]) < 0.1 ?", norm(errors - [0.017, 0.075]))
    @assert norm(errors - [0.017, 0.075]) < 0.1

    # test with user-supplied jacobian and weights
    fit = curve_fit(model, jacobian_model, xdata, ydata, 1 ./ yvars, [0.5, 0.5])
    println("norm(fit.param - [1.0, 2.0]) < 0.05 ? ", norm(fit.param - [1.0, 2.0]))
    @assert norm(fit.param - [1.0, 2.0]) < 0.05
    @test fit.converged

    curve_fit(model, jacobian_model, xdata, ydata, 1 ./ yvars, [0.5, 0.5]; tau=0.0001)
end


text:
function process_arg!(
  srcs_used, helpers_ex, ::Type{InnerJoin}, ex, index, primary)::Void
    # NOTE: Only equijoins are supported for now
    @assert ex.head == :call && ex.args[1] == :(==)
    lhs, rhs = ex.args[2], ex.args[3]
    f_ex1, ai1 = build_f_ex!(srcs_used, lhs, index, primary)
    f_ex2, ai2 = build_f_ex!(srcs_used, rhs, index, primary)
    push!(helpers_ex.args,
          Expr(:call, InnerJoin, Expr(:tuple,
                                      Expr(:tuple, f_ex1, f_ex2),
                                      Expr(:tuple, ai1, ai2))))
    return
end

# function gen_helper_ex{H<:JoinHelper}(::Type{H}, ex)::Expr
#     arg_parameters = Set{Symbol}()
#     ex.head != :call && throw(ArgumentError())
#     ex.args[1] != :(==) && throw(ArgumentError("Non-equijoins not supported."))
#     lhs, rhs = ex.args[2], ex.args[3]
#     f1_ex, arg_fields1 = build_kernel_ex!(lhs, arg_parameters)
#     f2_ex, arg_fields2 = build_kernel_ex!(rhs, arg_parameters)
#     return Expr(:call, H, f1_ex, f2_ex, arg_fields1, arg_fields2)
# end


text:
export read_particles

function read_particles( filename, mesh :: Mesh )

    nbpart = countlines(filename)
    println( " nbpart   : ", nbpart )
    println( " filename : ", filename )

    nx, ny = mesh.nx, mesh.ny
    dx, dy = mesh.dx, mesh.dy
    dimx   = mesh.xmax - mesh.xmin
    dimy   = mesh.ymax - mesh.ymin

    weight    = (dimx * dimy) / nbpart
    particles = Particles( nbpart, weight)

    open(filename) do f

        for (k,line) in enumerate(eachline(f))
            str_values = split(line)
            ix = parse(Int32,   str_values[1])
            iy = parse(Int32,   str_values[2])
            dpx = parse(Float64, str_values[3])
            dpy = parse(Float64, str_values[4])
            particles.v[1,k] = parse(Float64, str_values[5])
            particles.v[2,k] = parse(Float64, str_values[6])
            particles.x[1,k] = (dpx+ix) * dx
            particles.x[2,k] = (dpy+iy) * dy
        end

    end

    particles

end


text:
mutable struct Actor
    image::String
    surface::Ptr{SDL2.Surface}
    position::Rect
    scale::Vector{Float64}
    angle::Float64
    alpha::UInt8
    data::Dict{Symbol, Any}
end

"""
`Actor(image::String)`

Creates an Actor with the image given, which must be located in the `images` subdirectory.
"""
function Actor(image::String; kv...)
    sf=image_surface(image)
    w, h = size(sf)
    a = Actor(image, sf, Rect(0, 0, Int(w), Int(h)), [1.0, 1.0], 0, 255, Dict{Symbol,Any}())

    for (k, v) in kv
        setproperty!(a, k, v)
    end
    return a
end

"""
    TextActor(text::String, font_name::String; font_size=24, color=Int[255,255,0,255])

Creates an actor with text rendered using font font_name. Font should be located in fonts directory.
"""
function TextActor(text::String, font_name::String; font_size=24, color=Int[255,255,0,255], kv...)
    font = SDL2.TTF_OpenFont(file_path(font_name, :fonts), font_size)
    sf = SDL2.TTF_RenderText_Blended(font, text, SDL2.Color(color...))
    w, h = size(sf)
    a = Actor(
        text,
        sf,
        Rect(0, 0, Int(w), Int(h)),
        [1.,1.],
        0,
        255,
        Dict{Symbol,Any}()
    )

    for (k, v) in kv
        setproperty!(a, k, v)
    end
    return a
end

function Base.setproperty!(s::Actor, p::Symbol, x)
    if p == :image
        sf = image_surface(x)
        setfield!(s, :surface, sf)
    elseif hasfield(Actor, p)
        setfield!(s, p, convert(fieldtype(Actor, p), x))
    else
        position = getfield(s, :position)
        v = getproperty(position, p)
        if v != nothing
            setproperty!(position, p, x)
        else
            getfield(s, :data)[p] = x
        end
    end
end

function Base.getproperty(s::Actor, p::Symbol)
    if hasfield(Actor, p)
        getfield(s, p)
    else
        position = getfield(s, :position)
        v = getproperty(position, p)
        if v != nothing
            return v
        else
            data = getfield(s, :data)
            if haskey(data, p)
                return data[p]
            else
                @warn "Unknown data $p requested from Actor($(s.image))"
                return nothing
            end
        end
    end
end

"""
`draw(a::Actor)`

Draws the Actor on-screen at its current position.
"""
function draw(a::Actor)
    texture = SDL2.CreateTextureFromSurface(game[].screen.renderer, a.surface)
    r=a.position
    w′=floor(r.w * a.scale[1])
    h′=floor(r.h * a.scale[2])

    if (a.alpha < 255)
        SDL2.SetTextureBlendMode(texture, SDL2.BLENDMODE_BLEND)
        SDL2.SetTextureAlphaMod(texture, a.alpha)
    end

    SDL2.RenderCopyEx(
        game[].screen.renderer,
        texture,
        C_NULL,
        Ref(SDL2.Rect(r.x, r.y, w′, h′)),
        a.angle,
        C_NULL,
        UInt32(0)
    )
    SDL2.DestroyTexture(texture)
end

"""
`angle(a1::Actor, a2::Actor)`

Angle between the horizontal and the line between two actors. Value returned is in degrees.
"""
function Base.angle(a::Actor, target::Actor)
    angle(a, a.pos...)
end

"""
`angle(a::Actor, xy::Tuple{Number, Number})`

Angle between the horizontal of the line between an actor and a point in space. Value returned is in degrees.
"""
Base.angle(a::Actor, txy::Tuple) = angle(a, txy[1], txy[2])

"""
`angle(a::Actor, x::Number, y::Number)`

Angle between the horizontal of the line between an actor and a point in space. Value returned is in degrees.
"""
function Base.angle(a::Actor, tx, ty)
    myx, myy = a.pos
    dx = tx - myx
    dy = myy - ty
    return deg2rad(atan(dy/dx))
end

"""
`distance(a1::Actor, a2::Actor)`

Distance in pixels between two actors.
"""
function distance(a::Actor, target::Actor)
    distance(a, target.pos...)
end

"""
`distance(a::Actor, x::Number, y::Number)`

Distance in pixels between an actor and a point in space
"""
function distance(a::Actor, tx, ty)
    myx, myy = a.pos
    dx = tx - myx
    dy = ty - myy
    return sqrt(dx * dx + dy * dy)
end

atan2(y, x) = pi - pi/2 * (1 + sign(x)) * (1 - sign(y^2)) - pi/4 * (2 + sign(x)) * sign(y) -
                            sign(x*y) * atan((abs(x) - abs(y)) / (abs(x) + abs(y)))

function Base.size(s::Ptr{SDL2.Surface})
    ss = unsafe_load(s)
    return (ss.w, ss.h)
end

"""
```
collide(a, x::Integer, y::Integer)
collide(a, xy::Tuple{Integer, Integer})
```

Checks if a (a game object) is colliding with a point.
"""
function collide(a, x::Integer, y::Integer)
    a=rect(a)
    return a.x <= x < (a.x + a.w) &&
        a.y <= y < (a.y + a.h)
end

collide(a, pos::Tuple) = collide(a, pos[1], pos[2])

"""
`collide(a, b)`

Checks if a and b (both game objects) are colliding.
"""
function collide(a, b)
    a=rect(a)
    b=rect(b)
    return a.x < b.x + b.w &&
        a.y < b.y + b.h &&
        a.x + a.w > b.x &&
        a.y + a.h > b.y

end

rect(a::Actor) = a.position


text:
using FITSIO, DataFrames, FITSIO

function sxtevtlist2df(hdu::TableHDU)
#    f=FITS(evtfile)

    ti=read(hdu,"TIME")
    ccdframe = read(hdu, "CCDFrame")
        x= read(hdu, "X")
        y= read(hdu, "Y")
        rawx= read(hdu, "RAWX")
        rawy= read(hdu, "RAWY")
        detx= read(hdu, "DETX")
        dety= read(hdu, "DETY")
        phas= read(hdu, "PHAS")
        stts= read(hdu, "STATUS")
        pha= read(hdu, "PHA")
        grade= read(hdu, "GRADE")
        pixsabove= read(hdu, "PixsAbove")
        pi_channel= read(hdu, "PI")
        ra= read(hdu, "RA")
        decl= read(hdu, "Dec")
        phaso= read(hdu, "PHASO")
        evtdf=DataFrame(TIME=ti,CCDFrame=ccdframe, X=x,Y=y,RAWX=rawx,RAWY=rawy,DETX=detx,DETY=dety,PHAS1=phas[1,:],PHAS2=phas[2,:],PHAS3=phas[3,:],PHAS4=phas[4,:],PHAS5=phas[5,:],PHAS6=phas[6,:],PHAS7=phas[7,:],PHAS8=phas[8,:],PHAS9=phas[9,:],STATUS=stts,PHA=pha,GRADE=grade,PixsAbove=pixsabove,PI=pi_channel,RA=ra,Dec=decl,PHASO1=phaso[1,:],PHASO2=phaso[2,:],PHASO3=phaso[3,:],PHASO4=phaso[4,:],PHASO5=phaso[5,:],PHASO6=phaso[6,:],PHASO7=phaso[7,:],PHASO8=phaso[8,:],PHASO9=phaso[9,:])
    return evtdf
end

function sxtbadpixlist2df(hdu::TableHDU)
#    f=FITS(evtfile)
        rawx= read(hdu, "RAWX")
        rawy= read(hdu, "RAWY")
        typ= read(hdu, "TYPE")
        yextnt= read(hdu, "YEXTENT")
        badflag= read(hdu, "BADFLAG")
        badpixdf=DataFrame(RAWX=rawx,RAWY=rawy,TYPE=typ,YEXTENT=yextnt,BADFLAG=badflag)
    return badpixdf
end

max_stars_repo_name: gulabd/SXTMerger

max_issues_repo_name: gulabd/SXTMerger


text:
struct vcovCluster <: vcov
    cluster::Matrix
    n_cluster_vars::Int64
    ssc::String
    # defaults to HC1 for now
    function vcovCluster(cluster::Matrix)
        n_cluster_vars = size(cluster, 2)
        ssc = "HC1"
        new(cluster, n_cluster_vars, ssc)
    end
end

"""
    findsinglecluster(cluster_var::Vector, X::Matrix, resid::Vector)
If there's a single cluster variable this is a much quicker implementation.
"""
function findsinglecluster(cluster_var::Vector, X::Matrix, resid::Vector)
    unique_clusters = unique(cluster_var)
    n_clusters = length(unique(cluster_var))
    B_clusters = Vector{Matrix{Float64}}(undef, n_clusters)
    N, K = size(X, 1), size(X, 2)
    dof_adjustment = (n_clusters/(n_clusters - 1))*((N-1)/(N-K))
    for i in 1:n_clusters
        cluster_index = findall(cluster_var .== unique_clusters[i])
        X_cl = X[cluster_index, :]
        resid_cl = resid[cluster_index, :]
        B_cl = (X_cl' * resid_cl) * (resid_cl' * X_cl)
        B_clusters[i] = B_cl*dof_adjustment
    end
    B_hat = reduce(+, B_clusters)
    return B_hat
end

I_R(R) = findall(R .== 1)
"""
    createRset(cluster_obj::vcovCluster)

Creates a set of D-length vectors where D corresponds to clustering dimension
size. R = {r: r_d in {0, 1}, d = 1, 2,..., D} (page 9 Cameron, Gelbach, Miller)

These are the permutations of all the clustering dimensions we want to compare
across.
"""
function createRset(cluster_obj::vcovCluster)
    n_cluster_vars = cluster_obj.n_cluster_vars
    ones_matrix = ones(n_cluster_vars, n_cluster_vars)
    triangular_matrix = LinearAlgebra.UpperTriangular(ones_matrix)
    # it's vectors all the way down
    perm_matrix = Vector{Vector{Vector{Int}}}(undef, n_cluster_vars)
    for i in 1:n_cluster_vars
        perm_matrix[i] = unique(permutations(triangular_matrix[i, :]))
    end
    # unnest vec of vec
    R_set = reduce(vcat, perm_matrix)
    return R_set
end

"""
    createBtilde(cluster_matrix, indicators)
Creates a clustered B_tilde matrix as described in Cameron, Gelbach,
and Miller (2008)
(http://cameron.econ.ucdavis.edu/research/CGM_twoway_ALL_13may2008.pdf)

For now small sample correction defaults to HC1 (which is naughty).
"""
function createBtilde(X::Matrix,
                      resid::Vector,
                      N::Int64,
                      cluster_matrix,
                      indicators)
    K = size(X, 2)
    # Initialise B_r for each r in the datset.
    B_R = Vector{Matrix}(undef, length(indicators))
    for I in 1:length(indicators)
        # Only keep cluster dimensions we care about
        subset_cluster_matrix = cluster_matrix[:, indicators[I]]
        # Find each unique cluster tuple/pair in the dataset
        unique_clusters = unique(subset_cluster_matrix, dims = 1)
        cluster_G_size = size(unique_clusters, 1)
        B_G = Vector{Matrix}(undef, cluster_G_size)
        for i in 1:size(unique_clusters, 1)
            # iterate through each unique cluster tuple and select only Xs and
            # resids corresponding to that cluster tuple. all(cond, dims = 2)
            # means we check that each ROW of the matrix meets the condition.
            cluster_selector = all(subset_cluster_matrix .== unique_clusters[i:i, :], dims = 2)[:]
            X_g = X[cluster_selector, :]
            resid_g = resid[cluster_selector]
            # Create B_G(I_r, ij) using these Xs and resids
            B_G[i] = (X_g' * resid_g) * (resid_g' * X_g)
        end
        # Create odd even inclusion adjustment
        D = length(indicators[I])
        small_sample_correction = (cluster_G_size/(cluster_G_size - 1))*((N - 1)/(N-K))
        sum_adjustment = (-1)^(D+1)*small_sample_correction
        B_R[I] = sum_adjustment*reduce(+, B_G)
    end
    B_tilde = sum(B_R)
    return B_tilde
end

"""
    mat_posdef_fix(X::Matrix; tol = 1e-10)
Make matrix posdef if we have negative variances on the diagonal.

Proposed by Cameron, Gelbach, Miller (2011).
"""
function mat_posdef_fix(X::Matrix; tol = 1e-10)
    if any(diag(X) .< tol)
        e_vals, e_vecs = eigen(Symmetric(X))
        e_vals[e_vals .<= tol] .= tol
        X = e_vecs * Diagonal(e_vals) * e_vecs'
    end
    return X
end

"""
Clustered Standard errors inference.
"""
function inference(N::Int,
                   K::Int,
                   resid::Vector,
                   β::Vector,
                   XX_inv::Matrix,
                   X::Matrix,
                   vcov::vcovCluster)
    cluster_matrix = vcov.cluster

    n_cluster_vars = size(cluster_matrix, 2)
    if n_cluster_vars == 1
        B_tilde = findsinglecluster(cluster_matrix[:], X, resid)
    else
        Rset = createRset(vcov)
        indicators = I_R.(Rset)
        B_tilde = createBtilde(X, resid, N, cluster_matrix, indicators)
    end

    vcov_matrix = XX_inv * B_tilde * XX_inv
    vcov_matrix = mat_posdef_fix(vcov_matrix)
    se = sqrt.(diag(vcov_matrix))

    pval = 2 .* cdf.(TDist(N - K), -abs.(β ./ se))
    σ_sq = sum(resid.^2) / (N - K)
    return se, pval, σ_sq, vcov_matrix
end

max_stars_repo_name: EdJeeOnGitHub/jeeves

max_issues_repo_name: EdJeeOnGitHub/jeeves

max_forks_repo_name: EdJeeOnGitHub/jeeves


text:
using TOML, Artifacts, Base.BinaryPlatforms
include("./platform_augmentation.jl")
artifacts_toml = joinpath(dirname(@__DIR__), "Artifacts.toml")

# Get "target triplet" from ARGS, if given (defaulting to the host triplet otherwise)
target_triplet = get(ARGS, 1, Base.BinaryPlatforms.host_triplet())

# Augment this platform object with any special tags we require
platform = augment_platform!(HostPlatform(parse(Platform, target_triplet)))

# Select all downloadable artifacts that match that platform
artifacts = select_downloadable_artifacts(artifacts_toml; platform)

#Output the result to `stdout` as a TOML dictionary
TOML.print(stdout, artifacts)


text:
"""
A material for the illumination model (e.g. Phong illumination).
"""
abstract type Material end

"""
Data structure for a mesh material that is used to describe the light components of a [Phong reflection](https://en.wikipedia.org/wiki/Phong_reflection_model)
type model. All data is stored as RGBα for Red, Green, Blue and transparency.
"""
struct Phong{T,S<:Colorant} <: Material
    emission::S
    ambient::S
    diffuse::S
    specular::S
    shininess::T
end

"""
    material_to_opf_string(material::Phong)
    material_to_opf_string(material::Colorant)

Format a material into a Dict for OPF writting.
"""
function material_to_opf_string(material::Phong)
    Dict(
        :emission => colorant_to_string(material.emission),
        :ambient => colorant_to_string(material.ambient),
        :diffuse => colorant_to_string(material.diffuse),
        :specular => colorant_to_string(material.specular),
        :shininess => string(material.shininess)
    )
end

function material_to_opf_string(material::Colorant)
    # Here we use the same color for all Phong parameters
    Dict(
        :emission => join(Float64[0.0, 0.0, 0.0, 0.0], "\t"),
        :ambient => colorant_to_string(material),
        :diffuse => colorant_to_string(material),
        :specular => colorant_to_string(material),
        :shininess => "1.0"
    )
end

"""
    colorant_to_string(x)

Parse a geometry material for OPF writing.
"""
function colorant_to_string(x::T) where {T<:RGBA}
    join(Float64[x.r, x.g, x.b, x.alpha], "\t")
end

function colorant_to_string(x::T) where {T<:RGB}
    join(Float64[x.r, x.g, x.b, 1.0], "\t")
end

"""

    RefMesh(
        name::S
        mesh::SimpleMesh
        normals::N
        texture_coords::T
        material::M
        taper::Bool
    )

    RefMesh(name, mesh, material = RGB(220 / 255, 220 / 255, 220 / 255))

RefMesh type. Stores all information about a Mesh:
- `name::S`: the mesh name
- `mesh::SimpleMesh`: the actual mesh information -> points and topology
- `normals::Vector{Float64}`: the normals, given as a vector of x1,y1,z1,x2,y2,z2...
- `texture_coords::Vector{Float64}`: the texture coordinates (not used yet), idem, a vector
- `material::M`: the material, used to set the shading
- `taper::Bool`: `true` if tapering is enabled

The reference meshes are then transformed on each node of the MTG using a transformation matrix
to match the actual mesh.
"""
struct RefMesh{S<:Union{String,SubString},ME<:Meshes.SimpleMesh,M<:Union{Material,Colorant},N<:SVector,T<:Union{SVector,Nothing}}
    name::S
    mesh::ME
    normals::N
    texture_coords::T
    material::M
    taper::Bool
end

#! Make a method that computes the normals and texture_coords from the mesh

function RefMesh(name, mesh, material=RGB(220 / 255, 220 / 255, 220 / 255))
    RefMesh(
        name,
        mesh,
        normals_vertex(mesh),
        nothing,
        material,
        false
    )
end

"""
RefMeshes type. Data base that stores all [`RefMesh`](@ref) in an MTG. Usually stored in the
`:ref_meshes` attribute of the root node.
"""
mutable struct RefMeshes
    meshes::Vector{RefMesh}
end

Base.names(m::RefMeshes) = [i.name for i in m.meshes]

"""
    geometry(
        ref_mesh::M
        ref_mesh_index::Union{Int,Nothing}
        transformation::T
        dUp::S
        dDwn::S
        mesh::Union{SimpleMesh,Nothing}
    )

A Node geometry with the reference mesh, its transformation matrix and optionnally the
index of the reference mesh in the reference meshes data base (see notes) and the resulting
mesh (optional to save memory).

# Note

The ref_mesh usually points to a [`RefMesh`](@ref) stored in the `:ref_meshes` attribute of the
root node of the MTG.

Although optional, storing the index of the reference mesh (`ref_mesh_index`) in the database allows a faster
writing of the MTG as an OPF to disk.

If no transformation matrix is needed, you can use `I` from the Linear Algebra package (lazy)

The `transformation` field should be a `CoordinateTransformations.jl`'s transformation. In case
no transformation is needed, use `IdentityTransformation()`. If you already have the
transformation matrix, you can pass it to `LinearMap()`.
"""
mutable struct geometry{M<:RefMesh,S}
    ref_mesh::M
    ref_mesh_index::Union{Int,Nothing}
    transformation::Transformation #! replace by concrete types ?
    dUp::S
    dDwn::S
    mesh::Union{Meshes.SimpleMesh,Nothing}
end


text:
using Quiqbox

nuc = ["H", "H"]

nucCoords = [[-0.7,0.0,0.0], [0.7,0.0,0.0]]

grid = GridBox(1, 3.0)

gf1 = GaussFunc(0.7,1)

bs = genBasisFunc.(grid.box, Ref([gf1]))

pars = uniqueParams!(bs, filterMapping=true)

parsPartial = [pars[1], pars[4]]

Es, pars, grads = optimizeParams!(bs, parsPartial, nuc, nucCoords, maxSteps=20)

# # You can also use more advanced optimizers from other packages.
# using Flux # First do `Pkg.add("Flux")` if you haven't installed the package.
# using Flux.Optimise: update!
# optimizer = AMSGrad(0.001)
# GDm = (prs, grad) -> update!(optimizer, prs, grad)
# optimizeParams!(bs, parsPartial, nuc, nucCoords; GDmethod=GDm)


text:
@testset "1014.best-sightseeing-pair.jl" begin
    @test max_score_sightseeing_pair([8, 1, 5, 2, 6]) == 11
    @test max_score_sightseeing_pair([24, 9, 13, 15, 8, 4, 19, 30, 8, 3, 18, 30, 14, 15, 8,
                                      12, 11, 14, 7, 20, 28, 4, 11, 7, 8, 2, 1, 7, 2, 3]) == 56
end


text:
default_acceleration = CPU1()
default_solver_test = Feagin14()
default_solver_train = Tsit5(; thread=OrdinaryDiffEq.True())
default_sensealg = InterpolatingAdjoint(
    ;
    autodiff=true,
    chunk_size=0,
    autojacvec=ZygoteVJP(),
)
default_optimizer = Dict(
    FluxOptApp => Flux.AMSGrad(),
    OptimOptApp => BFGS(),
    SciMLOptApp => Optimisers.AMSGrad(),
)


text:
"""
http://adventofcode.com/2017/day/16
"""

function danceMove(programs::String, move::SubString{String})
    regex = r"(s|x|p)([0-9a-z]+)/?([0-9a-z]+)?"
    moveType,pos1,pos2 = match(regex, move).captures
    if moveType == "s"
        output = string(programs[end-parse(Int,pos1)+1:end], programs[1:end-parse(Int,pos1)])
    elseif moveType == "x"
        i,j = sort([parse(Int, pos1)+1, parse(Int, pos2)+1])
        output = string(programs[1:i-1], programs[j], programs[i+1:j-1], programs[i], programs[j+1:end])
    elseif moveType == "p"
        i,j = sort([i for (i,letter) in enumerate(programs) if string(letter) in [pos1,pos2]])
        output = string(programs[1:i-1], programs[j], programs[i+1:j-1], programs[i], programs[j+1:end])
    end
    return output
end

function makeDance(programs::String, moveList::Array{SubString{String},1})
    output = programs[1:end]
    for i in 1:length(moveList)
        move = moveList[i]
        output = danceMove(output, move)
    end
    return output
end

file = f = open("day16_input.txt")
s = readstring(f)
moves = split(s, ",")

programs = "abcdefghijklmnop"
programs_0 = "abcdefghijklmnop"

for i in 1:1000000000%42
    if programs == programs_0
        println(i)
    end
    programs = makeDance(programs, moves)
end

println(programs)


text:
"""
    binary_search(sorted::AbstractVector{T}, value::T)::Int where {T<:Real}

Perfom a binary search to return the index i of an sorted array
such that sorted[i-1] < value <= sorted[i]

# Examples
```jldoctest
julia> MonteCarloX.binary_search([1.,2.,3.,4.],2.5)
3
```
```jldoctest
julia> MonteCarloX.binary_search([1,2,3,4],2)
2
```
"""
function binary_search(sorted::AbstractVector{T}, value::T)::Int where {T <: Real}
    # catch lower-bound case that cannot be reached by binary search
    id = 1
    if value > sorted[1]
        index_l = 1
        index_r = length(sorted)
        while index_l < index_r - 1
            # this should be fine because build for integers!
            index_m = fld(index_l + index_r, 2)
            if sorted[index_m] < value
                index_l = index_m
            else
                index_r = index_m
            end
        end
        id = index_r
    end
    return id
end

"""
    log_sum(a::T,b::T)

Return result of logarithmic sum ``c = \\ln(A+B) = a + \\ln(1+e^{|b-a|})`` where ``C =
e^c = A+B = e^a + e^b``.

This is useful for sums that involve elements that span multiple orders of
magnitude, e.g., the partition sum that is required as normalization factor
during reweighting.

# Examples
```jldoctest
julia> exp(MonteCarloX.log_sum(log(2.), log(3.)))
5.000000000000001

```

"""
function log_sum(a::T, b::T)::T where T <: AbstractFloat
    if b < a
        return a + log(1 + exp(b - a))
    else
        return b + log(1 + exp(a - b))
    end
end

"""
    random_element(list_probabilities::Vector{T},rng::AbstractRNG)::Int where T<:AbstractFloat

Pick an index with probability defined by `list_probability` (which needs to be normalized).

#Remark
Deprecated for use of StatsBase.sample

# Examples
```jldoctest
julia> using Random

julia> rng = MersenneTwister(1000);

julia> MonteCarloX.random_element([0.1,0.2,0.3,0.4],rng)
4
julia> MonteCarloX.random_element([0.1,0.2,0.3,0.4],rng)
4
julia> MonteCarloX.random_element([0.1,0.2,0.3,0.4],rng)
3
julia> MonteCarloX.random_element([0.1,0.2,0.3,0.4],rng)
4
```
"""
function random_element(list_probabilities::Vector{Float64}, rng::AbstractRNG)::Int
    theta = rand(rng) * sum(list_probabilities)

    id = 1
    cumulated_prob = list_probabilities[id]
    while cumulated_prob < theta
        id += 1
        @inbounds cumulated_prob += list_probabilities[id]
    end

    return id
end

# TODO: iterate over n-dimensional StatsBase.Histogram
# for start, CartesianIndices(hist.weights) can give you a way to iterate over
# weights, and the cartesian index is your bin, if you want to convert the bin
# number back to edge value, you just need to access the corresponding element
# in hist.edges

# does this need performance boost? - maybe this is faseter when explicit for 1D
Base.getindex(h::AbstractHistogram{T,1}, x::Real) where {T} = getindex(h, (x,))

function Base.getindex(h::Histogram{T,N}, xs::NTuple{N,Real}) where {T,N}
    idx = StatsBase.binindex(h, xs)
    if checkbounds(Bool, h.weights, idx...)
        return @inbounds h.weights[idx...]
    else
        return missing
    end
end

Base.setindex!(h::AbstractHistogram{T,1}, value::Real, x::Real) where {T} = setindex!(h, value, (x,))

function Base.setindex!(h::Histogram{T,N}, value::Real, xs::NTuple{N,Real}) where {T,N}
    h.isdensity && error("Density histogram must have float-type weights")
    idx = StatsBase.binindex(h, xs)
    if checkbounds(Bool, h.weights, idx...)
        @inbounds h.weights[idx...] = value
    else
        return missing
    end
end


text:
abstract type Shape{N,T} <: FieldVector{N,T} end

struct Cube{T} <: Shape{3,T}
    a::T
end

struct Sphere{T} <: Shape{3,T}
    r::T
end

struct SphericalCap{T} <: Shape{3,T}
    a::T
    c::T
    function SphericalCap{T}(a, c) where {T}
        if a < 0 || c ≤ 0
            throw(ArgumentError("All lengths should be positive."))
        end
        return new(a, c)
    end
end

for shape in (:Ellipsoid, :EllipticCylinder, :Parallelepiped)
    @eval begin
        struct $shape{T} <: Shape{3,T}
            a::T
            b::T
            c::T
            function $shape{T}(a, b, c) where {T}
                if a ≤ 0 || b ≤ 0 || c ≤ 0
                    throw(ArgumentError("All lengths should be positive."))
                end
                return new(a, b, c)
            end
        end
    end
end

struct Cylinder{T} <: Shape{3,T}
    r::T
    c::T
    function Cylinder{T}(r, c) where {T}
        if r ≤ 0 || c ≤ 0
            throw(ArgumentError("All lengths should be positive."))
        end
        return new(r, c)
    end
end

struct HollowCylinder{T} <: Shape{3,T}
    R::T
    r::T
    c::T
    function HollowCylinder{T}(R, r, c) where {T}
        if r < 0 || c ≤ 0
            throw(ArgumentError("All lengths should be positive."))
        elseif R < r
            throw(ArgumentError("Outer radius `R` should be larger than `r`."))
        end
        return new(R, r, c)
    end
end

struct SquarePyramid{T,S} <: Shape{3,T}
    a::T
    b::T
    m::S
    function SquarePyramid{T,S}(a, b, m) where {T,S}
        if a ≤ 0 || b ≤ 0
            throw(ArgumentError("All lengths should be positive."))
        end
        @assert m == 2b / a
        return new(a, b, m)
    end
end

struct RectangularPyramid{T,S} <: Shape{3,T}
    a::T
    b::T
    c::T
    ma::S
    mb::S
    function RectangularPyramid{T,S}(a, b, c, ma, mb) where {T,S}
        if a ≤ 0 || b ≤ 0 || c ≤ 0
            throw(ArgumentError("All lengths should be positive."))
        end
        @assert ma == 2c / a
        @assert mb == 2c / b
        return new(a, b, c, ma, mb)
    end
end

struct TriangularToroid{T} <: Shape{3,T}
    r::T
    b::T
    c::T
    function TriangularToroid{T}(r, b, c) where {T}
        if r ≤ 0 || b ≤ 0 || c ≤ 0
            throw(ArgumentError("All lengths should be positive."))
        elseif r < b
            throw(ArgumentError("The radius must be equal or higher than " *
                                "the triangle base semi-length"))
        end
        return new(r, b, c)
    end
end

struct TruncatedSquarePyramid{T,R,S} <: Shape{3,T}
    a::T
    b::T
    r::R
    m::S
    function TruncatedSquarePyramid{T,R,S}(a, b, r, m) where {T,R,S}
        if a ≤ 0 || b ≤ 0
            throw(ArgumentError("All semi-lengths should be positive."))
        end
        @assert 0 < r ≤ 1
        @assert m == 2b / a
        return new(a, b, r, m)
    end
end

struct Ring{T} <: Shape{3,T}
    R::T
    a::T
    b::T
end

struct Torus{T} <: Shape{3,T}
    R::T
    r::T
end

struct Vec{T} <: FieldVector{3,T}
    x::T
    y::T
    z::T
end


text:
# returns the variable (or matrix), lagged p times,
# with the first p rows filled with ones (to avoid divide errors)
# remember to drop those rows before doing analysis
function lag(x::Array{Float64,2},p::Int64)
	n,k = size(x)
	lagged_x = [ones(p,k); x[1:n-p,:]]
end

function lag(x::Array{Float64,1},p::Int64)
	n = size(x,1)
	lagged_x = [ones(p); x[1:n-p]]
end

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
############################ stationary kernels ################################
# using LinearAlgebraExtensions: difference
# notation:
# x, y inputs
# τ = x-y difference of inputs
# r = norm(x-y) norm of τ
# r² = r^2
(k::StationaryKernel)(x, y) = k(difference(x, y)) # if the two argument signature is not defined, default to stationary
(k::IsotropicKernel)(x, y) = k(euclidean2(x, y)) # if the two argument signature is not defined, default to isotropic
(k::IsotropicKernel)(τ) = k(sum(abs2, τ)) # if only the scalar argument form is defined, must be isotropic

############################# constant kernel ##################################
# can be used to rescale existing kernels
# IDEA: Allow Matrix-valued constant
struct Constant{T} <: IsotropicKernel{T}
    c::T
    function Constant(c, check::Bool = true)
        if check && !ispsd(c)
            throw(DomainError("Constant is not positive semi-definite: $c"))
        end
        new{typeof(c)}(c)
    end
end
# isisotropic(::Constant) = true
# ismercer(k::Constant) = ispsd(k.c)
# Constant(c) = Constant{typeof(c)}(c)

# should type of constant field and r agree? what promotion is necessary?
# do we need the isotropic/ stationary evaluation, if we overwrite the mercer one?
(k::Constant)() = k.c
(k::Constant)(r²) = k.c # stationary / isotropic
(k::Constant)(x, y) = k.c # mercer

gramian(k::Constant, x::AbstractVector, y::AbstractVector) = Fill(k.c, length(x), length(y))

#################### standard exponentiated quadratic kernel ###################
struct ExponentiatedQuadratic{T} <: IsotropicKernel{T} end
const EQ{T} = ExponentiatedQuadratic{T}
@functor EQ
EQ() = EQ{Union{}}() # defaults to "bottom" type since it doesn't have any parameters
(k::EQ)(r²::Number) = exp(-r²/2)

########################## rational quadratic kernel ###########################
struct RationalQuadratic{T} <: IsotropicKernel{T}
    α::T # relative weighting of small and large length scales
    RationalQuadratic{T}(α) where T = (0 < α) ? new(α) : throw(DomainError("α not positive"))
end
const RQ = RationalQuadratic
@functor RQ
RQ(α::Real) = RQ{typeof(α)}(α)

(k::RQ)(r²::Number) = (1 + r² / (2*k.α))^-k.α

########################### exponential kernel #################################
struct Exponential{T} <: IsotropicKernel{T} end
const Exp = Exponential
Exp() = Exp{Union{}}()

(k::Exp)(r²::Number) = exp(-sqrt(r²))

############################ γ-exponential kernel ##############################
struct GammaExponential{T<:Real} <: IsotropicKernel{T}
    γ::T
    GammaExponential{T}(γ) where {T} = (0 ≤ γ ≤ 2) ? new(γ) : throw(DomainError("γ not in [0,2]"))
end
const γExp = GammaExponential
@functor γExp
γExp(γ::T) where T = γExp{T}(γ)

(k::γExp)(r²::Number) = exp(-r²^(k.γ/2) / 2)

########################### white noise kernel #################################
struct Delta{T} <: IsotropicKernel{T} end
@functor Delta
const δ = Delta
δ() = δ{Union{}}()

(k::δ)(r²) = all(iszero, r²) ? one(eltype(r²)) : zero(eltype(r²))
function (k::δ)(x, y)
    T = promote_type(eltype(x), eltype(y))
    (x == y) ? one(T) : zero(T) # IDEA: if we checked (x === y) could incorporate noise variance for vector inputs -> EquivDelta?
end
############################ Matern kernel #####################################
# IDEA: use rational types to dispatch to MaternP evaluation, i.e. 5//2 -> MaternP(3)
# seems k/2 are representable exactly in floating point?
struct Matern{T} <: IsotropicKernel{T}
    ν::T
    Matern{T}(ν) where T = (0 < ν) ? new(ν) : throw(DomainError("ν = $ν is negative"))
end
@functor Matern
Matern(ν::T) where {T} = Matern{T}(ν)

(k::Matern)(r²::Number) = Matern(r², k.ν)

# NOTE: this implementation only guarantees second order differentiability
# for higher orders, need to get higher order taylor expansion
@inline function Matern(r²::Number, ν::Real)
    ε = eps(typeof(r²))
    taylor_bound = (2 < ν) ? ε^(1/2) : ((1 < ν) ? ε : zero(ε))
    if r² < taylor_bound
        y = one(promote_type(typeof(r²), typeof(ν)))
        if ν > 1
            y += ν / (2*(1-ν)) * r² # first order
        end
        if ν > 2
            y += ν^2 / (8*(2 - 3ν + ν^2)) * (r²)^2 # second order
        end
        return y
    else
        r = sqrt(2ν*r²)
        2^(1-ν) / gamma(ν) * adbesselkxv(ν, r) # is equal to r^ν * besselk(ν, r)
    end
end

################# Matern kernel with ν = p + 1/2 where p ∈ ℕ ###################
struct MaternP{T, DT, CT} <: IsotropicKernel{T}
    p::Int
    derivatives::DT
    coefficients::CT
end

function MaternP(p::Int)
    0 > p && throw(DomainError("p = $p is negative"))
    d = MaternP_derivatives_at_zero(p)
    d = float(d)
    c = MaternP_coefficients(p)
    c = float(c)
    MaternP{Union{}, typeof(d), typeof(c)}(p, d, c)
end
MaternP(k::Matern) = MaternP(floor(Int, k.ν)) # project Matern to closest MaternP

(k::MaternP)(r²::Integer) = k(float(r²))
function (k::MaternP)(r²::Number)
    p = k.p
    r²_constant = (r² isa Taylor1 ? r².coeffs[1] : r²)
    taylor_bound = eps(typeof(r²_constant))^(1/p)
    use_taylor = r²_constant < taylor_bound
    if use_taylor # taylor_bound # around zero, use taylor expansion in r² to avoid singularity in derivative
        y = one(r²)
        r²i = r²
        for i in 1:p # iterating over r²^i allows for automatic differentiation (even though for a regular float, it would all be zero)
            y += k.derivatives[i] * r²i / factorial(i)
            r²i *= r²
        end
        return y
    else
        y = zero(r²)
        r = sqrt((2p+1)*r²)
        ri = one(r)
        for i in 1:p
            y += k.coefficients[i] * ri # (2r)^(p-i)
            ri *= 2r
        end
        y += ri # coefficient of (2r)^p is 1
        y *= exp(-r) / (factorial(2p) ÷ factorial(p))
    end
end

# naïve implementation, does not allow for differentiability at zero
# and recomputes coefficients based on factorials
@inline function MaternP(r²::Number, p::Int)
    val = zero(r²)
    r = sqrt((2p+1)*r²)
    for i in 0:p # IDEA: could table factorials and generate powers by iteration
        val += (factorial(p+i) ÷ (factorial(p-i) * factorial(i))) * (2r)^(p-i)
    end
    val *= exp(-r) / (factorial(2p) ÷ factorial(p))
end

# computes derivatives of k(r²) w.r.t. r² at zero
function MaternP_derivatives_at_zero(p::Int)
    r² = symbols(:r²) # uses SymEngine
    kr = MaternP(r², p)
    d = zeros(Rational{Int64}, p)
    for i in 1:p
        dkr = expand(diff(kr, r²))
        d[i] = dkr(0)
        kr = dkr
    end
    return d
end

function MaternP_coefficients(p::Int)
    coefficients = zeros(Rational{typeof(p)}, p)
    factorial_p = factorial(p)
    for i in 1:p
        coefficients[i] = (binomial(p, i) * (factorial(p+i) ÷ factorial_p))
    end
    return reverse(coefficients)
end

########################### cosine kernel ######################################
# interesting because it allows negative co-variances
# it is a valid stationary kernel,
# because it is the inverse Fourier transform of point measure at μ (delta distribution)
struct CosineKernel{T, V<:Union{T, AbstractVector{T}}} <: StationaryKernel{T}
    c::V
end
@functor CosineKernel
const Cosine = CosineKernel
const Cos = Cosine

# IDEA: trig-identity -> low-rank gramian
# NOTE: this is the only stationary non-isotropic kernel so far
input_trait(::CosineKernel) = StationaryLinearFunctionalInput() # dependent on dot(c, τ)
(k::CosineKernel)(c_dot_τ::Real) = cos(2π * c_dot_τ)
function (k::CosineKernel)(x, y)
    τ = difference(x, y)
    k(dot(k.c, τ))
end

####################### spectral mixture kernel ################################
# can be seen as product kernel of Constant, Cosine, ExponentiatedQuadratic
Spectral(w::Real, μ, l) = prod((w, Cosine(μ), ARD(EQ(), l))) # 2π^2/σ.^2)
SpectralMixture(w::AbstractVector, μ, l) = sum(Spectral.(w, μ, l))
const SM = SpectralMixture

############################ Cauchy Kernel #####################################
# there is something else in the literature with the same name ...
struct Cauchy{T} <: IsotropicKernel{T} end
@functor Cauchy
Cauchy() = Cauchy{Union{}}()
(k::Cauchy)(r²::Number) = inv(1+r²) # π is not necessary, we are not normalizing

# for spectroscopy
PseudoVoigt(α) = α*EQ() + (1-α)*Cauchy()

###################### Inverse Multi-Quadratic  ################################
# seems formally similar to Cauchy, Cauchy is equal to power of IMQ
struct InverseMultiQuadratic{T} <: IsotropicKernel{T}
    c::T
end
@functor InverseMultiQuadratic
(k::InverseMultiQuadratic)(r²::Number) = 1/√(r² + k.c^2)


text:
import Base.@deprecate
import Base.depwarn

# ci -> confint
@deprecate ci(bs::BootstrapSample, cim::ConfIntMethod) confint(bs, cim)
@deprecate ci(bs::BootstrapSample, cim::ConfIntMethod, i::Int) confint(bs, cim, i)
@deprecate ci(bs::BootstrapSample, sd::AbstractVector, cim::ConfIntMethod, i::Int) confint(bs, sd, cim, i)

# se -> stderror
@deprecate se(x::AbstractVector) stderror(x)
@deprecate se(bs::BootstrapSample) stderror(bs)
@deprecate se(bs::BootstrapSample, idx::Int) stderror(bs, idx)

# functions before data
@deprecate bootstrap(data, statistic::Function, sampling::BootstrapSampling) bootstrap(statistic, data, sampling)
@deprecate bootstrap(data, statistic::Function, model, sampling::BootstrapSampling) bootstrap(statistic, data, model, sampling)


text:
function solve()
    s = readline()
    ret = 0
    c = 0
    for w in s
        if w == 'R'
            c += 1
        else
            c = 0
        end
        ret = max(ret, c)
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
include("line_scatter.jl")

function subplots1()
    p1 = linescatter1()
    p2 = linescatter2()
    p = [p1 p2]
    p
end

function subplots2()
    p1 = linescatter1()
    p2 = linescatter2()
    p = [p1, p2]
    p
end

function subplots3()
    p1 = linescatter6()
    p2 = linescatter2()
    p3 = linescatter3()
    p4 = batman()
    p = [p1 p2; p3 p4]
    p.plot.layout["showlegend"] = false
    p.plot.layout["width"] = 1000
    p.plot.layout["height"] = 600
    p
end


text:
const STORAGE_API = Ref("https://document-storage-production-dot-remarkable-production.appspot.com")
const SERVICE_DISCOVERY_API = "https://service-manager-production-dot-remarkable-production.appspot.com"
const STORAGE_URL = "/service/json/1/document-storage"
const AUTH_API = "https://my.remarkable.com"
const UPDATE_STATUS = "upload/update-status"
const NEW_DEVICE = "/token/json/2/device/new"
const NEW_TOKEN = "/token/json/2/user/new"
const ITEM_LIST = "/document-storage/json/2/docs"
const ITEM = "/document-storage/json/2/"


text:
module module_RelNeq

using Hecke

struct RelNeq
  k::AnticNumberField
  K::AnticNumberField
  Kk::Hecke.NfRel{nf_elem}
  m_k_K::Map
  m_Kk_K::Map
  function RelNeq(k::AnticNumberField, Kk::Hecke.NfRel{nf_elem})
    k = base_ring(Kk)
    K, m_K_Kk = absolute_simple_field(Kk)
    m1 = inv(m_K_Kk)
    return new(k, K, Kk, m1, restrict(m1, k))
  end

  function RelNeq(k::AnticNumberField, K::AnticNumberField)
    kt, t = PolynomialRing(k, cached = false)
    fl, mp = Hecke.issubfield(k, K)
    Qt = parent(K.pol)
    h = gcd(gen(k) - evaluate(Qt(mp(gen(k))), t), evaluate(K.pol, t))
    Kk, _ = number_field(h)
    return new(k, K, Kk, mp, mp)
  end

  function RelNeq(m::Map)
    return RelNeq(domain(m), codomain(m))
  end
end

#= Via Lorenz... Let K/k be finite

  {N(I) = 1}
  ----------   < Cl_(Z_K * m)
  {N(a) = 1}

  For m the lorenz_module for n=K:k in Z_k: any unit u = 1 (m) is an n-th power.
  If I = J in Cl_(Z_K m) => I = a J and a = 1 (Z_K m). Furthermore, since N(I) = N(J) = 1Zk (the trivial ideal)
  N(a) is a unit. Since a = 1 (Z_K m) we have N(a) = 1 (m) so N(a) = eps^n, thus N(a/eps) = 1 and
  I = J in the norm-1-group.
  Since the map (injection) is only well defined on ideals coprime to m, we don't get everything.
  We're missing the contribution of finitely many explicit ideals, hence this can be fixed.

  In the case of k = Q, this amounts to the strict (narrow) class group (I think).

  This is easier to handle than the version below (free of lorenz_module).

  Next Q: Let N^1 be the subgroup of Cl_M defined above. Then

           X_M  the ray class field mod M = Z_K m
            |
            |
            X   the class field for N^1, thus the norm group is Cl_M/N^1
         Y  |   the norm group of Y is the norm of X_M over k.
         |  |   Possibly X = KY? To neat, but maybe a central/ ...? extension
         |  K
         | /
         |/
         k

 Now, finally the question: has X any (useful) Galois theoretic characterisation such as
   the minmal/ maximal field ...
   the compositum with ...
   is a version of the genus field ...
 Can one use this unknown characterisation to compute X (the degree of X) hence get a method to
 provably compute N^1?
 So far, to compute N^1: just try small prime ideals until the computation stablises. Then hope.
   Alternatively: use Jan Albert's thesis to get bound on the generation...

 For K/k normal this should reduce(?) to Scholz' paper (or Jehne's version) on Zahlknoten/ number knots
 In this case N(I) = 1 iff I = prod A^(1-sigma a) thus it can be computed from the Galois action on
 ideal (classes), then Cl/N^1 is fixed pointwise by Gal(K/k) (trivial action) (in fact this should be the
 largest/ smallest such quotient)
 hence the group extension (exact sequence)

   1 -> Cl/N^1 -> Gal(X/k) -> Gal(K/k) -> 1

 is split(?) direct(?) (nice?) (special?)

 My memory of Scholz is that there is a second field lurking around...
=#
function norm_1_subgroup(A::RelNeq)
  k = A.k
  K = A.K
  Kk = A.Kk

  d = lcm([numerator(discriminant(k)), numerator(discriminant(K)), numerator(norm(discriminant(Kk)))])

  @show I = Hecke.lorenz_module(k, degree(Kk), containing = discriminant(maximal_order(Kk)))
  Hecke.assure_2_normal(I)

  I_K = ideal(maximal_order(K), I.gen_one, maximal_order(K)(A.m_k_K(I.gen_two.elem_in_nf)))
  r, mr = ray_class_group(I_K, real_places(K))

  q, mq = quo(r, elem_type(r)[])

  S = PrimesSet(1, -1)
  gens = Set{NfOrdIdl}()
  gg = []

  max_stable = 2*ngens(r)
  stable = max_stable

  for p = S
    @show d, p, d%p
    if minimum(I) % p == 0
      continue
    end

    lp = prime_decomposition(maximal_order(k), p)
    for P = lp
      if d % p == 0
        @show "expensive", p
        lP = collect(factor(IdealSet(ZK)(A.m_k_K, P[1])))
      else
        lP = Hecke.prime_decomposition_nonindex(A.m_k_K, P[1])
      end
      f = [fmpz(div(degree(Q[1]), degree(P[1]))) for Q = lP]
      m = matrix(FlintZZ, 1, length(f), f)
      r, n = nullspace(m)

      decom = [mq(mr\Q[1]) for Q = lP]
      for i=1:r
        a = sum(n[j,i] * decom[j] for j = 1:length(decom))
        if iszero(a)
          stable -= 1
          continue
        end
        stable = max_stable

        q, nq = quo(q, [a])
        mq = mq*nq
        decom = [nq(x) for x = decom]
        push!(gens, P[1])
        push!(gg, FacElem(Dict((lP[j][1], n[j, i]) for j=1:length(decom) if n[j,i] != 0)))
      end
    end
    if stable <= 0
      break
    end
  end

  return mr, mq, gens, gg
end

#= The idea
  we have
                        1

                        |
                        V

        U_K        {N(a) in U_k}      {N(I) = 1}
1 -> ---------- -> -------------  ->  ----------  -> Cl
     {N(u) = 1}      {N(a) = 1}       {N(a) = 1}

                        | N
                        V

                       U_k

 So what we want is a group extension of a sub-group of a quotient of U_k
 by a subgroup of Cl:

                        1

                        |
                        V

                {N(a) in U_k}      {N(I) = 1}
      1  ->   ----------------  ->  ----------  -> Cl
              N(U_K){N(a) = 1}       {N(a) = 1}

                        | N
                        V

                       U_k
                      -----
                      N(U_K)

=#
mutable struct Norm1Group
  gens::Vector{Hecke.NfOrdFracIdl}
  rels
  A::RelNeq
  gC::Vector{Tuple{Hecke.NfOrdFracIdl, GrpAbFinGenElem}}
  sC::Tuple{GrpAbFinGen, Hecke.GrpAbFinGenMap}
  gU::Vector{Tuple{FacElem{nf_elem, AnticNumberField}, GrpAbFinGenElem}}
  sU::Tuple{GrpAbFinGen, Hecke.GrpAbFinGenMap}
  C::Any
  U::Any

  function Norm1Group(A::RelNeq)
    r = new()
    r.A = A
    r.C = class_group(A.K)
    c, mc = r.C
    r.sC = sub(c, elem_type(c)[])
    U, mU = Hecke.unit_group_fac_elem(maximal_order(A.K))
    u, mu = Hecke.unit_group_fac_elem(maximal_order(A.k))
    q, mq = quo(u, [mu\norm(A.m_k_K, mU(U[i])) for i=1:ngens(U)])
    r.U = q, pseudo_inv(mq)*mu
    r.sU = sub(u, elem_type(u)[])
    r.gC = [(ideal(maximal_order(A.K), 1)//1, 0*c[1])]
    r.gU = [(FacElem(A.k(1)), 0*q[1])]
    r.gens = []
    return r
  end
end

function Base.show(io::IO, N::Norm1Group)
  println(io, "Norm-1-class group")
  println(io, snf(N.sC[1])[1], " by ", snf(N.sU[1])[1])
  println(io, "currently, using $(length(N.gens)) generators")
end

function Hecke.isprincipal_fac_elem(A::FacElem{<:NfAbsOrdIdl})
  a,b = Hecke.reduce_ideal(A)
  # a*b == A
  fl, c = Hecke.isprincipal_fac_elem(a)
  if !fl
    return fl, c
  end
  return fl, c*b
end

function Hecke.isprincipal_fac_elem(A::FacElem{<:Hecke.NfOrdFracIdl})
  zk = order(base_ring(A))
  B = FacElem(Dict((numerator(x), v) for (x,v) = A.fac))
  den = Dict{nf_elem, fmpz}()
  for (x,v) = A.fac
    k = nf(zk)(denominator(x))
    if haskey(den, k)
      den[k] += v
    else
      den[k] = v
    end
  end
  #TODO: redude_ideal for FracIdl as well
  a,b = Hecke.reduce_ideal(B)
  # a*b == B = A*den
  fl, c = Hecke.isprincipal_fac_elem(a)
  if !fl
    return fl, c
  end
  return fl, c*b*inv(FacElem(den))
end

function Base.push!(N::Norm1Group, I::Hecke.NfOrdFracIdl)
  A = N.A
  @assert isone(norm(A.m_k_K, I))
  c, mc = N.C
  u, mu = N.U
  r = mc\numerator(I)
  fl, s = haspreimage(N.sC[2], r)
  if fl # found new relation
    J = FacElem(Dict((N.gC[i][1], s.coeff[1, i]) for i=1:ngens(N.sC[1])))
    J = I*inv(J)
    fl, g = Hecke.isprincipal_fac_elem(J)
    @assert fl
    ng = norm(A.m_k_K, g)
    @assert isunit(maximal_order(N.A.k)(evaluate(ng)))
    r = mu\ng
    fl, _ = haspreimage(N.sU[2], r)
    if fl
      return false # nothing new
    end
    push!(N.gens, I)
    push!(N.gU, (g, r))
    N.sU = sub(u, [x[2] for x = N.gU])
    return true
  else # found new generator
    push!(N.gens, I)
    push!(N.gC, (I, r))
    N.sC = sub(c, [x[2] for x = N.gC])
    return true
  end
end

function Hecke.order(N::Norm1Group)
  return order(N.sU[1]) * order(N.sC[1])
end

function order_bound(N::Norm1Group)
  return order(N.U[1]) * order(N.C[1])
end

Hecke.elem_type(::Type{Hecke.NfOrdFracIdlSet}) = Hecke.NfOrdFracIdl

function Hecke.evaluate(N::Norm1Group)
  # want the group extension (and the disc log and such)
  s1, ms1 = snf(N.sC[1])
  s2, ms2 = snf(N.sU[1])
  R = [rels(s2) zero_matrix(FlintZZ, ngens(s2), ngens(s1));
       zero_matrix(FlintZZ, ngens(s1), ngens(s2)) rels(s1)]

  for i = 1:ngens(s1)
    x = ms1(s1[i])
    I = FacElem(Dict((N.gC[j][1], x[j]) for j=1:ngens(N.sC[1])))
    I = I^order(x)
    fl, g = Hecke.isprincipal_fac_elem(I)
    @assert fl
    ng = norm(N.A.m_k_K, g)
    @assert isunit(maximal_order(N.A.k)(evaluate(ng)))
    r = N.U[2]\ng
    fl, x = haspreimage(N.sU[2], r)
    for j=1:ngens(s2)
      R[ngens(s2) + i, j] = -x[j]
    end
  end
  A = abelian_group(R)
  ZK = maximal_order(N.A.K)
  function exp(a::GrpAbFinGenElem)
    a1 = sub(a.coeff, 1:1, 1:ngens(s2))
    a2 = sub(a.coeff, 1:1, ngens(s2)+(1:ngens(s1)))
    b1 = ms2(s2(a1))
    b2 = ms1(s1(a2))
    I1 = FacElem(Dict((N.gC[i][1], b2[i]) for i=1:ngens(N.sC[1])))
    I2 = prod(ideal(ZK, N.gU[i][1]) ^ b1[i] for i=1:ngens(N.sU[1]))
    return I1*I2
  end

  function log(I::Hecke.NfOrdFracIdl)
    @assert isone(norm(N.A.m_k_K, I))
    r = N.C[2]\numerator(I)
    fl, s = haspreimage(N.sC[2], r)
    @assert fl
    J = FacElem(Dict((N.gC[i][1], s.coeff[1, i]) for i=1:ngens(N.sC[1])))
    J = I*inv(J)
    fl, g = Hecke.isprincipal_fac_elem(J)
    @assert fl
    ng = norm(N.A.m_k_K, g)
    @assert isunit(maximal_order(N.A.k)(evaluate(ng)))
    r = N.U[2]\ng
    fl, r = haspreimage(N.sU[2], r)
    @assert fl
    return A(hcat((ms2\r).coeff, (ms1\s).coeff))
  end
  return A, exp, log
end

function n1group(A::RelNeq, B::Int)
  K = A.K
  k = A.k
  mp = A.m_k_K

  S = PrimesSet(2, -1)
  max_stable = 20
  stable = max_stable
  zk = maximal_order(k)
  ZK = maximal_order(K)
  N = Norm1Group(A)

  function single_prime(P::NfAbsOrdIdl)
    p = minimum(P)
    if numerator(discriminant(K)) % p == 0 ||
       numerator(discriminant(k)) % p == 0
       @show "expensive", p
      lq = collect(factor(IdealSet(ZK)(A.m_k_K, P)))
    else
      lq = Hecke.prime_decomposition_nonindex(A.m_k_K, P)
    end
    f = matrix(FlintZZ, 1, length(lq), fmpz[div(degree(Q[1]), degree(P)) for Q = lq])
    r, n = nullspace(f)
    res = false
    for i = 1:r
      I = evaluate(FacElem(Dict((lq[j][1], n[j,i]) for j = 1:length(lq))), coprime = true)
      res = push!(N, I) || res
    end
    return res
  end

#TODO: missing: we NEED the ramified primes...
  for p = keys(factor(numerator(discriminant(ZK))).fac) #TODO: rel disc
    lp = prime_decomposition(zk, p)
    for P = lp
      if single_prime(P[1])
        stable = max_stable
      else
        stable -= 1
      end
    end
    if stable <= 0
      break
    end
  end

  d = lcm(numerator(discriminant(k)), numerator(discriminant(K)))

  for p = S
    if d % p == 0
      continue
    end
    lp = prime_decomposition(zk, p)
    for P = lp
      if single_prime(P[1])
        stable = max_stable
      else
        stable -= 1
      end
    end
    if stable <= 0
      break
    end
  end
  return N
end

function doit(f::fmpz_poly)
  K, a = number_field(f, cached = false)
  x = gen(parent(K.pol))
  k, b = number_field(x-1, cached = false)
  R = RelNeq(k, K)
  N = n1group(R, 10)
  C, _ = evaluate(N)
  return C
end

Zx, x = FlintZZ["x"]
function doit(n::String)
  f = open(n, "r")
  fo = open("$n.out", "w")
  i = 1
  while true
    @show l = readline(f)
    g = eval(parse(l))
    C = doit(g)
    C = snf(C)[1]
    println(fo, "$l -> $(C.snf)")
    println("$l -> $(C.snf)")
    i += 1
    if i % 10 == 0
      flush(fo)
    end
  end
end

end