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

# ╔═╡ 3f1cfd12-7b86-11eb-1371-c5795b87ef5b
begin
	using Pkg; Pkg.activate(".");
	function ingredients(path::String)
		# this is from the Julia source code (evalfile in base/loading.jl)
		# but with the modification that it returns the module instead of the last object
		name = Symbol(basename(path))
		m = Module(name)
		Core.eval(m,
	        Expr(:toplevel,
	             :(eval(x) = $(Expr(:core, :eval))($name, x)),
	             :(include(x) = $(Expr(:top, :include))($name, x)),
	             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
	             :(include($path))))
		m
	end
	using Plots, LaTeXStrings, LinearAlgebra, PlutoUI
	Tools = ingredients("tools.jl")
end;

# ╔═╡ 5a4c1582-7b86-11eb-2886-1547e93553dd
md"""
# §1 Approximation with Trigonometric Polynomials

## §1.1 Introduction

Trigonometric polynomials are functions of the form
```math
	t_N(x) = \sum_{k = -N}^N c_k e^{i k x}
```
Alternatively, we may write this as
```math
	t_N(x) = a_0 + \sum_{k = 1}^N \Big( a_k \cos(k x) + b_k \sin(k x) \Big),
```
but we will always use the much more convenient complex exponential notation.

Note that $t_N$ will be $2\pi$-periodic, i.e. we are now restricting ourselves to approximating periodic functions. But in return we gain a cleaner and simpler theory. Further, everything we do applies to some extent also to algebraic polynomials.

``N`` is called the degree of ``t_N``. The space of all trigonometric polynomials of degree ``N`` is denoted by
```math
\mathcal{T}_N := {\rm span}\big\{ x \mapsto \exp(i k x ) \,|\, k  = -N, \dots, N \big\}.
```
"""

# ╔═╡ 6030ee62-7c87-11eb-0699-67d75da87c97
md"""
Remember we said that approximation theory is about approximating complex functions with simple functions. Are trigonometric polynomials "simple" ? In other words, is ``e^{i k x}`` straightforward to implement? If ``z = e^{i x}`` is available then we can obtain ``e^{i k x} = z^k`` just by multiplication. Implementing ``e^{i x}`` is not difficult either, e.g., we can use a truncated Taylor series; see the exercise 1.1.1. But there are much more efficient implementations using rational functions. We will not pursue this but just accept that somebody has done that job for us.
"""

# ╔═╡ 80315186-7c87-11eb-0053-4bd9b0ada49f
md"""
### Function classes

The widest possible class of target functions ``f`` we will be concerned with are continuous functions, i.e., functions belonging to the space
```math
  C_{\rm per} := \big\{ f : \mathbb{R} \to \mathbb{R} \,|\, \text{ continuous, and } 2\pi\text{-periodic} \big\}.
```
Technically we could go to classes of integrable functions, ``L^p``, but this leads to further technicalities that less interesting.

We will assume higher regularity later on, such as differentiability of certain order, analyticity, etc. But for now it is only important for us to know that all ``f`` we encounter have well-defined point values.

The natural norm on ``C_{\rm per}`` is the max-norm (or, sup-norm, or ``L^\infty``-norm),
```math
	\| f \|_\infty := \max_{x \in [-\pi, \pi]} |f(x)|.
```
This is also a natural measure for the approximation error: throughout this course we will focus on understanding approximation errors in this norm, i.e., given a function ``f \in C_{\rm per}`` we aim to construct an approximation ``t_N \in \mathcal{T}_N`` such that
```math
   \|f - t_N \|_\infty
```
is small, and to understand the rate of this error as ``N \to \infty``.

Another very natural measure of the error is the ``L^2``-norm, ``\|f\|_2 := (\int_{-\pi}^\pi |f|^2 \,dx)^{1/2}``, which we will explore in §3 on least squares methods.
"""

# ╔═╡ 80158ed0-7b86-11eb-15c7-051e12dc4bd0
md"""
## § 1.2 Trigonometric Interpolation

A general degree ``N`` trigonometric polynomial,
```math
   t_N(x) = \sum_{k = -N}^N c_k e^{i k x}
```
has ``2N+1`` parameters ``c_{-N}, \dots, c_N``. How should we determine these? It would be ideal if we can prescribe exactly ``2N+1`` conditions. A natural and general approach is *interpolation*: given a target function ``f \in C_{\rm per}`` we demand that
```math
	t_N(x_j) = f(x_j), \qquad j = 0, \dots, 2N
```
where ``x_j`` are called the interpolation nodes.

How should they be chosen? It turns out  that equi-spaced nodes work very well. An intuitive justification for this choice is that in a periodic domain, all parts of the domain are "equal" and should be treated the same. By contrast in a finite interval one should *not* use equispaced nodes, but rather cluster them at the boundary (cf. Chebyshev interpolation [LN, Ch. 4]). Thus, we may choose the interpolation nodes
```math
	x_j = \frac{\pi j}{N}, \qquad j = 0, \dots, 2N.
```
(Note we could equally use equispaced notes on ``[-\pi, \pi]`` which often feels more natural, but the convention for certain algorithms -- see below -- is to use the nodes defined above.)
"""

# ╔═╡ 172b665d-2532-4847-ab2a-b4aaa36fc13e
md"""
Because of the periodic boundary condition the nodes ``x_0 = 0, x_{2N} = 2\pi`` are "equivalent", in the sense that ``f(x_0) = f(x_{2N})``, which is clearly a bad idea! A possibly way forward is to use instead the nodes ``\frac{2\pi j}{2N+1}, j= 1, \dots, 2N+1``. This might work (I haven't actually tested it!) but for algorithmic reasons it turns out that a much more convenient decision is to keep the nodes we defined above, but change the space of polynomials.

Namely, consider the two basis functions ``e^{iNx}`` and ``e^{-iNx}``. What values do they take on the grid?
```math
   e^{i N x_j} = e^{i N j \pi/N}  = e^{i\pi j}
  = (-1)^j = (-1)^{-j} = e^{-i \pi j} = e^{- i N \pi j/N} = e^{- i N x_j}.
```
That is, these two basis functions are identical on the grid. In other words we have found the kernel of the interpolation operator on the equispaced grid. There are multiple possible solution, but a simple one is to replace these basis functions with their mean, i.e.
```math
	\frac{e^{i N x} + e^{-i N x}}{2} = \cos(N x),
```
which of course takes again the same values on the grid. This results in a modified trigonometric polynomial space which contains the range of the interpolation operator,
```math
	\mathcal{T}_N' := {\rm span} \Big( \mathcal{T}_{N-1} \cup \{ \cos(Nx) \} \Big)
```
"""

# ╔═╡ 4fb28d6c-5de0-4b6d-9554-300bfafa4339
md"""
**Theorem:** The nodal interpolation operator ``I_N : C_{\rm per} \to \mathcal{T}_N'`` is well-defined, i.e., for all ``f \in C_{\rm per}`` there exists a unique ``I_N f := t_N \in \mathcal{T}_N`` such that
```math
	t_N(x_j) = f(x_j) \qquad \text{for } j = 1, \dots, 2N.
```

**Proof:** Straightforward exercise, e.g., a simple way to prove it is that it can be reduced to algebraic polynomial interpolation in ``\mathbb{C}``. It will also follow from results we will cover later.
"""

# ╔═╡ 7c7f18b4-854c-11eb-36f8-4f28cfdf4787
md"""
#### Implementation Notes

* To define the nodal interpolation operator we setup a linear system that specifies the coefficients. If ``I_N f(x) = \sum_{k = -N+1}^{N-1} \hat{F}_k e^{i k x} + \hat{F}_N \cos(Nx)`` then
```math
	\sum_{k = -N+1}^N \hat{F}_k e^{i k x_j} = F_j \qquad \text{for } j = 0, \dots, 2N-1,
```
where ``F_j = f(x_j)`` and ``\hat{F}_k`` are the coefficients. Knowing the coefficients we can then evaluate the interpolant.
* Since ``e^{i N x} = \cos(Nx)`` on the interpolation nodes, it doesn't matter which of the two we use to construct the interpolation operator.
* The ordering of the basis is in principle arbitrary. Here we use a convention used for fast algorithms (FFT),
```math
	(0, 1, \dots, N, -N+1, -N+1, \dots, -1)
```
This may look strange at first, but see more on this below!
"""

# ╔═╡ ad6c64a6-88e4-11eb-2c84-bf6926f78472
begin
	"interpolation nodes"
	xgrid(N) = range(0, 2*π - π/N, length = 2*N)

	"fourier coefficient indices"
	kgrid(N) = [0:N; -N+1:-1]
end;

# ╔═╡ 378d10e4-82f2-11eb-3d04-21f28b84d04a
let N = 3
	# implement the nodal interpolation operator
	A = [ exp(im * k * x) for k in kgrid(N), x in xgrid(N) ]
	# observe that A'A ~ diagonal
	real.(round.(A' * A, digits=12))
	# confirm it is an orthogonal matrix (up to scaling)! I.e. (A'/2N) = inv(A) !!
	# norm(A' * A - 2*N*I)
end

# ╔═╡ c7c4850e-82f2-11eb-16b9-a9e3b742249b
md"""
In addition to guaranteeing that ``I_N`` is well-defined we also see that the matrix ``A`` we need to invert is orthogonal (up to rescaling), which makes it very easy to invert it. We just need to multiply by ``A^H``, i.e. ``O(N^2)`` computational cost instead of ``O(N^3)`` for solving a full linear system via [Gaussian elimination](https://en.wikipedia.org/wiki/LU_decomposition).

These two operations ``F \mapsto \hat{F}`` and ``\hat{F} \mapsto F`` are called the discrete and inverse discrete fourier transforms. They can in fact be applied with ``O(N \log(N))`` computational cost, using the *fast fourier transform*. We will study this in the next lecture but for now use our naive implementation.
"""

# ╔═╡ 46d8f304-82f1-11eb-2675-a30b92a59ba2
begin
	"""
	construct the coefficients of the trigonometric interpolant
	"""
	function triginterp(f, N)
		X = xgrid(N)
		# nodal values at interpolation nodes
		F = f.(X)
		# system matrix
		A = [ exp(im * x * k) for k in kgrid(N), x in X ]
		# coefficients are given by F̂ = A' * F as discussed above!
		return (A' * F) / (2*N)
	end

	"""
	to evaluate a trigonometric polynomial just sum coefficients * basis
	we the take the real part because we assume the function we are
	approximating is real.
	"""
	evaltrig(x, F̂) = sum( real(F̂k * exp(im * x * k))
						  for (F̂k, k) in zip(F̂, kgrid(length(F̂) ÷ 2)) )

end

# ╔═╡ 39cbc3be-82f4-11eb-17a9-d15512292184
md"""
We are now ready to study some examples. Here is a selection of different functions we may want to approximate:
```math
\begin{aligned}
f_1(x) &= \sin(x) \\
f_2(x) &= \sin(4x) \\
f_3(x) &= |\sin(2x)| \\
f_4(x) &= |\sin(2x)|^3 \\
f_5(x) &= \exp(\cos(x)) \\
f_6(x) &= \frac{1}{1 + \sin(x)^2} \\
f_7(x) &= \frac{1}{1 + 10*\sin(x)^2} \\
f_8(x) &= \chi_{|x| < 2.5}  \exp\bigg(3 - \frac{3}{1-(x/2.5)^2}\bigg)
\end{aligned}
```
"""

# ╔═╡ 090507e4-82f5-11eb-0c51-f9a005e8d8d7
begin
	f1(x) = sin(x)
	f2(x) = sin(4*x)
	f3(x) = abs(sin(x))
	f4(x) = abs(sin(x))^3
	f5(x) = exp(cos(x))
	f6(x) = 1 / (1 + sin(x)^2)
	f7(x) = 1 / (1.0 + 10*sin(x)^2)
	g8(x) = x^2 < 2.49999^2 ? exp(3 - 3 / (1-(x/2.5)^2)) : 0.0
	f8(x) = g8(mod(x+π,2*π)-π)

	allf = [f1, f2, f3, f4, f5, f6, f7, f8]
	flabels = [L"f_{%$n}" for n in 1:8]

	xx = range(0, 2π, length=500)
	Ps = [ plot(xx, f.(xx), label = flabels[n], lw=2) for (n, f) in enumerate(allf) ]
	plot(Ps..., layout = (4,2))

end

# ╔═╡ e0309958-834b-11eb-041d-b7e8ad93a57b
md"""
We can just visually play around a bit with different target functions and degrees.....

Choose the function: $( @bind _fidx Slider(1:8) )

Choose the degree: $( @bind _N Slider(3:50) )
"""

# ╔═╡ e2603810-82f5-11eb-381b-7193adb06df2
let f = allf[_fidx], N = _N   # play with the functions and with the parameters!
	X = xgrid(N)
	xp = range(0, 2π, length=200)
	F̂ = triginterp(f, N)
	plot(xp, f.(xp), lw=6, label = L"f_{%$_fidx}", size = (450, 250))
	plot!(xp, evaltrig.(xp, Ref(F̂)), lw=3, label = L"I_{%$N} f")
	plot!(X, f.(X), lw=0, m=:o, ms=4, c = 2, label = "",
		 legend = :outertopright)
	title!(latexstring("\\mathrm{Trigonometric~Interpolant:} f_$_fidx, N = $N"))
end

# ╔═╡ f752a0d8-834b-11eb-1fb0-9b84e7b92c3f
md"""
....but this is not very informative. We need a more disciplined approach to measuring accuracy and evaluating the quality of the approximation scheme.

We will now study how the errors ``\| f - t_N \|_\infty`` behave asymptotically as ``N \to \infty``. This asymptotic behavious is called the *convergence rate*.
"""

# ╔═╡ 17eddcf2-835b-11eb-0335-35c09eb30e68
"""
implementation of a basic error function
"""
function triginterperror(f, N; Nerr = 1_362)
	xerr = range(0, 2π, length=Nerr)    # compute error on this grid
	F̂ = triginterp(f, N)                # trigonometric interpolant
	return norm(f.(xerr) - evaltrig.(xerr, Ref(F̂)), Inf)  # max-error on xerr grid
end;

# ╔═╡ 54ad3b38-834c-11eb-082c-cfa02a9a43ec
let NN = [1:2:5; (2).^(3:7)], flist = allf, flabels = flabels
	P = plot( xaxis  = (L"N~{\rm (degree)}", ),
    	      yaxis  = (:log, L"\Vert f - I_N f ~\Vert_{\infty}"),
        	  legend = :outertopright,
			  xscale = :log10, yscale = :log10,
			  size = (500, 300))
	for (f, lab) in zip(flist, flabels)
    	err = triginterperror.(f, NN)
		plot!(P, NN, err, lw=2, m=:o, ms=3, label = lab)
	end
	hline!(P, [1e-15], c = :red, lw=3, label = L"\epsilon")
end

# ╔═╡ 35b22f54-835e-11eb-11c7-215ac7ffcf5a
md"""
Straight lines on a log-log scale correspond to an *algebraic rate of convergence*
and the slope tells us what the rate is:
```math
	\log({\rm err}(N)) \sim p \log(N)
\qquad \Leftrightarrow  \qquad
	{\rm err}(N) \sim N^p
```
"""

# ╔═╡ 090729c0-835e-11eb-2957-617e3585eb86
let NN = (2).^(2:7), flist = allf[[3,4]], flabels = flabels[[3,4]]
	P = plot( xaxis  = (L"N~{\rm (degree)}", ),
    	      yaxis  = (:log, L"\Vert f - I_N f ~\Vert_{\infty}"),
        	  legend = :outertopright,
			  xscale = :log10, yscale = :log10,
			  size = (400, 250))
	for (f, lab) in zip(flist, flabels)
    	err = triginterperror.(f, NN)
		plot!(P, NN, err, lw=2, m=:o, ms=3, label = lab)
	end
	P
end

# ╔═╡ 8df5ebe4-835e-11eb-0fda-1314a49462e3
md"""
Sublinear behaviour on a log-log plot means super-algebraic convergence. There are many examples of that kind, but the simplest (and typical) is exponentioal convergence,
```math
	{\rm err}(N) \sim e^{- \alpha N}
	\qquad \Leftrightarrow \qquad
	\log({\rm err}(N)) \sim - \alpha N
```
i.e. the "rate" ``\alpha`` can be seen on a lin-log plot.
"""

# ╔═╡ e6df4f7a-835e-11eb-2fbd-15acdc141f47
let NN = 4:4:40, flist = allf[5:8], flabels = flabels[5:8]
	# check f8 by switching to NN = 10:10:100
	P = plot( xaxis  = (L"N~{\rm (degree)}", ),
    	      yaxis  = (:log, L"\Vert f - I_N f ~\Vert_{\infty}"),
        	  legend = :outertopright,
			  yscale = :log10,
			  size = (400, 250))
	for (f, lab) in zip(flist, flabels)
    	err = triginterperror.(f, NN)
		plot!(P, NN, err, lw=2, m=:o, ms=3, label = lab)
	end
	P
end

# ╔═╡ 67d305a4-7d30-11eb-2bac-4d26774d640c
md"""
## § 1.4 Approximation Error Estimates

We have seen how to construct trigonometric polynomial approximants and have observed several different kinds of convergence behaviour. Our last goal for the first set of lectures is to explain these observations rigorously. This material is contained in my lecture notes (with further references to the literature). I won't reproduce it here in much detail, but only summarize the main results.

Here we continue to focus on max-norm approximation, but will cover 2-norm approximations in the exercises.
"""

# ╔═╡ 63d298d8-836f-11eb-0a86-9bc0a3665311
md"""
### § 1.4.1 Best approximation and near-best approximation

We will momentarily formulate several results in terms of *best approximation*, i.e., we will show results along the lines
```math
   \inf_{t_N \in \mathcal{T}_N} \|f - t_N \|_\infty \leq \epsilon_N,
```
where ``\epsilon_N`` is some rate of convergence. Or equivalently, *there exist ``t_N \in \mathcal{T}_N`` such that* ``\| f - t_N\|_\infty \leq \epsilon_N``. But it will not be clear from those results how to construct such best or quasi-best approximants.

Here, a general principle comes into play that saves us: Let ``I_N f`` denote the degree-``N`` trigonometric interpolant of ``f``, then one can prove that
```math
	\| I_N f \|_\infty \leq \Lambda_N \| f \|_\infty \qquad \forall f \in C_{\rm per}
```
The constant ``\Lambda_N`` is called the Lebesgue constant and can be precisely estimated:

**Theorem:** ``\Lambda_N \leq \frac{2}{\pi} \log(N+1) + 1``.

**Proof:** See [LN, Thm 3.26] for a proof of a slightly weaker result. The proof is not difficult but quite dry, involving delicate estimates of the Dirichlet kernel. Please study this for further reading.

In practise we will just write ``\Lambda_N \lesssim \log N``. With this in hand we have the following argument at our disposal: Let ``t_N \in \mathcal{T}_N'`` be *any* trigonometric polynomial, then ``I_N t_N = t_N``, and hence,
```math
\begin{aligned}
	\|f - I_N f \|_\infty
	&\leq
	\|f - t_N\|_\infty + \| t_N - I_N f \|_\infty \\
    &= \|f - t_N\|_\infty + \| I_N (t_N - f) \|_\infty \\
	&\leq (1 + \Lambda_N) \|f - t_N \|_\infty.
\end{aligned}
```
Now, taking the infimum over all ``t_N \in \mathcal{T}_N'`` we obtain that
```math
   \|f - I_N f \|_\infty \leq (1+\Lambda_N) \inf_{t_N \in \mathcal{T}_N'} \|f - t_N \|_\infty.
```
Thus, the interpolation error is within a ``\log N`` factor of the best approximation error. In practise this cannot be observed, and indeed in many scenarios this factor can in fact be removed entirely. Because it is so mild, we won't worry about it.

**Remark:** The foregoing argument is extremely general. The only property that we used here is that the interpolation operator ``I_N : C_{\rm per} \to \mathcal{T}_N'`` is a bounded operator with respect to the ``\|\cdot\|_\infty`` norm. It is therefore easy to generalize and indeed used extensively throughout numerical analysis.

"""

# ╔═╡ f265a3ca-836e-11eb-2a3a-b5eba6541b40
md"""
### § 1.4.2 Jackson Theorems [LN, Sec. 3.3]

A function ``f \in C_{\rm per}`` has modulus of continuity ``\omega : \mathbb{R} \to \mathbb{R}`` if
```math
|f(x) - f(x')| \leq \omega(|x - x'|) \qquad \forall x, x' \in \mathbb{R}.
```

**First Jackson Theorem:** If ``f \in C_{\rm per}`` has modulus of continuity ``\omega``, then there exists ``t_N \in \mathcal{T}_N`` such that
```math
	\| f - t_N \|_\infty \omega(N^{-1}).
```

The second Jackson theorem extends the result to functions that are continuously differentiable. We say that ``f \in C^p_{\rm per}`` if ``f \in C_{\rm per}``, is ``p`` times continuously differentiable on ``\mathbb{R}`` and all derivatives up to ``f^{(p)} \in C_{\rm per}``.

**Second Jackson Theorem:** Assume that ``f \in C_{\rm per}^p`` and ``f^{(p)}`` has modulus of continuity ``\omega``, then there exists ``t_N \in \mathcal{T}_N`` such that
```math
	\| f - t_N \|_\infty N^{-p} \omega(N^{-1}).
```

A few comments and special cases:
* The smoother ``f`` the faster the convergence.
* These results are qualitatively sharp, i.e. Jackson's construction has the optimal rate for these function classes.
* If ``f`` is Lipschitz then ``\omega(r) = L r`` where ``L`` is the Lipschitz constant. An example of a Lipschitz continuous function is ``f(x) = |\sin(x)|``.
* If ``f`` is ``\sigma``-Hölder then ``\omega(r) = C r^\sigma``. An example of a Hölder continuous function is ``f(x) = |\sin(x)|^\sigma``.
"""

# ╔═╡ a5f2dc0a-859a-11eb-0fa7-6b32c67bf953
md"""
We can now explain the first set of numerical results:
* ``f_3`` is Lipschitz continuous, hence the rate is ``N^{-1}``
* ``f_4 \in C^2`` with ``f_4''`` Lischitz hence the rate is ``N^{-3}``
"""

# ╔═╡ ee4409ac-859a-11eb-1692-af3bf2dcd711
let NN = (2).^(2:7), flist = allf[[3,4]], flabels = flabels[[3,4]]
	P = plot( xaxis  = (L"N~{\rm (degree)}", ),
    	      yaxis  = (:log, L"\Vert f - I_N f ~\Vert_{\infty}"),
        	  legend = :outertopright,
			  xscale = :log10, yscale = :log10,
			  size = (400, 250))
	for (f, lab) in zip(flist, flabels)
    	err = triginterperror.(f, NN)
		plot!(P, NN, err, lw=2, m=:o, ms=3, label = lab)
	end
	nn = NN[3:end]
	plot!(nn, 1.3*nn.^(-1), lw=2, c=:black, ls=:dash, label = L"N^{-1}, N^{-3}")
	plot!(nn, 3*nn.^(-3), lw=2, c = :black, ls=:dash, label = "")
	P
end

# ╔═╡ 7b57d3ee-8525-11eb-29fe-3912c6b3a866
md"""
### § 1.4.3 Paley Wiener Theorems [LN, Sec. 2.3 and 3.4]

First we review some standard results of complex analysis that we cannot reproduce here in detail. They will be covered in any introductory textbook or lecture notes on complex analysis.

A function ``f : \mathbb{R} \to \mathbb{R}`` is called analytic in a point ``x_0`` if it is represented by a Taylor series, i.e.,
```math
	f(x) = \sum_{n = 0}^\infty c_n (x - x_0)^n
```
and the series has radius of convergence ``r > 0``. If ``f \in C_{\rm per}`` is analytic in ``[-\pi, \pi]`` (and hence in all of ``\mathbb{R}``) then using the Taylor series one can extend it to a function on a strip in the complex plane: Let
```math
\Omega_\alpha := \{ z \in \mathbb{C} : |{\rm Im} z| < \alpha \}
```
then there exists ``\alpha > 0`` and unique analytic ``f : \Omega_\alpha \to \mathbb{C}`` which agrees with the original ``f`` in ``\mathbb{R}``.

**Payley Wiener Theorem:** Let ``f \in C_{\rm per}`` have an analytic extension to ``\bar{\Omega}_\alpha``, then
```math
	f(x) = \sum_{k \in \mathbb{Z}} \hat{f}_k e^{i k x}
	\qquad \text{with} \qquad |\hat{f}_k| \leq M_\alpha e^{- \alpha |k|}
```
where ``M_\alpha = \sup_{z \in \bar{\Omega}_\alpha} |f(z)|``.

**Corollary:**  Let ``f \in C_{\rm per}`` have an analytic extension to ``\bar{\Omega}_\alpha``, then there exists ``t_N \in \mathcal{T}_N`` such that
```math
	\|f - t_N\|_\infty \leq M_\alpha e^{-\alpha N}
```

The most important question for us is "what is the maximal extension"? I.e., what is the largest value ``\alpha`` such that ``f`` can be extended to an analytic function in ``\Omega_\alpha``? We will study this on the two examples ``f_6, f_7``.
"""

# ╔═╡ 37a559b6-859b-11eb-1723-19d70fe7d43c
md"""
We are now ready to explain the remaining numerical results:

* ``f_5`` is entire, i.e., the error is bounded by ``C_\alpha e^{-\alpha N}`` for all ``\alpha > 0``. This is called *superexponential* convergence.

* ``f_8 \in C^\infty`` but is not analytic. The point where it becomes zero has a Taylor expansion that evaluates to zero and hence has a zero convergence radius.

* ``f_6, f_7`` are both analytic. To determine the region of analyticity of a function ``(1+c^2 \sin^2(x))^{-1}`` we simply need to find where it has a pole.
"""

# ╔═╡ d90a3c54-85b4-11eb-0f40-5906527dd6e6
md"""
Use this slider to adjust the parameter ``c \in [1, 20]`` in the target function ``1 / (1+ c \sin^2(z))``.

$(@bind _c Slider(1:20))
"""

# ╔═╡ 09fc0678-03d9-45bd-944f-4594571f6bfa

# ╔═╡ 9c4c1bea-85a8-11eb-24bd-11e4c39fb5dc
let f = z -> 1 / (1 + _c * sin(z)^2)
	xp = range(-π, π, length=100)
	yp = range(-2, 2, length=100)
	contourf(xp, yp, (x,y) -> log(abs(f(x + im * y))),
		    xlabel = L"x = Re(z)", ylabel = L"y = Im(z)", size = (500, 400),
	 		title = latexstring("|(1+ $_c \\sin^2(z))^{-1}|"))
	hline!([0.0, asinh(1/sqrt(_c)), - asinh(1/sqrt(_c))], lw=2, c=:red, label = L"\mathbb{R}, \pm \alpha")
end

# ╔═╡ 991a71e2-85a8-11eb-1a78-0f4ec0e83a8e
md"""
To carry this out rigorously we need to solve ``c^2 \sin^2(z) = -1``, where ``c = 1`` for ``f_6`` and ``c = \sqrt{10}`` for ``f_7``. Thus we need to solve
```math
	\sin(z) = \pm \frac{i}{c}
```
This is not difficult but needs a bit of complex analysis. With a bit of work one can prove that the poles occur at
```math
	z = n \pi \pm i  \sinh^{-1}(1/c) \sim n\pi \pm i \frac{1}{c} \qquad \text{as } c \to \infty.
```
This gives us the *open* region of analyticity ``\Omega_\alpha`` with ``\alpha = \sinh^{-1}(1/c)``. To apply the Paley Wiener theorem we need to make ``\alpha`` just a little bit smaller, and we get
```math
	\| f - t_N \|_\infty \leq C_\epsilon e^{ - (\alpha - \epsilon) N } \qquad \forall \epsilon > 0.
```
But in practise we just ignore this ``\epsilon`` correction:
"""

# ╔═╡ 641fb412-85aa-11eb-095e-e1bf03a017f1
let NN = 6:6:60, flist = allf[6:7], flabels = flabels[6:7]
	P = plot( xaxis  = (L"N~{\rm (degree)}", ),
    	      yaxis  = (:log, L"\Vert f - I_N f ~\Vert_{\infty}"),
        	  legend = :outertopright,
			  yscale = :log10,
			  size = (450, 250))
	for (f, lab) in zip(flist, flabels)
    	err = triginterperror.(f, NN)
		plot!(P, NN, err, lw=2, m=:o, ms=3, label = lab)
	end
	nn1 = [20, 35]; nn2 = [30, 50]
	plot!(nn1, 10*exp.( - asinh(1.0) * nn1), lw=2, ls=:dash, c=:black, label = L"\exp(-\alpha_i N)")
	plot!(nn2, 10*exp.( - asinh(1/sqrt(10.0)) * nn2), lw=2, ls=:dash, c=:black, label = "")
end

# ╔═╡ 55d9bffc-85ab-11eb-0c50-83abd402e48a
md"""
This completes the first set of lectures on:
* trigonometric polynomials
* interpolation with trigonometric polynomials
* rates of convergence
* implementation and numerical investigation of convergence rates
* summary of rigorous results (Jackson, Paley-Wiener)
* matching theory and numerical results

The next topics will be
* fast evaluation of the trigonometric interpolant via the fast fourier transform
* numerical solution of BVPs and and IBVPs using trigonometric polynomials
"""

# ╔═╡ Cell order:
# ╟─3f1cfd12-7b86-11eb-1371-c5795b87ef5b
# ╟─5a4c1582-7b86-11eb-2886-1547e93553dd
# ╟─6030ee62-7c87-11eb-0699-67d75da87c97
# ╟─80315186-7c87-11eb-0053-4bd9b0ada49f
# ╟─80158ed0-7b86-11eb-15c7-051e12dc4bd0
# ╟─172b665d-2532-4847-ab2a-b4aaa36fc13e
# ╟─4fb28d6c-5de0-4b6d-9554-300bfafa4339
# ╟─7c7f18b4-854c-11eb-36f8-4f28cfdf4787
# ╠═ad6c64a6-88e4-11eb-2c84-bf6926f78472
# ╠═378d10e4-82f2-11eb-3d04-21f28b84d04a
# ╟─c7c4850e-82f2-11eb-16b9-a9e3b742249b
# ╠═46d8f304-82f1-11eb-2675-a30b92a59ba2
# ╟─39cbc3be-82f4-11eb-17a9-d15512292184
# ╠═090507e4-82f5-11eb-0c51-f9a005e8d8d7
# ╟─e0309958-834b-11eb-041d-b7e8ad93a57b
# ╟─e2603810-82f5-11eb-381b-7193adb06df2
# ╟─f752a0d8-834b-11eb-1fb0-9b84e7b92c3f
# ╠═17eddcf2-835b-11eb-0335-35c09eb30e68
# ╠═54ad3b38-834c-11eb-082c-cfa02a9a43ec
# ╟─35b22f54-835e-11eb-11c7-215ac7ffcf5a
# ╟─090729c0-835e-11eb-2957-617e3585eb86
# ╟─8df5ebe4-835e-11eb-0fda-1314a49462e3
# ╟─e6df4f7a-835e-11eb-2fbd-15acdc141f47
# ╟─67d305a4-7d30-11eb-2bac-4d26774d640c
# ╟─63d298d8-836f-11eb-0a86-9bc0a3665311
# ╟─f265a3ca-836e-11eb-2a3a-b5eba6541b40
# ╟─a5f2dc0a-859a-11eb-0fa7-6b32c67bf953
# ╟─ee4409ac-859a-11eb-1692-af3bf2dcd711
# ╟─7b57d3ee-8525-11eb-29fe-3912c6b3a866
# ╟─37a559b6-859b-11eb-1723-19d70fe7d43c
# ╟─d90a3c54-85b4-11eb-0f40-5906527dd6e6
# ╠═09fc0678-03d9-45bd-944f-4594571f6bfa
# ╟─9c4c1bea-85a8-11eb-24bd-11e4c39fb5dc
# ╟─991a71e2-85a8-11eb-1a78-0f4ec0e83a8e
# ╟─641fb412-85aa-11eb-095e-e1bf03a017f1
# ╟─55d9bffc-85ab-11eb-0c50-83abd402e48a

max_stars_repo_name: cortner/ApxThyApp

max_issues_repo_name: cortner/ApxThyApp

max_forks_repo_name: cortner/ApxThyApp


text:
@testset "Linear Regression" begin

    X₃₃ = [rand(0:2, 3) for _ in 1:10]
    X₃₃_test = [rand(0:2, 3) for _ in 1:100]

    B₉ = [rand(Bool, 9) for _ in 1:10]
    B₉_test = [rand(Bool, 9) for _ in 1:2000]

    training_configs = [
        Dict(
            :info => "No binarization on ⟦0,2⟧³",
            :train => X₃₃,
            :test => X₃₃_test,
            :encoding => :none,
            :binarization => :none,
        ),
        Dict(
            :info => "Domain Wall binarization on ⟦0,2⟧³",
            :train => X₃₃,
            :test => X₃₃_test,
            :encoding => :none,
            :binarization => :domain_wall,
        ),
        Dict(
            :info => "One-Hot binarization on ⟦0,2⟧³",
            :train => X₃₃,
            :test => X₃₃_test,
            :encoding => :none,
            :binarization => :one_hot,
        ),
        Dict(
            :info => "Domain Wall pre-encoded on ⟦0,2⟧³",
            :train => B₉,
            :test => B₉_test,
            :encoding => :domain_wall,
            :binarization => :none,
        ),
        Dict(
            :info => "One-Hot pre-encoded on ⟦0,2⟧³",
            :train => B₉,
            :test => B₉_test,
            :encoding => :one_hot,
            :binarization => :none,
        ),
    ]

    function all_different(x, encoding)
        encoding == :none && (return allunique(x))
        isv = if encoding == :one_hot
            mapreduce(i -> is_valid(x[i:i+2], Val(encoding)), *, 1:3:9)
        else
            mapreduce(i -> is_valid(x[i:i+1], Val(encoding)), *, 1:2:6)
        end
        if isv
            b = all_different(debinarize(x; binarization = encoding), :none)
            return b ? -1. : 1.
        else
            return length(x)
        end
    end

    function all_different(x, encoding, binarization)
        return all_different(x, encoding == :none ? binarization : encoding)
    end

    for config in training_configs
        @testset "$(config[:info])" begin
            println("\nTest for $(config[:info])")
            penalty = x -> all_different(x, config[:encoding], config[:binarization])
            train(
                config[:train], penalty;
                binarization = config[:binarization], X_test = config[:test]
            )
        end
    end

    # for config in training_configs
    #     # @testset "$(config[:info])" begin
    #         println("\nTest for $(config[:info])")
    #         penalty = x -> all_different(x, config[:encoding], config[:binarization])
    #         train(config[:train], penalty; binarization = config[:binarization])
    #     # end
    # end

    # @testset "All Different 3×3" begin
    #     println("\nTest for All Different 3×3")
    #     all_different(x) = allunique(x)

    #     binary_all_different(x) = all_different(integerize(x))

    #     opt_all_different(x) = 0.

    #     function sat_all_different(X)
    #         if count(x -> x == sqrt(length(X)), integerize(X)) > 0
    #             return sqrt(length(X))
    #         elseif binary_all_different(X)
    #             return -1.
    #         else
    #             return 1.
    #         end
    #     end

    #     for ss in set_sizes, ds in domain_sizes, st in set_types
    #         str = "ALL_DIFFERENT $ds×$ds: set size = $ss, set type = $(st[1])"
    #         println("\n\t$str\n")
    #         X_train = st[2](ds, ss)
    #         X_check = first(set_types)[2](ds, last(set_sizes))
    #         train(X_train, sat_all_different, opt_all_different; X_check)
    #     end
    # end

    # @testset "Ordered" begin
    #     println("\nTest for Ordered 3×3")
    #     ordered(x) = issorted(x)

    #     binary_ordered(x) = ordered(integerize(x))

    #     opt_ordered(x) = 0.

    #     function sat_ordered(X)
    #         if count(x -> x == sqrt(length(X)), integerize(X)) > 0
    #             return sqrt(length(X))
    #         elseif binary_ordered(X)
    #             return -1.
    #         else
    #             return 1.
    #         end
    #     end

    #     for ss in set_sizes, ds in domain_sizes, st in set_types
    #         str = "ORDERED $ds×$ds: set size = $ss, set type = $(st[1])"
    #         println("\n\t$str\n")
    #         X_train = st[2](ds, ss)
    #         X_check = first(set_types)[2](ds, last(set_sizes))
    #         train(X_train, sat_ordered, opt_ordered; X_check)
    #     end
    # end

    # @testset "LinearSum" begin
    #     function param_set(ds)
    #         N = ds^2
    #         return [
    #             0,
    #             ds,
    #             round(Int,1.5*ds),
    #             round(Int, .25*N),
    #             round(Int, .5*N),
    #             round(Int, .75*N),
    #             N,
    #         ]
    #     end

    #     for ss in set_sizes, ds in domain_sizes, st in set_types, σ in param_set(ds)
    #         linear_sum(x) = sum(x) == σ

    #         binary_linear_sum(x) = linear_sum(integerize(x))

    #         opt_linear_sum(x) = 0.

    #         function sat_linear_sum(X)
    #             if count(x -> x == sqrt(length(X)), integerize(X)) > 0
    #                 return sqrt(length(X))
    #             elseif binary_linear_sum(X)
    #                 return -1.
    #             else
    #                 return abs(sum(integerize(X)) - σ)
    #             end
    #         end

    #         QP = QUBO_linear_sum(ds, σ)

    #         qubo_handmade(X) = transpose(X) * QP * X + σ^2

    #         str = "LINEAR_SUM(=$σ) $ds×$ds: set size = $ss, set type = $(st[1])"
    #         println("\n\t$str\n")
    #         X_train = st[2](ds, ss)
    #         X_check = first(set_types)[2](ds, last(set_sizes))
    #         train(X_train, qubo_handmade, opt_linear_sum; X_check)
    #     end
    # end
end

ext: jl


text:
export AbstractApproximator, AbstractVApproximator, AbstractQApproximator

"""
An approximator is a functional object to estimate a state.
Two typical kinds of approximators are
[`AbstractVApproximator`](@ref) and [`AbstractQApproximator`](@ref).
"""
abstract type AbstractApproximator end

"""
A collection of approximators to estimate the value of a state.
"""
abstract type AbstractVApproximator <: AbstractApproximator end

"""
A collection of approximators to estimate the values of actions given a state.
"""
abstract type AbstractQApproximator <: AbstractApproximator end


text:
using ReversedStrings
using Test

@testset "ReversedStrings.jl" begin
    @testset "prevind, nextind" begin
        s = "abc"
        rs = ReversedString(s)
        for n in 1:5
            i = 0
            @test_throws BoundsError prevind(s,i,n)
            @test_throws BoundsError prevind(rs,i,n)
            for i in 1:4
                @test prevind(s,i,n) == prevind(rs,i,n)
            end
            for i in 0:3
                @test nextind(s,i,n) == nextind(rs,i,n)
            end
            i = 5
            @test_throws BoundsError nextind(s,i,n)
            @test_throws BoundsError nextind(rs,i,n)
        end

        s="aän"
        rs = ReversedString(s)
        rs |> collect

        @test s[nextind(s,1)] == rs[nextind(rs,1)]
    end
    @test [ c for c in reverse("JuliaCon") ] == [ c for c in reversed("JuliaCon") ]
    @test reversed(reversed("JuliaCon")) == "JuliaCon"
    @test ncodeunits(reverse("JuliaCon")) == ncodeunits(reversed("JuliaCon"))
    @test isvalid(reverse("JuliaCon"),1) == isvalid(reversed("JuliaCon"),1)
    @test isvalid(reverse("JuliaCon"),10) == isvalid(reversed("JuliaCon"),10)

    @test reverse("JuliaCon")[2:4] == reversed("JuliaCon")[2:4]

    @test SubString(reverse("JuliaCon"),2:4) == SubString(reversed("JuliaCon"),2:4)
    @test SubString(reverse("JuliaCon"),2,4) == SubString(reversed("JuliaCon"),2,4)
end


text:
using FountainCodes, Test, SparseArrays

function test_r10_degree()

    # test the values at which the degree increments
    f = [0, 10241, 491582, 712794, 831695, 948446, 1032189]
    correct = [1, 2, 3, 4, 10, 11, 40]
    for (v, d) in zip(f, correct)
        ans = FountainCodes.r10_degree(v)
        if ans != d
            error("expected r10_degree($v) to be $d, but got $ans.")
        end
    end

    # test the values one before the degree increments
    f = [10241, 491582, 712794, 831695, 948446, 1032189, 1048576]
    correct = [1, 2, 3, 4, 10, 11, 40]
    for (v, d) in zip(f, correct)
        ans = FountainCodes.r10_degree(v-1)
        if ans != d
            error("expected r10_degree($(v-1)) to be $d, but got $ans.")
        end
    end
    return true
end
@test test_r10_degree()

function test_ldpc_constraints()
    code = R10(10)
    A = FountainCodes.ldpc_constraint_matrix(code)

    # 1st LDPC constraint
    correct = sparsevec([1, 6, 7, 8, 11], ones(GF256, 5), dimension(code))
    ans = A[:, 1]
    @test ans == correct

    # 7th LDPC constraint
    correct = sparsevec([5, 6, 7, 10, 17], ones(GF256, 5), dimension(code))
    ans = A[:, 7]
    @test ans == correct
end
test_ldpc_constraints()

function test_hdpc_constraints()
    code = R10(10)
    A = FountainCodes.hdpc_constraint_matrix(code)

    # 1st HDPC constraint
    Is = [1, 2, 4, 5, 8, 10, 11, 15, 18]
    Vs = ones(GF256, length(Is))
    correct = sparsevec(Is, Vs, dimension(code))
    ans = A[:, 1]
    @test ans == correct

    # 4th HDPC constraint
    Is = [2, 3, 4, 5, 6, 7, 14, 15, 16, 17, 21]
    Vs = ones(Bool, length(Is))
    correct = sparsevec(Is, Vs, dimension(code))
    ans = A[:, 4]
    @test ans == correct

    # 6th HDPC constraint
    Is = [11, 12, 13, 14, 15, 16, 17, 23]
    Vs = ones(Bool, length(Is))
    correct = sparsevec(Is, Vs, dimension(code))
    ans = A[:, 6]
    @test ans == correct
    return
end
test_hdpc_constraints()

"""Test that decoding with ESIs 0 through K-1 succeeds"""
function test_precode(K::Integer)
    rng = MersenneTwister(123)
    code = R10(K)
    src = rand_nonzero(rng, GF256, K)

    # compute the intermediate symbols
    A = constraint_matrix(code, 0:K-1)
    pre = zeros(GF256, size(A, 2))
    pre[end-K+1:end] .= src
    int = decode(A, pre)

    # test that the first K coded symbols are equal to the src symbols, i.e., that the code is
    # systematic
    G = generator_matrix(code, 0:K-1)
    enc = G'*int
    @test enc[1:K] == src
end
for K in 4:100
    test_precode(K)
end
test_precode(7000)
test_precode(8192)

max_stars_repo_name: severinson/RaptorCodes

max_issues_repo_name: severinson/RaptorCodes

max_forks_repo_name: severinson/RaptorCodes


text:
function decode(str::String)::Array{String}
    tokenformat = r"\\\d+|."
    result = String[]
    idx = 1
    while idx ≤ length(str)
        token = match( tokenformat, str, idx )
        idx += length(token.match)
        if length(token.match) == 1
            push!( result, token.match )
        else
            n = parse(Int,token.match[2:end])
            push!( result, str[idx:min(idx+n-1,end)] )
            idx += n
        end
    end
    result
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using Pkg; Pkg.activate(pwd())
using Revise, BenchmarkTools, Printf
using Random, Distances

const PASS = 1
const BET = 2
const NUM_ACTIONS = 2

mutable struct Knode
    infoset::String
    regretsum::Vector{Float64}
    strategy::Vector{Float64}
    strategysum::Vector{Float64}
end

Knode(i::String) = Knode(i, zeros(NUM_ACTIONS), fill(1 / NUM_ACTIONS, NUM_ACTIONS), zeros(NUM_ACTIONS))

function getstrategy!(n::Knode, weight::Float64)
    normalsum = 0.0
    for a = 1:NUM_ACTIONS
        n.strategy[a] = n.regretsum[a] > 0 ? n.regretsum[a] : 0.0
        normalsum += n.strategy[a]
    end
    for a = 1:NUM_ACTIONS
        if normalsum > 0
            n.strategy[a] /= normalsum
        else
            n.strategy[a] = 1.0 / NUM_ACTIONS
        end
        n.strategysum[a] += weight * n.strategy[a]
    end
    n.strategy
end

function getstrategy_average(n::Knode)
    avgstrategy = Vector{Float64}(undef, NUM_ACTIONS)
    normalsum = 0.0
    for a = 1:NUM_ACTIONS
        normalsum += n.strategysum[a]
    end
    for a = 1:NUM_ACTIONS
        if normalsum > 0
            avgstrategy[a] = n.strategysum[a] / normalsum
        else
            avgstrategy[a] = 1.0 / NUM_ACTIONS
        end
    end
    avgstrategy
end

function infostring(n::Knode)
    @sprintf("%4s: %s", n.infoset, string(getstrategy_average(n)))
end

function nextnode!(nodedict::Dict{String, Knode}, infoset::String)
    if haskey(nodedict, infoset)
        n = nodedict[infoset]
    else
        n = Knode(infoset)
        nodedict[infoset] = n
    end
    n
end

function cfr!(cards::Vector{Int64}, h::String, p1::Float64, p2::Float64, nodedict::Dict{String,Knode})
    plays = length(h)
    player = rem(plays, 2) + 1
    opponent = 2 - player + 1
    if plays > 1
        terminal = (h[plays] == 'p')
        doublebet = h[(plays - 1):plays] == "bb"
        higher = cards[player] > cards[opponent]
        if terminal
            if h == "pp"
                payoff = higher ? 1 : -1
                # println("h: $h ", cards[1:2], " Player $player gets $payoff")
                return payoff
            else
                payoff = 1
                # println("h: $h ", cards[1:2], " Player $player gets $payoff")
                return payoff
            end
        elseif doublebet
            payoff = higher ? 2 : -2
            # println("h: $h ", cards[1:2], " Player $player gets $payoff")
            return payoff
        end
    end
    infoset = string(cards[player], h)
    n = nextnode!(nodedict, infoset)
    strategy = getstrategy!(n, player == 1 ? p1 : p2)
    # println("regretsum a: ", n.infoset, n.regretsum)
    util = Vector{Float64}(undef, NUM_ACTIONS)
    nodeutil = 0.0
    for a = 1:NUM_ACTIONS # recur cfr
        hnext = string(h, a == 1 ? "p" : "b")
        if player == 1
            util[a] = -1 * cfr!(cards, hnext, p1 * strategy[a], p2, nodedict)
        else
            util[a] = -1 * cfr!(cards, hnext, p1, p2 * strategy[a], nodedict)
        end
        nodeutil += strategy[a] * util[a]
    end
    for a = 1:NUM_ACTIONS # regret
        temp = util[a]
        # println(n.infoset, " util = $temp, nodeutil = $nodeutil")
        regret = util[a] - nodeutil
        n.regretsum[a] += (player == 1 ? p1 : p2) * regret
        # @show strategy[a]
    end
    # println("regretsum b: ", n.infoset, n.regretsum)
    return nodeutil
end

function train(iter::Int64)
    println("STARTING...")
    cards = [1, 2, 3]
    util = 0.0
    nodedict = Dict{String, Knode}()
    vals = Vector{Float64}(undef, iter)
    for i = 1:iter
        print("\rITER $i")
        shuffle!(cards)
        vals[i] = cfr!(cards, "", 1.0, 1.0, nodedict)
        util += vals[i]
    end
    println("Average game value: ", util / iter)
    for n in values(nodedict)
        println(infostring(n))
    end
end

function kuhncfr(reps::Int64)
    Random.seed!(23894719)
    train(reps)
end

@btime kuhncfr(10_000)

max_stars_repo_name: ajkeith/Cyber-Air-Defense

max_issues_repo_name: ajkeith/Cyber-Air-Defense


text:
function capitalize(cur_string::AbstractString)
  cur_string = join(map(x -> uppercasefirst(x), split(cur_string, "/")), "/")

  cur_string
end


text:
using MethodAnalysis
using OrderedCollections
using CodeTracking

const fully_specialized = OrderedSet{Core.MethodInstance}()
const nonspecialized = OrderedSet{Core.MethodInstance}()

function categorizemi(mi::Core.MethodInstance)
    bin = all(isconcretetype, Base.unwrap_unionall(mi.specTypes).parameters) ?
        fully_specialized : nonspecialized
    push!(bin, mi)
    return nothing
end

categorizemi(x) = nothing

visit(categorizemi)
println("Found ", length(fully_specialized), " MethodInstances and ", length(nonspecialized), " other MethodInstances")

# We're especially interested in the nonspecialized ones.
# Find specializations that have a TypeVar parameter
function hastv(typ)
    isa(typ, UnionAll) && return true
    if isa(typ, DataType)
        for p in typ.parameters
            hastv(p) && return true
        end
    end
    return false
end
mitv = filter(mi->hastv(mi.specTypes), nonspecialized)
println("There are ", length(mitv), " MethodInstances that have a TypeVar in specTypes")

# Let's analyze a specific case
m = which(similar, (Vector,))
_mi = Ref{Any}(nothing)
visit(m.specializations) do mi
    _mi[] isa Core.MethodInstance && return nothing
    hastv(mi.specTypes) && (_mi[] = mi)
    return nothing
end
mi = _mi[]
@assert(mi isa Core.MethodInstance)
println("\n## Performing an analysis of ", mi, '\n')
# Find the last callers with TypeVar specializations
callers = Core.MethodInstance[]
visit_backedges(mi) do caller
    hastv(caller.specTypes) || return false
    if !isdefined(caller, :backedges)
        push!(callers, caller)
        return false
    end
    foundone = false
    for edge in caller.backedges
        if !hastv(edge.specTypes)
            push!(callers, caller)
            foundone = true
            break
        end
    end
    return !foundone
end
# Let's look at the code of these callers
for caller in callers
    for be in caller.backedges
        f, t = call_type(be.specTypes)
        println("caller: ", be.def, " with argtypes ", t)
        code_warntype(f, t)
        println()
    end
end


text:
using DisplayAs
using Plots
using ToyStencils

us = ToyStencils.example()

heatmap(us; color = cgrad(:hawaii)) |> DisplayAs.PNG


text:
day18_input = [[[[[4, 6], 4], [1, 7]], [[[1, 6], [8, 4]], [1, 1]]], [[[[8, 5], [9, 2]], 1], [[2, 5], [[9, 4], [5, 9]]]], [[[[7, 3], 0], [8, 9]], 6], [[6, [[7, 2], [6, 2]]], [[[9, 8], 9], [9, 6]]], [2, [[[9, 2], 6], [[5, 3], [6, 7]]]], [[[5, [9, 6]], 0], [[[2, 8], [7, 0]], [7, [4, 4]]]], [[[[5, 0], 2], [0, 1]], 4], [2, [8, 8]], [[[[2, 5], [6, 8]], [[9, 8], 4]], [[[2, 3], [5, 8]], [9, 5]]], [[[[0, 7], [9, 4]], [[1, 0], 9]], [[[8, 8], [7, 2]], [3, [6, 5]]]], [[[[3, 2], 8], 1], [[4, [3, 4]], [[6, 5], [0, 6]]]], [[[7, 8], 8], [0, [5, 2]]], [[3, [3, 3]], [[[6, 9], [1, 1]], [6, [2, 9]]]], [[[[9, 7], [6, 8]], 4], [[[8, 2], [2, 9]], [8, [6, 2]]]], [[[[7, 3], 2], [9, 6]], [[[1, 7], [0, 0]], [4, 9]]], [[8, [7, [1, 0]]], 7], [[[7, [5, 1]], 0], [[8, [5, 3]], 4]], [1, [[[2, 6], 2], [1, [6, 0]]]], [[[5, 8], [[9, 1], 1]], [[3, [5, 0]], 5]], [[[[1, 5], [4, 9]], 8], [[7, 0], 6]], [9, [[0, [1, 0]], 6]], [[[[6, 8], 6], 9], [[7, 3], 2]], [[9, [[8, 7], 4]], [[[4, 0], [9, 0]], [8, 1]]], [[[2, [4, 4]], [7, [0, 1]]], [8, [[8, 6], [4, 0]]]], [0, 9], [[[[1, 8], [7, 4]], [[5, 0], [6, 1]]], [5, 7]], [[[[8, 2], [9, 2]], [8, [8, 4]]], [0, 4]], [[[[0, 7], [5, 8]], 3], 6], [[[7, [3, 4]], [3, [1, 5]]], 2], [[[1, [4, 2]], 5], [[1, 2], 1]], [[[[8, 2], [0, 9]], 1], [[[9, 0], [3, 5]], [8, [8, 0]]]], [[[0, 5], [1, [3, 3]]], [[[1, 0], [5, 2]], [7, 5]]], [[[4, [7, 3]], [0, 9]], [[2, 0], 8]], [[[[2, 2], 8], [7, 1]], 5], [[1, [[3, 8], 7]], [[7, [5, 8]], [4, [1, 7]]]], [[[[2, 7], 4], [8, [9, 1]]], [[5, 2], [4, 3]]], [[2, [7, 2]], [[8, [0, 8]], [0, [4, 2]]]], [[6, [6, [7, 4]]], [[7, [2, 0]], [[8, 2], 8]]], [[[7, [1, 7]], [[4, 1], 4]], [1, [4, 6]]], [1, [[1, 0], [[0, 3], [6, 9]]]], [[[[8, 6], 0], [[2, 8], [3, 0]]], [[[8, 2], 7], [[3, 0], 5]]], [[[[2, 8], 4], [2, [0, 7]]], [[3, [1, 2]], [[8, 0], [4, 2]]]], [1, 8], [[5, 6], 6], [[[[1, 0], [3, 6]], [[4, 0], 1]], [0, 7]], [[[5, [9, 6]], [7, [1, 2]]], 2], [[[6, 4], [[5, 6], [1, 8]]], [[[9, 0], [7, 7]], [[5, 8], [6, 8]]]], [8, 5], [5, [[[6, 8], 8], 0]], [[[[5, 7], [0, 0]], [6, [0, 0]]], [[[5, 5], 3], [[1, 1], [3, 4]]]], [[[4, 0], [[8, 6], 2]], [[3, [3, 1]], [[2, 8], [7, 2]]]], [[[8, 7], [[5, 5], [5, 3]]], 4], [[[[5, 4], 1], [3, 4]], [3, 5]], [[[6, 5], [[6, 3], 6]], 4], [[[[2, 2], [7, 1]], [6, 6]], [[8, [8, 7]], [[1, 6], [3, 0]]]], [[4, [[5, 0], [7, 4]]], [3, 1]], [[[3, [5, 8]], 5], [1, [[9, 6], 3]]], [[0, [[3, 0], [8, 7]]], [[1, 3], 3]], [5, [[3, [3, 3]], [3, 6]]], [[[[7, 3], 8], 3], [2, [[9, 8], 2]]], [[[2, 4], [[1, 2], 5]], [[[1, 2], [6, 0]], 3]], [[9, [[1, 1], [1, 7]]], [1, [2, [9, 1]]]], [[[5, [0, 0]], 5], [6, [0, 1]]], [[3, [[6, 5], 7]], [[7, 8], 3]], [[5, [2, 6]], 8], [[6, [0, [3, 0]]], [1, 2]], [3, [[[3, 7], 2], [[4, 0], 6]]], [[[8, [2, 7]], [4, 1]], [[2, [4, 2]], 3]], [[3, 2], [[[8, 8], [8, 6]], [[5, 3], 1]]], [1, [2, [[3, 2], [2, 9]]]], [8, [[9, 1], [[8, 4], [9, 9]]]], [[[4, [4, 6]], [1, 8]], [[7, 7], [[7, 4], 3]]], [[[8, 2], [[9, 7], [0, 8]]], [[4, 4], [[6, 1], 5]]], [[[3, [6, 6]], [[8, 6], [3, 7]]], [[7, 9], [[5, 3], 8]]], [[[8, 9], [8, 6]], [[[3, 3], [2, 9]], [[6, 6], 9]]], [8, [[[3, 0], 5], 2]], [[[[1, 3], 1], [[1, 9], 4]], [7, [3, 1]]], [[[[9, 3], 3], [[6, 8], 7]], [[[2, 0], 3], [8, [3, 6]]]], [[[[7, 1], [8, 1]], [[4, 6], [5, 9]]], [[[4, 5], 3], 5]], [6, [[3, [0, 0]], [6, 6]]], [[[[8, 8], [7, 6]], 3], [[[7, 7], [1, 1]], [[1, 8], [1, 4]]]], [[9, [8, [3, 4]]], [[6, 2], [1, 5]]], [[5, [3, 3]], [5, [0, 5]]], [[[[8, 9], 5], [1, 9]], [[5, [2, 8]], [[6, 4], [9, 4]]]], [2, 6], [[[[1, 4], 8], 5], [5, [0, [1, 7]]]], [[[[1, 0], [9, 9]], [0, 9]], [[[5, 4], [1, 6]], [9, [6, 7]]]], [[[7, 1], 5], [[3, 2], 5]], [[9, [[8, 8], [7, 0]]], [5, [3, [1, 3]]]], [[[[5, 2], [7, 5]], [4, [6, 7]]], [[[8, 1], 6], [2, [6, 6]]]], [[[5, 7], [6, [8, 2]]], [8, 2]], [[[[5, 7], 8], [[9, 8], 2]], [[2, 8], [[7, 6], 3]]], [[1, [[1, 6], 1]], [0, [[5, 9], [9, 1]]]], [[[[1, 4], [5, 0]], [[5, 5], [9, 3]]], [[6, 4], [4, [4, 6]]]], [7, [[5, [4, 8]], [[5, 9], 2]]], [[[[2, 9], [1, 8]], [4, 2]], 0], [[5, [[0, 9], [3, 7]]], [2, [6, [4, 8]]]], [[0, [5, 5]], 0], [[[5, 0], [[0, 5], 8]], [6, [[8, 7], [6, 5]]]], [[[5, [8, 2]], [8, 4]], [[6, 2], [8, [7, 0]]]]
]


text:
if basename(pwd()) == "aoc"
    cd("2017/17")
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
"""
    AbstractEddyViscosityClosure{TD, Iso} <: AbstractScalarDiffusivity{TD, Iso}

Abstract supertype for turbulence closures that are defined by an isotropic viscosity
and isotropic diffusivities.
"""
abstract type AbstractEddyViscosityClosure{TD, Iso} <: AbstractScalarDiffusivity{TD, Iso} end

@inline viscosity(closure::AbstractEddyViscosityClosure, diffusivities, args...) = diffusivities.νₑ
@inline diffusivity(::AbstractEddyViscosityClosure, ::Val{tracer_index}, diffusivities, args...) where tracer_index = diffusivities.κₑ[tracer_index]


text:
using MultiAssign
using Test

@testset "MultiAssign.jl" begin
    @multiassign a = 1
    @test a == 1
    @multiassign x, y, z = zeros(3)
    @test x == y == z == zeros(3)
end


text:
using Plots

@userplot orbitanim
@recipe function f!(oa::orbitanim)
    x, y = oa.args
    title --> "Trajectory"
    xaxis --> ("x", (-1.5,1.5))
    yaxis --> ("y", (-1.5, 1.5))
    markersize --> 2
    markercolor --> "red"
    seriestype --> :scatter
    legend --> :none
    framestyle --> :grid
    aspect_ratio --> 1
    dpi --> 300
    [x], [y]
end

function main()

    # IC
    topo = 10
    ω = 2π * topo / 60
    x0 = 1
    v0 = -x0*ω
    u0 = 1

    # equations
    x(t) = x0*cos(ω*t) + u0*sin(ω*t) / ω
    y(t) = v0*sin(ω*t) / ω

    # time frame
    t = 0:0.1:6
    n = length(t)

    xpos = x.(t)
    ypos = y.(t)

    # animation
    anim = @animate for i in 1:n

        # scatter!([xpos[i]], [ypos[i]], legend=:none, title="Trajectory",
        #         markersize=2, markercolor="blue",
        #         xaxis=("x", (-1.5,1.5)), yaxis=("y", (-1.5, 1.5)),
        #         framestyle=:grid, aspect_ratio=1, dpi=300)

        orbitanim!(xpos[i], ypos[i])
        scatter!((0, 0), markercolor="black")

    end

    gif(anim, "anim.gif", fps=15)

end

main()

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
export Dense,
       DenseTensor,
       array,
       matrix,
       vector,
       contract,
       outer,
       read,
       permutedims!!,
       write,
       ⊗

#
# Dense storage
#

struct Dense{ElT,VecT<:AbstractVector} <: TensorStorage{ElT}
  data::VecT
  function Dense{ElT,VecT}(data) where {ElT,VecT<:AbstractVector{ElT}}
    return new{ElT,VecT}(data)
  end
end

function Dense(data::VecT) where {VecT<:AbstractVector{ElT}} where {ElT}
  return Dense{ElT,VecT}(data)
end

function Dense{ElR}(data::AbstractVector{ElT}) where {ElR,ElT}
  ElT == ElR ? Dense(data) : Dense(ElR.(data))
end

Dense{ElT}(dim::Integer) where {ElT} = Dense(zeros(ElT,dim))

Dense{ElT}(::UndefInitializer,
           dim::Integer) where {ElT} = Dense(Vector{ElT}(undef,dim))

Dense(::Type{ElT},
      dim::Integer) where {ElT} = Dense{ElT}(dim)

Dense(x::ElT,
      dim::Integer) where {ElT<:Number} = Dense(fill(x,dim))

Dense(dim::Integer) = Dense(Float64,dim)

Dense(::Type{ElT},
      ::UndefInitializer,
      dim::Integer) where {ElT} = Dense{ElT}(undef,dim)

Dense(::UndefInitializer,dim::Integer) = Dense(Float64,undef,dim)

Dense{ElT}() where {ElT} = Dense(ElT[])
Dense(::Type{ElT}) where {ElT} = Dense{ElT}()

Base.copy(D::Dense) = Dense(copy(data(D)))

Base.complex(::Type{Dense{ElT,Vector{ElT}}}) where {ElT} = Dense{complex(ElT),Vector{complex(ElT)}}

Base.similar(D::Dense) = Dense(similar(data(D)))

Base.similar(D::Dense,length::Int) = Dense(similar(data(D),length))
Base.similar(::Type{<:Dense{ElT,VecT}},length::Int) where {ElT,VecT} = Dense(similar(VecT,length))

Base.similar(D::Dense,::Type{T}) where {T<:Number} = Dense(similar(data(D),T))

# TODO: should this do something different for SubArray?
Base.zeros(::Type{<:Dense{ElT}},dim::Int) where {ElT} = Dense{ElT}(dim)

function Base.promote_rule(::Type{<:Dense{ElT1,VecT1}},
                           ::Type{<:Dense{ElT2,VecT2}}) where {ElT1,VecT1,
                                                               ElT2,VecT2}
  ElR = promote_type(ElT1,ElT2)
  VecR = promote_type(VecT1,VecT2)
  return Dense{ElR,VecR}
end

# This is to get around the issue in Julia that:
# promote_type(Vector{ComplexF32},Vector{Float64}) == Vector{T} where T
function Base.promote_rule(::Type{<:Dense{ElT1,Vector{ElT1}}},
                           ::Type{<:Dense{ElT2,Vector{ElT2}}}) where {ElT1,ElT2}
  ElR = promote_type(ElT1,ElT2)
  VecR = Vector{ElR}
  return Dense{ElR,VecR}
end

# This is for type promotion for Scalar*Dense
function Base.promote_rule(::Type{<:Dense{ElT1,Vector{ElT1}}},
                           ::Type{ElT2}) where {ElT1,
                                                ElT2<:Number}
  ElR = promote_type(ElT1,ElT2)
  VecR = Vector{ElR}
  return Dense{ElR,VecR}
end

Base.convert(::Type{<:Dense{ElR,VecR}},
             D::Dense) where {ElR,VecR} = Dense(convert(VecR,data(D)))

# Make generic to all storage types?
Base.:*(D::Dense,x::Number) = Dense(x*data(D))
Base.:*(x::Number,D::Dense) = D*x

#
# DenseTensor (Tensor using Dense storage)
#

const DenseTensor{ElT,N,StoreT,IndsT} = Tensor{ElT,N,StoreT,IndsT} where {StoreT<:Dense}

DenseTensor(::Type{ElT},
            inds::Dims) where {ElT} = tensor(Dense(ElT,dim(inds)),inds)

DenseTensor(::Type{ElT},
            inds::Int...) where {ElT} = DenseTensor(ElT,inds)

DenseTensor(inds::Dims) = tensor(Dense(dim(inds)),inds)

DenseTensor(inds::Int...) = DenseTensor(inds)

DenseTensor(::Type{ElT},
            ::UndefInitializer,
            inds) where {ElT} = tensor(Dense(ElT,undef,dim(inds)),inds)

DenseTensor(::Type{ElT},
            ::UndefInitializer,
            inds::Int...) where {ElT} = DenseTensor(ElT,undef,inds)

DenseTensor(::UndefInitializer,
            inds::Dims) = tensor(Dense(undef,dim(inds)),inds)

DenseTensor(::UndefInitializer,
            inds::Int...) = DenseTensor(undef,inds)

# For convenience, direct Tensor constructors default to Dense
Tensor(::Type{ElT},inds::Dims) where {ElT} = DenseTensor(ElT,inds)
Tensor(::Type{ElT},inds::Int...) where {ElT} = DenseTensor(ElT,inds...)
Tensor(inds::Dims) = DenseTensor(inds)
Tensor(inds::Int...) = DenseTensor(inds...)

Tensor(::Type{ElT},
       ::UndefInitializer,
       inds::Dims) where {ElT} = DenseTensor(ElT,undef,inds)
Tensor(::Type{ElT},
       ::UndefInitializer,
       inds::Int...) where {ElT} = DenseTensor(ElT,undef,inds...)
Tensor(::UndefInitializer,
       inds::Dims) = DenseTensor(undef,inds)
Tensor(::UndefInitializer,
       inds::Int...) = DenseTensor(undef,inds...)

Tensor(A::Array{<:Number,N},
       inds::Dims{N}) where {N} = tensor(Dense(vec(A)),inds)

# Basic functionality for AbstractArray interface
Base.IndexStyle(::Type{<:DenseTensor}) = IndexLinear()

# Override CartesianIndices iteration to iterate
# linearly through the Dense storage (faster)
Base.iterate(T::DenseTensor,args...) = iterate(store(T),args...)

function Base.similar(::Type{<:DenseTensor{ElT}},
                      inds) where {ElT}
  return DenseTensor(ElT,undef,inds)
end

# To fix method ambiguity with similar(::AbstractArray,::Type)
function Base.similar(T::DenseTensor,
                      ::Type{ElT}) where {ElT}
  return tensor(similar(store(T),ElT), inds(T))
end

# To fix method ambiguity with similar(::AbstractArray,::Tuple)
function Base.similar(::Type{<:DenseTensor{ElT}},
                      inds::Dims) where {ElT}
  return DenseTensor(ElT,undef,inds)
end

Base.similar(T::DenseTensor,inds) = similar(typeof(T),inds)

# To fix method ambiguity with similar(::AbstractArray,::Tuple)
Base.similar(T::DenseTensor,inds::Dims) = similar(typeof(T),inds)

#
# Single index
#

Base.@propagate_inbounds function _getindex(T::DenseTensor{<:Number,N},
                                            I::CartesianIndex{N}) where {N}
  return store(T)[@inbounds LinearIndices(T)[CartesianIndex(I)]]
end

Base.@propagate_inbounds function Base.getindex(T::DenseTensor{<:Number,N},
                                                I::Vararg{Int,N}) where {N}
  return _getindex(T,CartesianIndex(I))
end

Base.@propagate_inbounds function Base.getindex(T::DenseTensor{<:Number,N},
                                                I::CartesianIndex{N}) where {N}
  return _getindex(T,I)
end

Base.@propagate_inbounds function _setindex!(T::DenseTensor{<:Number,N},
                                             x::Number,
                                             I::CartesianIndex{N}) where {N}
  store(T)[@inbounds LinearIndices(T)[CartesianIndex(I)]] = x
  return T
end

Base.@propagate_inbounds function Base.setindex!(T::DenseTensor{<:Number,N},
                                                 x::Number,
                                                 I::Vararg{Int,N}) where {N}
  _setindex!(T,x,CartesianIndex(I))
  return T
end

Base.@propagate_inbounds function Base.setindex!(T::DenseTensor{<:Number,N},
                                                 x::Number,
                                                 I::CartesianIndex{N}) where {N}
  _setindex!(T,x,I)
  return T
end

#
# Linear indexing
#

Base.@propagate_inbounds Base.getindex(T::DenseTensor,i::Int) = store(T)[i]

Base.@propagate_inbounds Base.setindex!(T::DenseTensor,v,i::Int) = (store(T)[i] = v; T)

#
# Slicing
# TODO: this doesn't allow colon right now
# Create a DenseView that stores a Dense and an offset?
#

Base.@propagate_inbounds function _getindex(T::DenseTensor{ElT,N},
                                            I::CartesianIndices{N}) where {ElT,N}
  storeR = Dense(vec(@view array(T)[I]))
  indsR = Tuple(I[end]-I[1]+CartesianIndex(ntuple(_->1,Val(N))))
  return tensor(storeR, indsR)
end

Base.@propagate_inbounds function Base.getindex(T::DenseTensor{ElT,N},
                                                I...) where {ElT,N}
  return _getindex(T,CartesianIndices(I))
end

# Reshape a DenseTensor using the specified dimensions
# This returns a view into the same Tensor data
function Base.reshape(T::DenseTensor, dims)
  dim(T)==dim(dims) || error("Total new dimension must be the same as the old dimension")
  return tensor(store(T), dims)
end

# This version fixes method ambiguity with AbstractArray reshape
function Base.reshape(T::DenseTensor, dims::Dims)
  dim(T)==dim(dims) || error("Total new dimension must be the same as the old dimension")
  return tensor(store(T), dims)
end

function Base.reshape(T::DenseTensor, dims::Int...)
  return tensor(store(T), tuple(dims...))
end

# Create an Array that is a view of the Dense Tensor
# Useful for using Base Array functions
array(T::DenseTensor) = reshape(data(store(T)),dims(inds(T)))

function Base.Array{ElT,N}(T::DenseTensor{ElT,N}) where {ElT,N}
  return copy(array(T))
end

function Base.Array(T::DenseTensor{ElT,N}) where {ElT,N}
  return Array{ElT,N}(T)
end

# TODO: call permutedims!(R,T,perm,(r,t)->t)?
function Base.permutedims!(R::DenseTensor{<:Number,N},
                           T::DenseTensor{<:Number,N},
                           perm::NTuple{N,Int}) where {N}
  permutedims!(array(R),array(T),perm)
  return R
end

function apply!(R::DenseTensor,
                T::DenseTensor,
                f::Function=(r,t)->t)
  RA = array(R)
  TA = array(T)
  RA .= f.(RA,TA)
  return R
end

# Version that may overwrite the result or promote
# and return the result
# TODO: move to tensor.jl?
function permutedims!!(R::Tensor,
                       T::Tensor,
                       perm::NTuple{N,Int},
                       f::Function=(r,t)->t) where {N}
  if !is_trivial_permutation(perm)
    R = permutedims!(R,T,perm,f)
  else
    R = apply!(R,T,f)
  end
  return R
end

# TODO: move to tensor.jl?
function Base.permutedims(T::Tensor{<:Number,N},
                          perm::NTuple{N,Int}) where {N}
  Tp = similar(T,permute(inds(T),perm))
  Tp = permutedims!!(Tp,T,perm)
  return Tp
end

Base.permutedims(T::Tensor,
                 perm::Tuple{Vararg{Int}}) = error("Permutation size must match tensor order")

# TODO: move to tensor.jl?
function Base.:*(x::Number,
                 T::Tensor)
  return tensor(x*store(T),inds(T))
end
Base.:*(T::Tensor, x::Number) = x*T

# For use in custom permutedims!
using Base.Cartesian: @nexprs,
                      @ntuple,
                      @nloops

#
# A generalized permutedims!(P,B,perm) that also allows
# a function to be applied elementwise
# TODO: benchmark to make sure it is similar to Base.permutedims!
#
# Based off of the permutedims! implementation in Julia's base:
# https://github.com/JuliaLang/julia/blob/91151ab871c7e7d6689d1cfa793c12062d37d6b6/base/multidimensional.jl#L1355
#
@generated function Base.permutedims!(TTP::DenseTensor{<:Number,N},
                                      TT::DenseTensor{<:Number,N},
                                      perm,
                                      f::Function) where {N}
  quote
    TP = array(TTP)
    T = array(TT)
    Base.checkdims_perm(TP, T, perm)

    #calculates all the strides
    native_strides = Base.size_to_strides(1, size(T)...)
    strides_1 = 0
    @nexprs $N d->(strides_{d+1} = native_strides[perm[d]])

    #Creates offset, because indexing starts at 1
    offset = 1 - sum(@ntuple $N d->strides_{d+1})

    ind = 1
    @nexprs 1 d->(counts_{$N+1} = strides_{$N+1}) # a trick to set counts_($N+1)
    @nloops($N, i, TP,
            d->(counts_d = strides_d), # PRE
            d->(counts_{d+1} += strides_{d+1}), # POST
            begin # BODY
                sumc = sum(@ntuple $N d->counts_{d+1})
                @inbounds TP[ind] = f(TP[ind],T[sumc+offset])
                ind += 1
            end)

    return TTP
  end
end
function Base.permutedims!(TP::DenseTensor{<:Number,0},
                           T::DenseTensor{<:Number,0},
                           perm,
                           f::Function)
  TP[] = f(TP[],T[])
  return TP
end

function outer!(R::DenseTensor,
                T1::DenseTensor,
                T2::DenseTensor)
  v1 = vec(T1)
  v2 = vec(T2)
  RM = reshape(R,dim(v1),dim(v2))
  RM .= v1 .* transpose(v2)
  return R
end

function outer!!(R::Tensor,
                 T1::Tensor,
                 T2::Tensor)
  outer!(R,T1,T2)
  return R
end

# TODO: call outer!!, make this generic
function outer(T1::DenseTensor{ElT1},
               T2::DenseTensor{ElT2}) where {ElT1,ElT2}
  array_outer = vec(array(T1))*transpose(vec(array(T2)))
  inds_outer = unioninds(inds(T1),inds(T2))
  return tensor(Dense{promote_type(ElT1,ElT2)}(vec(array_outer)),
                inds_outer)
end
const ⊗ = outer

# TODO: move to tensor.jl?
function contraction_output_type(::Type{TensorT1},
                                 ::Type{TensorT2},
                                 ::Type{IndsR}) where {TensorT1<:Tensor,
                                                       TensorT2<:Tensor,
                                                       IndsR}
  return similar_type(promote_type(TensorT1,TensorT2),IndsR)
end

function contraction_output(::TensorT1,
                            ::TensorT2,
                            indsR::IndsR) where {TensorT1<:DenseTensor,
                                                 TensorT2<:DenseTensor,
                                                 IndsR}
  TensorR = contraction_output_type(TensorT1,TensorT2,IndsR)
  return similar(TensorR,indsR)
end

# TODO: move to tensor.jl?
function contraction_output(T1::Tensor,
                            labelsT1,
                            T2::Tensor,
                            labelsT2,
                            labelsR)
  indsR = contract_inds(inds(T1),labelsT1,inds(T2),labelsT2,labelsR)
  R = contraction_output(T1,T2,indsR)
  return R
end

# TODO: move to tensor.jl?
function contract(T1::Tensor{<:Any,N1},
                  labelsT1,
                  T2::Tensor{<:Any,N2},
                  labelsT2,
                  labelsR = contract_labels(labelsT1,labelsT2)) where {N1,N2}
  # TODO: put the contract_inds logic into contraction_output,
  # call like R = contraction_ouput(T1,labelsT1,T2,labelsT2)
  #indsR = contract_inds(inds(T1),labelsT1,inds(T2),labelsT2,labelsR)
  R = contraction_output(T1,labelsT1,T2,labelsT2,labelsR)
  # contract!! version here since the output R may not
  # be mutable (like UniformDiag)
  R = contract!!(R,labelsR,T1,labelsT1,T2,labelsT2)
  return R
end

# Move to tensor.jl? Is this generic for all storage types?
function contract!!(R::Tensor{<:Number,NR},
                    labelsR::NTuple{NR},
                    T1::Tensor{<:Number,N1},
                    labelsT1::NTuple{N1},
                    T2::Tensor{<:Number,N2},
                    labelsT2::NTuple{N2},
                    α::Number=1, β::Number=0) where {NR,N1,N2}
  if N1==0
    (α ≠ 1 || β ≠ 0) &&
      error("contract!! not yet implemented for scalar ITensor with non-trivial α and β")
    # TODO: replace with an add! function?
    # What about doing `R .= T1[] .* PermutedDimsArray(T2,perm)`?
    perm = getperm(labelsR,labelsT2)
    R = permutedims!!(R,T2,perm,(r,t2)->T1[]*t2)
  elseif N2==0
    (α ≠ 1 || β ≠ 0) &&
      error("contract!! not yet implemented for scalar ITensor with non-trivial α and β")
    perm = getperm(labelsR,labelsT1)
    R = permutedims!!(R,T1,perm,(r,t1)->T2[]*t1)
  elseif N1+N2==NR
    (α ≠ 1 || β ≠ 0) &&
      error("contract!! not yet implemented for outer product tensor contraction with non-trivial α and β")
    # TODO: permute T1 and T2 appropriately first (can be more efficient
    # then permuting the result of T1⊗T2)
    # TODO: implement the in-place version directly
    R = outer!!(R,T1,T2)
    labelsRp = (labelsT1..., labelsT2...)
    perm = getperm(labelsR, labelsRp)
    if !is_trivial_permutation(perm)
      R = permutedims!!(R, copy(R), perm)
    end
  else
    if α ≠ 1 || β ≠ 0
      R = _contract!!(R, labelsR,
                      T1, labelsT1,
                      T2, labelsT2,
                      α, β)
    else
      R = _contract!!(R, labelsR,
                      T1, labelsT1,
                      T2, labelsT2)
    end
  end
  return R
end

# Move to tensor.jl? Overload this function
# for immutable storage types
function _contract!!(R::Tensor, labelsR,
                     T1::Tensor, labelsT1,
                     T2::Tensor, labelsT2,
                     α::Number=1, β::Number=0)
  if α ≠ 1 || β ≠ 0
    contract!(R, labelsR,
              T1, labelsT1,
              T2, labelsT2,
              α, β)
  else
    contract!(R, labelsR,
              T1, labelsT1,
              T2, labelsT2)
  end
  return R
end

function contract!(R::DenseTensor{<:Number,NR},
                   labelsR,
                   T1::DenseTensor{ElT1,N1},
                   labelsT1,
                   T2::DenseTensor{ElT2,N2},
                   labelsT2,
                   α::Number=1,
                   β::Number=0) where {ElT1,ElT2,N1,N2,NR}
  if N1+N2==NR
    outer!(R,T1,T2)
    labelsRp = (labelsT1..., labelsT2...)
    perm = getperm(labelsR,labelsRp)
    if !is_trivial_permutation(perm)
      permutedims!(R, copy(R), perm)
    end
    return R
  end

  props = ContractionProperties(labelsT1,labelsT2,labelsR)
  compute_contraction_properties!(props,T1,T2,R)

  if ElT1 != ElT2
    # TODO: use promote instead
    # T1,T2 = promote(T1,T2)

    ElR = promote_type(ElT1,ElT2)
    if ElT1 != ElR
      # TODO: get this working
      # T1 = ElR.(T1)
      T1 = one(ElR) * T1
    end
    if ElT2 != ElR
      # TODO: get this working
      # T2 = ElR.(T2)
      T2 = one(ElR) * T2
    end
  end

  _contract!(R,T1,T2,props,α,β)
  return R
end

function _contract!(CT::DenseTensor{El,NC},
                    AT::DenseTensor{El,NA},
                    BT::DenseTensor{El,NB},
                    props::ContractionProperties,
                    α::Number=one(El),
                    β::Number=zero(El)) where {El,NC,NA,NB}
  # TODO: directly use Tensor instead of Array
  C = array(CT)
  A = array(AT)
  B = array(BT)

  tA = 'N'
  if props.permuteA
    AM = reshape(permutedims(A,NTuple{NA,Int}(props.PA)),
                 props.dmid, props.dleft)
    tA = 'T'
  else
    #A doesn't have to be permuted
    if Atrans(props)
      AM = reshape(A,props.dmid,props.dleft)
      tA = 'T'
    else
      AM = reshape(A,props.dleft,props.dmid)
    end
  end

  tB = 'N'
  if props.permuteB
    BM = reshape(permutedims(B,NTuple{NB,Int}(props.PB)),
                 props.dmid, props.dright)
  else
    if Btrans(props)
      BM = reshape(B,props.dright,props.dmid)
      tB = 'T'
    else
      BM = reshape(B,props.dmid,props.dright)
    end
  end

  # TODO: this logic may be wrong
  if props.permuteC
    # Need to copy here since we will be permuting
    # into C later
    CM = reshape(copy(C), props.dleft, props.dright)
  else
    if Ctrans(props)
      CM = reshape(C,props.dleft,props.dright)
      (AM,BM) = (BM,AM)
      if tA==tB
        tA = tB = (tA == 'T' ? 'N' : 'T')
      end
    else
      CM = reshape(C,props.dleft,props.dright)
    end
  end

  BLAS.gemm!(tA,tB,El(α),AM,BM,El(β),CM)

  if props.permuteC
    permutedims!(C,reshape(CM,props.newCrange...),
                 NTuple{NC,Int}(props.PC))
  end
  return C
end

"""
permute_reshape(T::Tensor,pos...)

Takes a permutation that is split up into tuples. Index positions
within the tuples are combined.

For example:

permute_reshape(T,(3,2),1)

First T is permuted as `permutedims(3,2,1)`, then reshaped such
that the original indices 3 and 2 are combined.
"""
function permute_reshape(T::DenseTensor{ElT,NT,IndsT},
                         pos::Vararg{<:Any,N}) where {ElT,NT,IndsT,N}
  perm = flatten(pos...)

  length(perm)≠NT && error("Index positions must add up to order of Tensor ($N)")
  isperm(perm) || error("Index positions must be a permutation")

  dimsT = dims(T)
  indsT = inds(T)
  if !is_trivial_permutation(perm)
    T = permutedims(T,perm)
  end
  N==NT && return T
  newdims = MVector(ntuple(_->eltype(IndsT)(1),Val(N)))
  for i ∈ 1:N
    if length(pos[i])==1
      # No reshape needed, just use the
      # original index
      newdims[i] = indsT[pos[i][1]]
    else
      newdim_i = 1
      for p ∈ pos[i]
        newdim_i *= dimsT[p]
      end
      newdims[i] = eltype(IndsT)(newdim_i)
    end
  end
  newinds = similar_type(IndsT,Val{N})(Tuple(newdims))
  return reshape(T,newinds)
end

# svd of an order-n tensor according to positions Lpos
# and Rpos
function LinearAlgebra.svd(T::DenseTensor{<:Number,N,IndsT},
                           Lpos::NTuple{NL,Int},
                           Rpos::NTuple{NR,Int};
                           kwargs...) where {N,IndsT,NL,NR}
  M = permute_reshape(T,Lpos,Rpos)
  UM,S,VM,spec = svd(M;kwargs...)
  u = ind(UM,2)
  v = ind(VM,2)

  Linds = similar_type(IndsT,Val{NL})(ntuple(i->inds(T)[Lpos[i]],Val(NL)))
  Uinds = push(Linds,u)

  # TODO: do these positions need to be reversed?
  Rinds = similar_type(IndsT,Val{NR})(ntuple(i->inds(T)[Rpos[i]],Val(NR)))
  Vinds = push(Rinds,v)

  U = reshape(UM,Uinds)
  V = reshape(VM,Vinds)

  return U,S,V,spec
end

# qr decomposition of an order-n tensor according to
# positions Lpos and Rpos
function LinearAlgebra.qr(T::DenseTensor{<:Number,N,IndsT},
                          Lpos::NTuple{NL,Int},
                          Rpos::NTuple{NR,Int};kwargs...) where {N,IndsT,NL,NR}
  M = permute_reshape(T,Lpos,Rpos)
  QM,RM = qr(M;kwargs...)
  q = ind(QM,2)
  r = ind(RM,1)
  # TODO: simplify this by permuting inds(T) by (Lpos,Rpos)
  # then grab Linds,Rinds
  Linds = similar_type(IndsT,Val{NL})(ntuple(i->inds(T)[Lpos[i]],Val(NL)))
  Qinds = push(Linds,r)
  Q = reshape(QM,Qinds)
  Rinds = similar_type(IndsT,Val{NR})(ntuple(i->inds(T)[Rpos[i]],Val(NR)))
  Rinds = pushfirst(Rinds,r)
  R = reshape(RM,Rinds)
  return Q,R
end

# polar decomposition of an order-n tensor according to positions Lpos
# and Rpos
function polar(T::DenseTensor{<:Number,N,IndsT},
               Lpos::NTuple{NL,Int},
               Rpos::NTuple{NR,Int}) where {N,IndsT,NL,NR}
  M = permute_reshape(T,Lpos,Rpos)
  UM,PM = polar(M)

  # TODO: turn these into functions
  Linds = similar_type(IndsT,Val{NL})(ntuple(i->inds(T)[Lpos[i]],Val(NL)))
  Rinds = similar_type(IndsT,Val{NR})(ntuple(i->inds(T)[Rpos[i]],Val(NR)))

  # Use sim to create "similar" indices, in case
  # the indices have identifiers. If not this should
  # act as an identity operator
  simRinds = sim(Rinds)
  Uinds = (Linds..., simRinds...)
  Pinds = (simRinds..., Rinds...)

  U = reshape(UM,Uinds)
  P = reshape(PM,Pinds)
  return U,P
end

function LinearAlgebra.exp(T::DenseTensor{ElT,N},
                           Lpos::NTuple{NL,Int},
                           Rpos::NTuple{NR,Int};
                           ishermitian::Bool=false) where {ElT,N,
                                                           NL,NR}
  M = permute_reshape(T,Lpos,Rpos)
  indsTp = permute(inds(T), (Lpos...,Rpos...))
  if ishermitian
    expM = exp(Hermitian(matrix(M)))
    return tensor(Dense{ElT}(vec(expM)), indsTp)
  else
    expM = exp(M)
    return reshape(expM,indsTp)
  end
end

function HDF5.write(parent::Union{HDF5File, HDF5Group},
                    name::String,
                    D::Store) where {Store <: Dense}
  g = g_create(parent,name)
  attrs(g)["type"] = "Dense{$(eltype(Store))}"
  attrs(g)["version"] = 1
  if eltype(D) != Nothing
    write(g,"data",D.data)
  end
end

function HDF5.read(parent::Union{HDF5File,HDF5Group},
                   name::AbstractString,
                   ::Type{Store}) where {Store <: Dense}
  g = g_open(parent,name)
  ElT = eltype(Store)
  typestr = "Dense{$ElT}"
  if read(attrs(g)["type"]) != typestr
    error("HDF5 group or file does not contain $typestr data")
  end
  if ElT == Nothing
    return Dense{Nothing}()
  end
  data = read(g,"data")
  return Dense{ElT}(data)
end

function Base.show(io::IO,
                   mime::MIME"text/plain",
                   T::DenseTensor)
  summary(io,T)
  println(io)
  print_tensor(io,T)
  println(io)
end


text:
module ArgoData

include("MITprof.jl")
include("DownloadArgo.jl")

export MITprof, DownloadArgo

end # module


text:
using Gen
import Random

@gen function model()
    if ({:z} ~ bernoulli(0.5))
        m1 = ({:m1} ~ gamma(1, 1))
        m2 = ({:m2} ~ gamma(1, 1))
    else
        m = ({:m} ~ gamma(1, 1))
        (m1, m2) = (m, m)
    end
    {:y1} ~ normal(m1, 0.1)
    {:y2} ~ normal(m2, 0.1)
end

@gen function mean_random_walk_proposal(trace)
    if trace[:z]
        {:m1} ~ normal(trace[:m1], 0.1)
        {:m2} ~ normal(trace[:m2], 0.1)
    else
        {:m} ~ normal(trace[:m], 0.1)
    end
end

@gen function split_merge_proposal(trace)
    if trace[:z]
        # currently two means, switch to one
    else
        # currently one mean, switch to two
        {:u} ~ uniform_continuous(0, 1)
    end
end

function merge_mean(m1, m2)
    m = sqrt(m1 * m2)
    u = m1 / (m1 + m2)
    (m, u)
end

function split_mean(m, u)
    m1 = m * sqrt((u / (1 - u)))
    m2 = m * sqrt(((1 - u) / u))
    (m1, m2)
end

@transform split_merge_involution (model_in, aux_in) to (model_out, aux_out) begin
    if @read(model_in[:z], :discrete)
        # currently two means, switch to one
        @write(model_out[:z], false, :discrete)
        m1 = @read(model_in[:m1], :continuous)
        m2 = @read(model_in[:m2], :continuous)
        (m, u) = merge_mean(m1, m2)
        @write(model_out[:m], m, :continuous)
        @write(aux_out[:u], u, :continuous)
    else
        # currently one mean, switch to two
        @write(model_out[:z], true, :discrete)
        m = @read(model_in[:m], :continuous)
        u = @read(aux_in[:u], :continuous)
        (m1, m2) = split_mean(m, u)
        @write(model_out[:m1], m1, :continuous)
        @write(model_out[:m2], m2, :continuous)
    end
end

is_involution!(split_merge_involution)

function do_inference_simple(y1, y2)
    trace, = generate(model, (), choicemap((:y1, y1), (:y2, y2), (:z, false), (:m, 1.2)))
    zs = Bool[]
    m = Float64[]
    m1 = Float64[]
    m2 = Float64[]
    for iter=1:100
        trace, = mh(trace, select(:z))
        trace, = mh(trace, mean_random_walk_proposal, ())
        push!(zs, trace[:z])
        push!(m, trace[:z] ? NaN : trace[:m])
        push!(m1, trace[:z] ? trace[:m1] : NaN)
        push!(m2, trace[:z] ? trace[:m2] : NaN)
    end
    (zs, m, m1, m2)
end

function do_inference_rjmcmc(y1, y2)
    trace, = generate(model, (), choicemap((:y1, y1), (:y2, y2), (:z, false), (:m, 1.2)))
    zs = Bool[]
    m = Float64[]
    m1 = Float64[]
    m2 = Float64[]
    for iter=1:100
        trace, = mh(
            trace, split_merge_proposal, (), split_merge_involution; check=true)
        trace, = mh(trace, mean_random_walk_proposal, ())
        push!(zs, trace[:z])
        push!(m, trace[:z] ? NaN : trace[:m])
        push!(m1, trace[:z] ? trace[:m1] : NaN)
        push!(m2, trace[:z] ? trace[:m2] : NaN)
    end
    (zs, m, m1, m2)
end

y1, y2 = (1.0, 1.3)
(zs, m, m1, m2) = do_inference_rjmcmc(y1, y2)
println(zs)

using Plots

function make_plots()
    Random.seed!(2)

    y1, y2 = (1.0, 1.3)
    (zs, m, m1, m2) = do_inference_rjmcmc(y1, y2)
    p1 = plot(title="Involution MH (RJMCMC)", m, label="m")
    plot!(m1, label="m1")
    plot!(m2, label="m2")
    ylims!(0.5, 1.5)

    p2 = plot(zs, label="z", color="black")
    xlabel!("# MCMC moves")
    yticks!([0, 1], ["F", "T"])
    ylims!(-0.1, 1.1)

    y1, y2 = (1.0, 1.3)
    (zs, m, m1, m2) = do_inference_simple(y1, y2)
    p3 = plot(title="Selection MH", m, label="m")
    plot!(m1, label="m1")
    plot!(m2, label="m2")
    ylims!(0.5, 1.5)

    p4 = plot(zs, label="z", color="black")
    xlabel!("# MCMC moves")
    yticks!([0, 1], ["F", "T"])
    ylims!(-0.1, 1.1)

    plot(p1, p3, p2, p4)
    savefig("rjmcmc.png")
end

make_plots()

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
using ExecutableSpecifications:
    ExecutorEngine, ColorConsolePresenter, Driver,
    readstepdefinitions!, runfeatures!, issuccess,
    FromSourceExecutionEnvironment, NoExecutionEnvironment
import ExecutableSpecifications:
    findfileswithextension, readfile, fileexists

struct OSAL <: ExecutableSpecifications.OSAbstraction end
function findfileswithextension(::OSAL, path::String, extension::String)
    [joinpath(path, filename) for filename in readdir(path)
              if isfile(joinpath(path, filename)) &&
                 splitext(joinpath(path, filename))[2] == extension]
end

readfile(::OSAL, path::String) = read(path, String)
fileexists(::OSAL, path::String) = isfile(path)

"""
    runspec()

Execute all features found from the current directory, or another specified directory.
"""
function runspec(rootpath::String = ".")
    featurepath = joinpath(rootpath, "features")
    stepspath = joinpath(featurepath, "steps")
    execenvpath = joinpath(featurepath, "environment.jl")
    os = OSAL()

    executionenv = if fileexists(os, execenvpath)
        FromSourceExecutionEnvironment(readfile(os, execenvpath))
    else
        NoExecutionEnvironment()
    end

    engine = ExecutorEngine(ColorConsolePresenter(); executionenv=executionenv)
    driver = Driver(os, engine)

    readstepdefinitions!(driver, stepspath)
    resultaccumulator = runfeatures!(driver, featurepath)

    #
    # Present number of scenarios that succeeded and failed for each feature
    #
    results = featureresults(resultaccumulator)

    # Find the longest feature name, so we can align the result table.
    maxfeature = maximum(length(r.feature.header.description) for r in results)

    featureprefix = "  Feature: "
    printstyled(" " ^ (length(featureprefix) + maxfeature + 1), "| Success | Failure\n"; color=:white)
    for r in results
        linecolor = r.n_failure == 0 ? :green : :red
        printstyled(featureprefix, rpad(r.feature.header.description, maxfeature); color=linecolor)
        printstyled(" | "; color=:white)
        printstyled(rpad("$(r.n_success)", 7); color=:green)
        printstyled(" | "; color=:white)
        printstyled(rpad("$(r.n_failure)", 7), "\n"; color=linecolor)
    end

    println()

    istotalsuccess = issuccess(resultaccumulator)
    if istotalsuccess
        println("SUCCESS")
    else
        println("FAILURE")
    end

    istotalsuccess
end


text:
using Test
using Distributions
using PkgA
using Random

@testset "randchisq1" begin
  Random.seed!(0)
  x = [PkgA.randchisq1() for _ in 1:1_000_000]
  @test isapprox(mean(x), 1, rtol=0.01)
  @test isapprox(var(x), 2, rtol=0.01)
end


text:
__precompile__()

module Requires

include("init.jl")
include("getthing.jl")
include("require.jl")

end # module


text:
export OSProc, Context

abstract type Processor end

"""
OS process - contains pid returned by `addprocs`
"""
struct OSProc <: Processor
    pid::Int
end

"""
A context represents a set of processors to use
for a papply operation.
"""
mutable struct Context
    procs::Array{OSProc}
    log_sink::Any
    profile::Bool
    options
end

function Context(xs)
    Context(xs, NoOpLog(), false, nothing) # By default don't log events
end
Context(xs::Array{Int}) = Context(map(OSProc, xs))
Context(;nthreads=Threads.nthreads()) = Context([workers() for i=1:nthreads] |> Iterators.flatten |> collect)
procs(c::Context) = c.procs

"""
Write a log event
"""
function write_event(ctx::Context, event::Event)
    write_event(ctx.log_sink, event)
end


text:
module MaxHelpingHandTempleEyeTrackingMadeline

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Temple Eye (Small, Follow Madeline) (max480's Helping Hand)" => Ahorn.EntityPlacement(
        Maple.TempleEye,
        "point",
        Dict{String, Any}(
            "followMadeline" => true
        )
    )
)

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
using Revise
using JDF

a=gf()
b = Array(a.ORIGIN_STATE_ABR)
io = iow()
@which compress_then_write(eltype(b), b, io)
@time metadata = compress_then_write(eltype(b), b, io)
close(io)

io = ior()
@which column_loader(eltype(b), b, io)

@time okhl = column_loader(eltype(b), io, metadata)
close(io)

all(skipmissing(b) .== skipmissing(okhl))


text:
using Documenter

makedocs(
    modules = [polyline],
    format = :html,
    sitename = "polyline.jl",
    pages = Any[
        "Contents" => "contents.md",
        "index.md",
    ],
)

deploydocs(
    repo = "github.com/NikStoyanov/polyline.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing,
)


text:
#=
  https://edu.swi-prolog.org/mod/assign/view.php?id=254&forceview=1
  """
  Coin tosses

  http://cplint.eu/p/coin_tosses.swinb

  Coin tosses

  Consider a process where you repeatedly toss coins and record the results.

  After each toss, you continue tossing with probability 0.8.

  Write a predicate tosses(Coin,L) that, given an initial coin id Coin, returns in L the list of
  results of coin tosses obtained using the process above.

  Moreover, write a predicate length_tosses(N) that returns the number of coin tosses.

  Compute the probability of the sequence [h,h,h] using MCINTYRE.

  Compute the probability of the sequence [h,h,h,h,h] given that the subsequence [h,h,h] was observed
  using rejection sampling and Metropolis Hastings.

  Compute the probability of the sequences of 10 coin tosses using MCINTYRE.

  Compute the expected length of the sequences of coin tosses using MCINTYRE.
  """

  Note: The cplint version (corrected by  Fabrizio Riguzzi) at http://cplint.eu/p/coin_tosses_hakank_rzf.swinb
  give another solutions:
   - first problem: probability of the sequence [h,h,h] -> 0.015
   - second problem: probability of the sequence [h,h,h,h,h] given that the subsequence [h,h,h] -> about 0.03
   - third problem: probability of the sequences of 10 coin tosses -> about 0.026
   - fourth problem: expected length of the sequences of coin tosses -> about 4.

  Cf ~/blog/coin_tosses.blog
     ~/webppl/coin_tosses.wppl

   This is a recursive approach.

=#

#=

# using
#   chains = sample(model, MH(), MCMCThreads(), 10_000, num_chains)
#
Julia> @time include("coin_tosses2.jl")
Problem 1
Probability of problem 1:

Distributions of variable (num:20)
0.00000 =>   37429  (0.935725)
1.00000 =>    2571  (0.064275)

Problem 2
Probability of problem 2:

Distributions of variable (num:20)
0.00000 =>   36836  (0.920900)
1.00000 =>    3164  (0.079100)

Problem 3
Probability of problem 3:

Distributions of variable (num:20)
0.00000 =>   38931  (0.973275)
1.00000 =>    1069  (0.026725)

Problem 4
Probability of problem 4:

Distributions of variable (num:20)
1.00000 =>    7959  (0.198975)
2.00000 =>    6391  (0.159775)
3.00000 =>    5168  (0.129200)
4.00000 =>    4184  (0.104600)
5.00000 =>    3259  (0.081475)
6.00000 =>    2610  (0.065250)
7.00000 =>    2083  (0.052075)
8.00000 =>    1661  (0.041525)
9.00000 =>    1329  (0.033225)
10.00000 =>    1081  (0.027025)
11.00000 =>     872  (0.021800)
12.00000 =>     711  (0.017775)
13.00000 =>     525  (0.013125)
14.00000 =>     422  (0.010550)
15.00000 =>     373  (0.009325)
16.00000 =>     301  (0.007525)
17.00000 =>     211  (0.005275)
18.00000 =>     175  (0.004375)
19.00000 =>     155  (0.003875)
20.00000 =>      99  (0.002475)
Mean:4.9843

  4.923102 seconds (71.82 M allocations: 3.655 GiB, 11.41% gc time)

=#
using Turing, StatsPlots, DataFrames, StatsBase
include("jl_utils.jl")

@model function coin_tosses(problem)
    head = 1
    tail = 2

    len ~ DiscreteUniform(1,100)
    # We create a fixed length list which seems to be much easier
    # than dynamically adding elements to tosses, especially
    # since for problem 2 we always access to tosses[1:5]
    function toss(a)
        if length(a) == 0
            # Alway do the first throw
            t ~ Categorical([0.5,0.5])
            return toss(vcat(a, t))
        else
            if problem == 2 && length(a) <= 3
                # For problem 2 we know that there are 3 heads
                return toss(vcat(a, head))
            else
                # For other problems
                if rand(flip(0.8))
                    t ~ Categorical([0.5,0.5])
                    return toss(vcat(a, t))
                else
                    return a
                end
            end
        end
    end

    tosses = toss(Int8[])
    len = length(tosses)

    # * first problem: probability of the sequence [h,h,h] -> 0.015
    # * second problem: probability of the sequence [h,h,h,h,h] given that the subsequence [h,h,h] -> about 0.03
    # * third problem: probability of the sequences of 10 coin tosses -> about 0.026
    # * fourth problem: expected length of the sequences of coin tosses -> about 4.
    if problem == 1
        return len == 3 && tosses[1] == head && tosses[2] == head && tosses[3] == head
    elseif problem == 2
        return len == 5 && tosses[1] == head && tosses[2] == head && tosses[3] == head &&
                           tosses[4] == head && tosses[5] == head
    elseif problem == 3
        return len == 10
    else
        return len
    end
end

function run_problem(problem)
    println("Problem $problem")
    model = coin_tosses(problem)
    num_chains = 4

    # chains = sample(model, MH(), MCMCThreads(), 100_000, num_chains)
    # chains = sample(model, MH(), MCMCThreads(), 10_000, num_chains)
    # chains = sample(model, MH(), 10_000)

    # chains = sample(model, PG(15), MCMCThreads(), 1_000, num_chains)

    # chains = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chains)
    # chains = sample(model, SMC(1000), 10_000)
    chains = sample(model, IS(), 10_000)

    display(chains)
    println("avg length: $(mean(chains[:len]))")
    println("fit(length): ")
    cc = round.(Int64,vcat(chains[:len]...)) # Convert to integer
    d = fit_mle(Geometric{Int64}, cc)
    println("d: $d")
    # println("avg t: $(mean(chains[:t]))")
    println("Probability of problem $problem:")
    genq = generated_quantities(model,chains)

    show_var_dist_pct(genq,100)
    if problem == 4
        println("Mean:",mean(genq))
    end
end

for problem in 1:4
    run_problem(problem)
    println("\n")
end

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using SymbolicUtils
using SymbolicUtils: Sym, Term, FnType, Symbolic
using SymbolicCodegen: Sum

# get it?
function gensum(t,i,a,b)
    new_i = Sym{Int}(gensym(:i))
    new_t = substitute(t, Dict(i => new_i))
    return Sum(new_t, new_i, a, b)
end

using SymbolicUtils
const RW = SymbolicUtils.Rewriters

POSTRULES = [
    @rule getindex(UnitRange(~a,~b), ~i) => ~i - ~a + 1
    @rule (+(~~x))^2 => sum(a*b for a in (~~x) for b in (~~x))
]

PRERULES = [
    @acrule (~a + ~b)*(~c) => (~a) * (~c) + (~b) * (~c)
    @rule Sum(+(~~x), ~i, ~a, ~b) => sum([gensum(t, ~i, ~a, ~b) for t in (~~x)])
    @rule Sum(*(~~x), ~i, ~a, ~b) => SymbolicCodegen.tryfactor(~~x, ~i, ~a, ~b) # ifelse(!_contains(~x,~i) || !_contains(~y,~i), Sum(~x, ~i, ~a, ~b) * Sum(~y, ~i, ~a, ~b), nothing)
    @rule Sum(~x, ~i, ~a, ~b) => ifelse(~i ∈ atoms(~x), nothing, ((~b) - (~a) + 1) * (~x))
]

export rewrite

function rewrite(s)
    # TODO: Put this back once this issue is fixed:
    # https://github.com/JuliaSymbolics/SymbolicUtils.jl/issues/175

    s = symify(s)
    # s=  simplify(s ; polynorm=true)
    r1 = RW.Postwalk(RW.Fixpoint(RW.Chain(POSTRULES)))
    r2 = RW.Prewalk(RW.Fixpoint(RW.Chain(PRERULES)))
    s = RW.Chain([r1, r2])(s)

    s = simplify(s)
end


text:
# geometry setup in domain [0,1]^2
solid_left = 0.45
solid_right = 0.55
solid_top = 0.5
solid_bottom = 0.52

chip_left = 0.48
chip_right = 0.52
chip_top = 0.5
chip_bottom = 0.505

delta=1e-5
delta2=1e-5

filename = "mesh/CHT_2D.stl"
file_format = "stl"
mesh = Mesh(filename, file_format = file_format, degree=2)
mesh = Mesh(mesh.nodes ./ 0.030500000342726707, mesh.elems, -1, 2)
h = 0.1

# m = 100
# n = 100
# h = 1/n
# mesh = Mesh(m, n, h, degree=2)

nnode = mesh.nnode
nedge = mesh.nedge
ndof = mesh.ndof
nelem = mesh.nelem
ngauss = get_ngauss(mesh)

# compute solid indices and chip indices
solid_fem_idx = Array{Int64, 1}([])
solid_fvm_idx = Array{Int64, 1}([])
chip_fem_idx = Array{Int64, 1}([])
# chip_fvm_idx = Array{Int64, 1}([])
chip_fem_top_idx = Array{Int64, 1}([])
bd = Array{Int64, 1}([])
fvm_bd = Array{Int64, 1}([])

for j = 1:nnode
    nodex, nodey = mesh.nodes[j, 1], mesh.nodes[j, 2]
    if nodex >= solid_left-delta2 && nodex <= solid_right+delta2 && nodey >= solid_top-delta2 && nodey <= solid_bottom+delta2
        # print(i, j)
        global solid_fem_idx = [solid_fem_idx; j]
        if nodex >= chip_left-delta2 && nodex <= chip_right+delta2 && nodey >= chip_top-delta2 && nodey <= chip_bottom+delta2
            global chip_fem_idx = [chip_fem_idx; j]
            if nodey <= chip_top+delta2
                global chip_fem_top_idx = [chip_fem_top_idx; j]
            end
        end
    end
    if abs(nodex-0.0) <= delta || abs(nodex-1.0) <= delta || abs(nodey-0.0) <= delta || abs(nodey-1.0) <= delta
        global bd = [bd; j]
    end
end

# save a copy of solid_fem_idx without edges for plotting
chip_fem_idx_nodes = copy(chip_fem_idx)

# fix chip_fem_top_idx
if size(chip_fem_top_idx, 1) == 0
    chip_fem_top_idx = chip_fem_idx
end

for j = 1:nedge
    edgex, edgey = (mesh.nodes[mesh.edges[j, 1], :] .+ mesh.nodes[mesh.edges[j, 2], :]) ./ 2
    if edgex >= solid_left-delta2 && edgex <= solid_right+delta2 && edgey >= solid_top-delta2 && edgey <= solid_bottom+delta2
        # print(i, j)
        global solid_fem_idx = [solid_fem_idx; nnode + j]
        if edgex >= chip_left-delta2 && edgex <= chip_right+delta2 && edgey >= chip_top-delta2 && edgey <= chip_bottom+delta2
            global chip_fem_idx = [chip_fem_idx; nnode + j]
            if edgex >= chip_left-delta2 && edgex <= chip_right+delta2 && edgey >= chip_top-delta2 && edgey <= chip_top+delta2
                global chip_fem_top_idx = [chip_fem_top_idx; nnode + j]
            end
        end
    end
    if abs(edgex-0.0) <= delta || abs(edgex-1.0) <= delta || abs(edgey-0.0) <= delta || abs(edgey-1.0) <= delta
        global bd = [bd; nnode + j]
    end
end

gaussxy = gauss_nodes(mesh)

for i in 1:nelem
    gaussx, gaussy = mean(gaussxy[6*i-5: 6*i, 1]), mean(gaussxy[6*i-5: 6*i, 2])
    if gaussx >= solid_left-delta2 && gaussx <= solid_right+delta2 && gaussy >= solid_top-delta2 && gaussy <= solid_bottom+delta2
        global solid_fvm_idx = [solid_fvm_idx; i]
    end
    if abs(gaussy - 0.0) <= h/2 + delta && ( abs(gaussx - 0.0) <= h + delta || abs(gaussx - 1.0) <= h + delta )
        global fvm_bd = [fvm_bd; i]
    end
end


text:
using sb, LinearAlgebra, StaticArrays

"""

Structure to hold data for ECEID dynamics

"""
mutable struct ECEIDOps
    No :: Int64
    ρ  :: SArray{Tuple{2,2},Complex{Float64},2,4}
    Cs :: Vector{SArray{Tuple{2,2},Complex{Float64},2,4}}
    Cc :: Vector{SArray{Tuple{2,2},Complex{Float64},2,4}}
    As :: Vector{SArray{Tuple{2,2},Complex{Float64},2,4}}
    Ac :: Vector{SArray{Tuple{2,2},Complex{Float64},2,4}}
    N  :: Vector{Float64}
    p  :: Vector{Float64}
    q  :: Vector{Float64}
    function ECEIDOps(ρ0, sbm :: SBModel)
        Cs = Vector{SArray{Tuple{2,2},Complex{Float64},2,4}}(undef,sbm.No)
        Cc = Vector{SArray{Tuple{2,2},Complex{Float64},2,4}}(undef,sbm.No)
        As = Vector{SArray{Tuple{2,2},Complex{Float64},2,4}}(undef,sbm.No)
        Ac = Vector{SArray{Tuple{2,2},Complex{Float64},2,4}}(undef,sbm.No)
        N = sb.thermalns(sbm)
        p = zeros(Float64,sbm.No)
        q = zeros(Float64,sbm.No)
        for i in 1:sbm.No
            Cs[i] = SA[0.0im 0.0im; 0.0im 0.0im]
            Cc[i] = SA[0.0im 0.0im; 0.0im 0.0im]
            Ac[i] = SA[0.0im 0.0im; 0.0im 0.0im]
            As[i] = SA[0.0im 0.0im; 0.0im 0.0im]
        end
        new(sbm.No,ρ0,Cs,Cc,As,Ac,N,p,q)
    end

    return nothing
end

"""

Leapfrog bootstrap of data (ECEID)

"""
function ECEIDbootstrap!(ops :: ECEIDOps, oldops :: ECEIDOps,
                      dotops :: ECEIDOps, dt :: Float64)
    oldops.ρ = ops.ρ - dt * dotops.ρ
    Threads.@threads for i in 1:ops.No
        oldops.Cs[i] = ops.Cs[i] - dt * dotops.Cs[i]
        oldops.Cc[i] = ops.Cc[i] - dt * dotops.Cc[i]
        oldops.As[i] = ops.As[i] - dt * dotops.As[i]
        oldops.Ac[i] = ops.Ac[i] - dt * dotops.Ac[i]
        oldops.N[i]  = ops.N[i] - dt * dotops.N[i]
    end

    return nothing
end

"""

Leapfrog step forward (ECEID)

"""
function ECEIDforward!(ops :: ECEIDOps, oldops :: ECEIDOps,
                    dotops :: ECEIDOps, dt :: Float64)
    oldops.ρ = oldops.ρ + 2 * dt * dotops.ρ
    Threads.@threads for i in 1:ops.No
        oldops.Cs[i] = oldops.Cs[i] + 2 * dt * dotops.Cs[i]
        oldops.Cc[i] = oldops.Cc[i] + 2 * dt * dotops.Cc[i]
        oldops.As[i] = oldops.As[i] + 2 * dt * dotops.As[i]
        oldops.Ac[i] = oldops.Ac[i] + 2 * dt * dotops.Ac[i]
        oldops.N[i] = oldops.N[i] + 2 * dt * dotops.N[i]
    end

    return nothing
end

"""

Leapfrog ECEID calculate EOM RHS

"""
function ECEIDcalcdots!(dotops :: ECEIDOps, ops :: ECEIDOps,
                    dt :: Float64, sbm :: SBModel; thermostat = false)

    dotops.ρ = - IOHBAR * comm(H(sbm,ops.q,ops.p),ops.ρ)
    for i in 1:ops.No
        μ = (1.0im*ops.Cc[i] - ops.As[i]) / sbm.ωs[i]
        dotops.ρ += IOHBAR * comm(F(sbm,ops.q,i),μ)
    end
    Threads.@threads for i in 1:ops.No
        dotops.Cc[i] = - IOHBAR * comm(H(sbm,ops.q,ops.p),ops.Cc[i]) + sbm.ωs[i] * ops.Cs[i] + (ops.N[i] + 0.5) * comm(F(sbm,ops.q,i),ops.ρ)
        dotops.Cs[i] = - IOHBAR * comm(H(sbm,ops.q,ops.p),ops.Cs[i]) - sbm.ωs[i] * ops.Cc[i]
        dotops.Ac[i] = - IOHBAR * comm(H(sbm,ops.q,ops.p),ops.Ac[i]) + sbm.ωs[i] * ops.As[i] + 0.5 * acomm(F(sbm,ops.q,i),ops.ρ)
        dotops.As[i] = - IOHBAR * comm(H(sbm,ops.q,ops.p),ops.As[i]) - sbm.ωs[i] * ops.Ac[i]
        λ = (1.0im*ops.Cs[i] + ops.Ac[i])
        if !thermostat
            dotops.N[i] = real(tr(F(sbm,ops.q,i) * λ)) / HBAR / sbm.ωs[i]
        else
            dotops.N[i] = 0.0
        end
    end

    return nothing
end

"""

ECEID total energy

"""
function ECEIDenergy(ops :: ECEIDOps, sbm :: SBModel)
    eenergy = real(tr(H(sbm,ops.q,ops.p) * ops.ρ))
    oenergy = 0.0
    intenergy = 0.0
    for i in 1:ops.No
        oenergy += HBAR * sbm.ωs[i] * (ops.N[i] + 0.5)
        μ = (im*ops.Cc[i]-ops.As[i])/sbm.ωs[i]
        intenergy -= real(tr(F(sbm,ops.q,i) * μ))
    end
    return eenergy+oenergy+intenergy
end

"""

This function integrates the ECEID EOMs starting from the
initial condition provided in ops for nsteps of dt.
The store! callback function is called at every timestep.

"""
function RunECEID!(sbm,ops,nsteps,dt,store!,storage; thermostat = false)

    dotops = ECEIDOps(SA[0.0im 0.0im; 0.0im 0.0im],sbm)
    oldops = ECEIDOps(SA[0.0im 0.0im; 0.0im 0.0im],sbm)

    # Bootstrap the integrator
    sb.ECEIDcalcdots!(dotops, ops , dt , sbm; thermostat = thermostat)
    sb.ECEIDbootstrap!(ops,oldops,dotops,dt)

    # call results storage function at t = 0
    t = 0.0
    store!(storage,t,ops,sbm)

    # integrate
    for i in 1:nsteps-1
        t = i*dt
        sb.ECEIDcalcdots!(dotops,ops,dt,sbm; thermostat = thermostat)
        sb.ECEIDforward!(ops,oldops,dotops,dt)
        (ops,oldops) = (oldops,ops)
        store!(storage,t,ops,sbm)
    end

    return nothing
end