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

using Markdown
using InteractiveUtils

# ╔═╡ ffed97e2-dc07-48bb-870a-737a115e69dc
begin
	using Pkg
	cd(joinpath(dirname(@__FILE__),".."))
    Pkg.activate(pwd())
    using NativeSVG
	using Plots
	using LaTeXStrings
	using GeneralQP
	using NLopt
end

# ╔═╡ 9f5f8730-0188-11eb-209f-e771579d35ce
md"# Sequential Quadratic Programming"

# ╔═╡ c30dbe90-0188-11eb-0c2b-ed5e93a86b0b
md"""## Quadratic Programming

An optimization problem with a quadratic objective function and linear
constraints is called a _quadratic program_. Problems of this
type are important in their own right, and they also arise a subproblems
in methods for general constrained optimization such as sequential
quadratic programming (this chapter) and interior-point methods (lecture 10).

The general quadratic program can be stated as

```math
\begin{aligned}
\min_{\vec{x}}\, & f\left(\vec{x}\right)\overset{\vartriangle}{=}\frac{1}{2}\vec{x}^\mathsf{T}Q\vec{x}-\vec{c}^\mathsf{T}\vec{x}\\
\textrm{subject to}\, & \begin{cases}
A_{\textrm{eq}}\vec{x}=\vec{b}_{\textrm{eq}}\,,\\
A_{\textrm{in}}\vec{x}\leq\vec{b}_{\textrm{in}}\,,
\end{cases}
\end{aligned}
```

where $Q$ is a symmetric $n\times n$ matrix, $\vec{c}\in\mathbb R^{n}$,
$A_{\textrm{eq}}$ is a $m\times n$ matrix, $\vec{b}_{\textrm{eq}}\in\mathbb R^{m}$,
$A_{\textrm{in}}$ is a $p\times n$ matrix and $\vec{b}_{\textrm{in}}\in\mathbb R^{p}$.
If the Hessian matrix $Q$ is positive semidefinite, we have a convex
quadratic program. Non-convex quadratic programs, in which $Q$ is
an indefinite matrix, can be more challenging because they have several
stationary points and local minima."""

# ╔═╡ e0be4cbe-0188-11eb-1991-5b001cf2434a
md"""## Equality Constraints

We begin our discussion of algorithms for quadratic programming by considering
the case in which only equality constraints are present. We consider
to following equality-constrained quadratic problem
```math
\begin{aligned}
\min_{\vec{x}}\, & f\left(\vec{x}\right)\overset{\vartriangle}{=}\frac{1}{2} \vec{x}^\mathsf{T}Q\vec{x}- \vec{c}^\mathsf{T}\vec{x}\\
\textrm{subject to}\, & A\vec{x}=\vec{b}\,,\nonumber
\end{aligned}
```
where $A$ is the $m\times n$ Jacobian of constraints and $\vec{b}\in\mathbb R^{m}$.
We assume that $A$ has rank $m$ so that the constraints are consistent.

The First Order Necessary Condition (FONC) uses the Langrangian function
```math
\mathcal L\left(\vec x, \vec \lambda\right)=\vec{x}^\mathsf{T}Q\vec{x}- \vec{c}^\mathsf{T}\vec{x} + \vec\lambda^\mathsf{T} \left(A\vec{x}-\vec{b}\right)
```
and gives the following condition
```math
\begin{aligned}
\nabla_{\vec x}\mathcal L\left(\vec x ^\star, \vec \lambda ^\star\right)&=\vec 0\,,\\
\nabla_{\vec \lambda}\mathcal L\left(\vec x ^\star, \vec \lambda ^\star\right)&=\vec 0\,,
\end{aligned}
```
So the FONC for $\vec{x}^{\star}$ to be a solution of the quadratic
problem yields a vector $\vec{\lambda}^{\star}$ such
that the following system of equations is satisfied:
```math
\begin{pmatrix}Q & A^\mathsf{T}\\
A & 0
\end{pmatrix}\begin{pmatrix}\vec{x}^{\star}\\
\vec{\lambda}^{\star}
\end{pmatrix}=\begin{pmatrix}\vec{c}\\
\vec{b}
\end{pmatrix}\,.
```
This system can be solved directly by factorization. An alternative
is to use an iterative method.

Consider the quadratic programming problem
```math
\begin{aligned}
\min_{\vec{x}}\, & \frac{1}{2} \vec{x}^\mathsf{T}\begin{pmatrix}6 & 2 & 1\\
2 & 5 & 2\\
1 & 2 & 4
\end{pmatrix}\vec{x}- \begin{pmatrix}8\\
3\\
3
\end{pmatrix}^\mathsf{T}\vec{x}\\
\textrm{subject to}\, & \begin{pmatrix}1 & 0 & 1\\
0 & 1 & 1
\end{pmatrix}\vec{x}=\begin{pmatrix}3\\
0
\end{pmatrix}\,.
\end{aligned}
```"""

# ╔═╡ 56fd4bc0-0189-11eb-04c2-0d000ded6e97
let
	Q = [6 2 1
     	 2 5 2
         1 2 4]
	c = [8, 3, 3]
	A = [1 0 1
		 0 1 1]
	b = [3, 0]
	sol = [Q transpose(A)
		   A zeros(2,2)] \ [c; b]
end

# ╔═╡ 80edfd80-0189-11eb-1687-01a14b31ad39
md"""The solution $\vec{x}^{\star}$ and optimal Lagrange multiplier vector $\vec \lambda^\star$ are given by
```math
\begin{aligned}
\vec{x}^{\star} & = \begin{pmatrix}2 & -1 & 1\end{pmatrix}^\mathsf{T}\,,\\
\vec{\lambda}^{\star} & = \begin{pmatrix}-3 & 2\end{pmatrix}^\mathsf{T}\,.
\end{aligned}
```"""

# ╔═╡ 8f340330-0189-11eb-27c8-7ba7c0954ac7
md"""## Active Set Method

We now describe active-set methods for solving quadratic programs
containing both equality and inequality constraints. We consider only
the convex case, in which the matrix ``Q`` is positive semidefinite.

If the contents of the optimal active set ``J\left(\vec x^\star\right)`` were known in advance, we could find the solution ``\vec{x}^{\star}``
by applying the technique for equality constrained quadratic programs
to the problem
```math
\begin{aligned}
\min_{\vec{x}}\, & f\left(\vec{x}\right)\overset{\vartriangle}{=}\frac{1}{2} \vec{x}^\mathsf{T}Q\vec{x}- \vec{c}^\mathsf{T}\vec{x}\\
\textrm{subject to}\, & \begin{cases}
A_{\textrm{eq}}\vec{x}=\vec{b}_{\textrm{eq}}\,,\\
 \vec{a}_{\textrm{in},j}^\mathsf{T}\vec{x}=b_{\textrm{in},j}\,, & \forall j\in J\left(\vec{x}^{\star}\right)\,,
\end{cases}
\end{aligned}
```
where ``\vec{a}_{\textrm{in},j}`` is the ``j``th row in the matrix ``A_{\textrm{in}}``
and ``b_{\textrm{in},j}`` is the ``j``th element of the vector ``\vec{b}_{\textrm{in}}``.
Of course, we usually do not have prior knowledge of ``J\left(\vec{x}^{\star}\right)``
and determination of this set is the main challenge facing algorithms
for inequality-constrained quadratic programs.

Active-set methods find a step from one iterate to the next by solving
a quadratic subproblem in which some of the inequality constraints,
and all the equality constraints are imposed as equalities. This subset
is referred to as the _working set_ and is denoted at the ``k``th
iterate by ``W_{k}``. An important requirement we impose on ``W_{k}``
is that the gradients ``\vec{a}_{\textrm{eq},i}``, ``i=1,\dots,m`` and
``\vec{a}_{\textrm{in},j}``, ``j\in W_{k}`` are linearly independent,
even when the full set of active constraints at that point has linearly
dependent gradients.

Given an iterate ``\vec{x}^{\left(k\right)}`` and the working set ``W_{k}``,
we first check whether ``\vec{x}^{\left(k\right)}`` minimizes the quadratic
function ``f`` in the subspace defined by the working set. If not,
we compute a step ``\vec{d}^{\left(k\right)}`` by solving an equality-constrained
quadratic subproblem in which the inequality constraints corresponding
to the working set ``W_{k}`` are regarded as equalities and the other
inequality constraints are temporarily disregarded. To express this
subproblem in terms of the step ``\vec{d}^{\left(k\right)}``, we define
```math
\vec{d}^{\left(k\right)}=\vec{x}^{\left(k+1\right)}-\vec{x}^{\left(k\right)}\,,\quad\vec{g}_{k}=Q\vec{x}^{\left(k\right)}-\vec{c}\,.
```
By substituting for ``\vec{x}^{\left(k+1\right)}`` into the objective
function, we find that
```math
f\left(\vec{x}^{\left(k+1\right)}\right)=f\left(\vec{x}^{\left(k\right)}+\vec{d}^{\left(k\right)}\right)=\frac{1}{2}\left(\vec{d}^{\left(k\right)}\right)^\mathsf{T}Q\vec{d}^{\left(k\right)}+\vec{g}_{k}^\mathsf{T}\vec{d}^{\left(k\right)}+\rho_{k}\,,
```
where ``\rho_{k}=\frac{1}{2} \left(\vec{x}^{\left(k\right)}\right)^\mathsf{T}Q\vec{x}^{\left(k\right)}- \vec{c}^\mathsf{T}\vec{x}^{\left(k\right)}``
is independent of ``\vec{d}^{\left(k\right)}``. Since we can drop ``\rho_{k}``
from the objective function without changing the solution of the problem,
we can write the subproblem to be solved at the ``k``th iteration as
follows
```math
\begin{aligned}
\min_{\vec{d}^{\left(k)\right)}}\, & \frac{1}{2} \left(\vec{d}^{\left(k\right)}\right)^\mathsf{T}Q\vec{d}^{\left(k\right)}+ \vec{g}_{k}^\mathsf{T}\vec{d}^{\left(k\right)}\label{eq:quadratic_subproblem}\\
\textrm{subject to}\, & \begin{cases}
A_{\textrm{eq}}\vec{d}^{\left(k\right)}=\vec{0}\,,\\
 \vec{a}_{\textrm{in},j}^\mathsf{T}\vec{d}^{\left(k\right)}=0\,, & \forall j\in W_{k}\,.
\end{cases}
\end{aligned}
```
Note that for each ``j\in W_{k}``, the value of `` \vec{a}_{\textrm{in},j}^\mathsf{T}\vec{x}^{\left(k\right)}``
does not change as we move along ``\vec{d}^{\left(k\right)}``, since
we have `` \vec{a}_{\textrm{in},j}^\mathsf{T}\left(\vec{x}^{\left(k\right)}+\alpha\vec{d}^{\left(k\right)}\right)= \vec{a}_{\textrm{in},j}^\mathsf{T}\vec{x}^{\left(k\right)}=b_{\textrm{in},j}``
for all ``\alpha``. Since the constraints in ``W_{k}`` were satisfied
at ``\vec{x}^{\left(k\right)}``, they are also satisfied at ``\vec{x}^{\left(k\right)}+\alpha\vec{d}^{\left(k\right)}``,
for any value of ``\alpha``.

Supposing for the moment that the optimal ``\vec{d}^{\left(k\right)}``
is nonzero, we need to decide how far to move along this direction.
If ``\vec{x}^{\left(k\right)}+\vec{d}^{\left(k\right)}`` is feasible
with respect to all the constraints, we set ``\vec{x}^{\left(k+1\right)}=\vec{x}^{\left(k\right)}+\vec{d}^{\left(k\right)}``.
Otherwise, we set
```math
\vec{x}^{\left(k+1\right)}=\vec{x}^{\left(k\right)}+\alpha_{k}\vec{d}^{\left(k\right)}\,,
```
where the step-length parameter ``\alpha_{k}`` is chosen to be the
largest value in the range ``\left[0,1\right]`` for which all constraints
are satisfied. We can derive an explicit definition of ``\alpha_{k}``
by considering what happens to the constraints ``j\notin W_{k}``, since
the constraints ``j\in W_{k}`` will certainly be satisfied regardless
of the choice of ``\alpha_{k}``. If `` \vec{a}_{\textrm{in},j}^\mathsf{T}\vec{d}^{\left(k\right)}\leq0``
for some ``j\notin W_{k}``, then for all ``\alpha_{k}\geq0``, we have
`` \vec{a}_{\textrm{in},j}^\mathsf{T}\left(\vec{x}^{\left(k\right)}+\alpha_{k}\vec{d}^{\left(k\right)}\right)\leq \vec{a}_{\textrm{in},j}^\mathsf{T}\vec{x}^{\left(k\right)}\leq b_{\textrm{in},j}``.
Hence, constraint ``j`` will be satisfied for all nonnegative choices
of the step-length parameter. Whenever `` \vec{a}_{\textrm{in},j}^\mathsf{T}\vec{d}^{\left(k\right)}>0``
for some ``j\notin W_{k}``, however, we have that `` \vec{a}_{\textrm{in},j}^\mathsf{T}\left(\vec{x}^{\left(k\right)}+\alpha_{k}\vec{d}^{\left(k\right)}\right)\leq b_{\textrm{in},j}``
only if
```math
\alpha_{k}\leq\frac{b_{\textrm{in},j}- \vec{a}_{\textrm{in},j}^\mathsf{T}\vec{x}^{\left(k\right)}}{ \vec{a}_{\textrm{in},j}^\mathsf{T}\vec{d}^{\left(k\right)}}\,.
```
To maximize the decrease in ``f``, we want ``\alpha_{k}`` to be as large
as possible in ``\left[0,1\right]`` subject to retaining feasibility,
so we obtain the following definition
```math
\alpha_{k}\overset{\textrm{def}}{=}\min\left\{ 1,\min_{j\notin W_{k}, \vec{a}_{\textrm{in},j}^\mathsf{T}\vec{d}^{\left(k\right)}>0}\frac{b_{\textrm{in},j}- \vec{a}_{\textrm{in},j}^\mathsf{T}\vec{x}^{\left(k\right)}}{ \vec{a}_{\textrm{in},j}^\mathsf{T}\vec{d}^{\left(k\right)}}\right\} \,.
```
We call the constraints ``j`` for which this minimum is achieved the
_blocking constraints_. Note that it is quite possible for ``\alpha_{k}``
to be zero, because we could have `` \vec{a}_{\textrm{in},j}^\mathsf{T}\vec{d}^{\left(k\right)}>0``
for some constraint ``j`` that is active at ``\vec{x}^{\left(k\right)}``
but not a member of the current working set ``W_{k}``.

If ``\alpha_{k}<1``, that is, the step along ``\vec{d}_{k}`` was blocked
by some constraint not in ``W_{k}``, a new working set ``W_{k+1}`` is
constructed by adding one of the blocking constraints to ``W_{k}``.
We continue to iterate in this manner, adding constraints to the working
set until the subproblem has solution ``\vec{d}^{\circ}=\vec{0}``.
Since ``\vec{d}^{\circ}=\vec{0}`` satisfy the optimality condition, we have that
```math
\sum_{i=1}^{m}\lambda_{i}^{\circ}\vec{a}_{\textrm{eq},i}+\sum_{j\in W^{\circ}}\mu_{j}^{\circ}\vec{a}_{\textrm{in},j}=-\vec{g}^{\circ}=-Q\vec{x}^{\circ}+\vec{c}\,,
```
for some Lagrange multipliers ``\lambda_{i}^{\circ}``, ``i=1,\dots,m``
and ``\mu_{j}^{\circ}``, ``j\in W^{\circ}``. It follows that ``\vec{x}^{\circ}``,
``\vec{\lambda}^{\circ}`` and ``\vec{\mu}^{\circ}`` satisfy the second
KKT condition, if we define the multipliers corresponding to the inequality
constraints not in the working set to be zero. Because of the control
imposed on the step length, ``\vec{x}^{\circ}`` is also feasible with
respect to all the constraints, so the third, fourth and fifth KKT
conditions are satisfied at this point.

We now examine the signs of the KKT multipliers in the working set,
that is, the indices ``j\in W^{\circ}``. If these conditions are all
nonnegative, the first KKT condition is also satisfied, so we conclude
that ``\vec{x}^{\circ}`` is a KKT point for the original problem. In
fact, since ``Q`` is positive semidefinite, we have that ``\vec{x}^{\circ}``
is a global minimum.

If, on the other hand, on or more of the multipliers ``\mu_{j}^{\circ}``,
``j\in W^{\circ}``, are negative, the first KKT condition is not satisfied
and the objective function ``f`` may be decreased by dropping one of
these constraints. Thus, we remove an index ``j`` corresponding to
one of the negative multipliers from the working set and solve a new
subproblem for the next step. While any index ``j`` for which ``\mu_{j}^{\circ}<0``
usually will yield in a direction ``\vec{d}`` along which the algorithm
can make progress, the most negative multiplier is often chosen in
practice. This choice is motived by a sensitivity analysis, which
shows that the rate of decrease in the objective function when one
constraint is removed, is proportional to the magnitude of the Lagrange
multiplier for that constraint."""

# ╔═╡ db088c90-0189-11eb-338d-81af9469d556
md"""## Example

Apply the active-set method to the following problem:
```math
\begin{aligned}
\min_{\vec{x}}\, & f\left(\vec{x}\right)=\left(x_{1}-1\right)^{2}+\left(x_{2}-2.5\right)^{2}\\
\textrm{subject to}\, & \begin{cases}
-x_{1}+2x_{2}-2\leq0\,,\\
x_{1}+2x_{2}-6\leq0\,.\\
x_{1}-2x_{2}-2\leq0\,,\\
-x_{1}\leq0\,,\\
-x_{2}\leq0\,.
\end{cases}
\end{aligned}
```"""

# ╔═╡ 6bdfb630-018a-11eb-1010-b3c332121581
begin
	Q = Float64[2 0
                0 2]
	c = Float64[2, 5]
	A = Float64[-1  2
				 1  2
				 1 -2
				-1  0
				 0 -1]
	b = Float64[2, 6, 2, 0, 0]

	x = -1:0.05:5
	y = -1:0.05:5
	z = Surface((x,y)->(0.5 .* [x y]*Q*[x;y] .- transpose(c)*[x;y] .+ 0.5 .* transpose(c)*Q*c)[1], x, y)
	contour(x, y, z, levels=35)
	plot!(x, (2 .+ x) ./ 2, linestyle=:dash, label=L"-x_1+2x_2-2\le 0")
	plot!(x, (6 .- x) ./ 2, linestyle=:dash, label=L"x_1+2x_2-6\le 0")
	plot!(x, (2 .- x) ./ -2, linestyle=:dash, label=L"x_1-2x_2-2\le 0")
	plot!([0,2,4,2,0,0],[0,0,1,2,1,0], linewidth=2, label="domain")
end

# ╔═╡ 88c16630-018b-11eb-3ea1-a555d69b2e26
md"""We refer the constraints, in order, by indices $1$ through $5$.
For this problem it is easy to determine a feasible initial point;
say $\vec{x}^{\left(0\right)}=\begin{pmatrix}2 & 0\end{pmatrix}\mathsf{T}$.
Constraints $3$ and $5$ are active at this point, and we set $W_{0}=\left\{ 3,5\right\} $.
Note that we could just as validly have chosen $W_{0}=\left\{ 5\right\} $
or $W_{0}=\left\{ 3\right\} $ or even $W_{0}=\emptyset$; each choice
would lead the algorithm to perform somewhat differently."""

# ╔═╡ 4a460a20-018e-11eb-3da1-4d37893ead18
x₀ = [2, 0]

# ╔═╡ 8df4e190-018b-11eb-2eae-1d8f6c4e6248
let x = x₀
	g = Q*x - c
	ain = [reshape(A[3,:], 1, 2); reshape(A[5,:], 1, 2)]
	sol = [Q transpose(ain)
		   ain zeros(2,2)] \ [-g; 0; 0]
end

# ╔═╡ a8042d20-018b-11eb-0f74-3b70e587f20d
md"""Since $\vec{x}^{\left(0\right)}$ lies on a vertex of the feasible
region, it is obviously a minimizer of the objective function $f$
with respect to the working set $W_{0}$; that is, the solution of
the subproblem with $k=0$ is $\vec{d}^{\left(0\right)}=\vec{0}$.
We can then find the multipliers $\mu_{3}^{\circ}$ and $\mu_{5}^{\circ}$
associated with the active constraints. Substitution of the data from
our problem yields
```math
\begin{pmatrix}-1\\
2
\end{pmatrix}\lambda_{3}^{\circ}+\begin{pmatrix}0\\
1
\end{pmatrix}\lambda_{5}^{\circ}=\begin{pmatrix}2\\
-5
\end{pmatrix}\,,
```
which has solution $\lambda_{3}^{\circ}=-2$ and $\lambda_{5}^{\circ}=-1$.

We now remove constraint $3$ from the working set, because it has
the most negative multiplier, and set $W_{1}=\{5\}$."""

# ╔═╡ b9b7e250-018b-11eb-1eca-51dc8250a3eb
sol₁ = let x = x₀
	g = Q*x - c
	ain = reshape(A[5,:], 1, 2)
	sol = [Q transpose(ain)
		   ain zeros(1,1)] \ [-g; 0]
end

# ╔═╡ e14bb850-018b-11eb-0bb0-15ebdb2f7fcf
md"""We begin iteration $1$ by finding the solution of the subproblem for $k=1$, which is
$\vec{d}^{\left(1\right)}= \begin{pmatrix}-1 & 0\end{pmatrix}^\mathsf{T}$."""

# ╔═╡ ea6595a0-018b-11eb-07dc-f394137f59f4
α₁, d₁ = let x = x₀
	d = sol₁[1:2]
	α = min(1.0, [(reshape(A[j,:], 1, 2) * d)[1,1] ≤ 0 ? 1.0 : ((b[j] .- reshape(A[j,:], 1, 2) * x) / (reshape(A[j,:], 1, 2) * d))[1,1] for j in (1,2,3,4)]...)
	α, d
end

# ╔═╡ 0d2eeb40-018c-11eb-0fdc-1b6668c5941f
md"""The step-length formula yields $\alpha_{1}=1$, and the new iterate
is $\vec{x}^{\left(2\right)}=\begin{pmatrix}1 & 0\end{pmatrix}^\mathsf{T}$. There are no blocking constraints, so that ``W_{2}=W_{1}=\left\{ 5\right\}``"""

# ╔═╡ ee702bf0-018c-11eb-1ec9-6b481c119add
x₁ = x₀ + α₁ .* d₁

# ╔═╡ 2cb44400-018d-11eb-0138-c78f6a18b8f8
sol₂ = let x = x₁
	g = Q*x - c
	ain = reshape(A[5,:], 1, 2)
	sol = [Q transpose(ain)
		   ain zeros(1,1)] \ [-g; 0]
end

# ╔═╡ 456d9c82-018d-11eb-0e81-bf7346539a32
md"""We find at the start of iteration $2$ that the solution of the
subproblem is $\vec{d}^{\left(2\right)}=\vec{0}$. We deduce that
the Lagrange multiplier for the lone working constraint is $\lambda_{5}^{\circ}=-5$,
so we drop the working set to obtain $W_{3}=\emptyset$."""

# ╔═╡ 84712210-018f-11eb-0bef-53cdfb93e0e2
x₂ = x₁

# ╔═╡ 96ecfeb2-018e-11eb-2dda-994575ee6011
sol₃ = let x = x₂
	g = Q*x - c
	sol = Q \ (-g)
end

# ╔═╡ 2ddf0860-018d-11eb-2884-7dacdda18597
md"""Iteration $3$ starts by solving the unconstrained problem, to obtain
the solution $\vec{d}^{\left(3\right)}=\begin{pmatrix}0 & 2.5\end{pmatrix}^\mathsf{T}$."""

# ╔═╡ c71e2f00-018e-11eb-1c30-a54d6d4d04ff
α₃, d₃ = let x = x₂
	d = sol₃[1:2]
	α = min(1.0, [(reshape(A[j,:], 1, 2) * d)[1,1] ≤ 0 ? 1.0 : ((b[j] .- reshape(A[j,:], 1, 2) * x) / (reshape(A[j,:], 1, 2) * d))[1,1] for j in (1,2,3,4,5)]...)
	α, d
end

# ╔═╡ e470e5c0-018e-11eb-229b-e521909f84e3
md"""The step-length formula yields a step length of $\alpha_{3}=0.6$
and a new iterate $\vec{x}^{\left(4\right)}=\begin{pmatrix}1 & 1.5\end{pmatrix}^\mathsf{T}$. There
is a single blocking constraint (constraint $1$), so we obtain ``W_{4}=\left\{ 1\right\}``."""

# ╔═╡ 596ba4f0-018f-11eb-321e-03ec8a6d18db
x₃ = x₂ + α₃ .* d₃

# ╔═╡ 91ada700-018f-11eb-37d0-757a263bae99
sol₄ = let x=x₃
	g = Q*x - c
	ain = reshape(A[1,:], 1, 2)
	sol = [Q transpose(ain)
		   ain zeros(1,1)] \ [-g; 0]
end

# ╔═╡ abfbf760-018f-11eb-18dc-bd573cc8be06
md"""The solution of the subproblem for $k=4$ is then $\vec{d}^{\left(4\right)}=\begin{pmatrix}0.4 & 0.2\end{pmatrix}^\mathsf{T}$."""

# ╔═╡ ba0f5450-018f-11eb-2188-4327aa9c02cc
α₄, d₄ = let x = x₃
	d = sol₄[1:2]
	α = min(1.0, [(reshape(A[j,:], 1, 2) * d)[1,1] ≤ 0 ? 1.0 : ((b[j] .- reshape(A[j,:], 1, 2) * x) / (reshape(A[j,:], 1, 2) * d))[1,1] for j in (2,3,4,5)]...)
	α, d
end

# ╔═╡ d8c62ef0-018f-11eb-21ce-ebe037f9d91d
md"""The new step-length is $1$. There are no blocking constraints
on this step, so the next working set in unchanged: $W_{5}=\left\{ 1\right\} $. The new iterate is $\vec{x}^{\left(5\right)}=\begin{pmatrix}1.4 & 1.7\end{pmatrix}^\mathsf{T}$."""

# ╔═╡ e59a7c30-018f-11eb-2b94-7b8dff2a6383
x₄ = x₃ + α₄ .* d₄

# ╔═╡ f5f705d0-018f-11eb-2462-e7e2265b700c
sol₅ = let x = x₄
	g = Q*x - c
	ain = reshape(A[1,:], 1, 2)
	sol = [Q transpose(ain)
		   ain zeros(1,1)] \ [-g; 0]
end

# ╔═╡ 1254d540-0190-11eb-3066-b5a853eef0c0
md"""Finally, we solve the subproblem for $k=5$ to obtain a solution $\vec{d}^{\left(5\right)}=\vec{0}$.
We find a multiplier $\mu_{1}^{\circ}=0.8$, so we have found the
solution. Wet set $\vec{x}^{\star}=\begin{pmatrix}1.4 & 1.7\end{pmatrix}^\mathsf{T}$
and terminate."""

# ╔═╡ 02787820-0190-11eb-24ed-c5296c91df6b
let
	x = -1:0.05:5
	y = -1:0.05:5
	z = Surface((x,y)->(0.5 .* [x y]*Q*[x;y] .- transpose(c)*[x;y] .+ 0.5 .* transpose(c)*Q*c)[1], x, y)
	contour(x, y, z, levels=35)
	plot!(x, (2 .+ x) ./ 2, linestyle=:dash, label=L"-x_1+2x_2-2\le 0")
	plot!(x, (6 .- x) ./ 2, linestyle=:dash, label=L"x_1+2x_2-6\le 0")
	plot!(x, (2 .- x) ./ -2, linestyle=:dash, label=L"x_1-2x_2-2\le 0")
	plot!([0,2,4,2,0,0],[0,0,1,2,1,0], linewidth=2, label="domain")
	plot!([2,2,1,1,1,1.4], [0,0,0,0,1.5,1.7], linewidth=2, label="iterates", markershape=:circle)
end

# ╔═╡ aa606210-018a-11eb-29bd-5d4f9ae0966e
md"""## Julia

The package `GeneralQP` solves the quadratic programming problems using the function `solve(Q, -c, A, b, x₀)`."""

# ╔═╡ f3007460-018a-11eb-3fa2-8753aceaf89a
sol = GeneralQP.solve(Q, -c, A, b, Float64[2,0])

# ╔═╡ 24c67580-018b-11eb-3eb7-47a2e2ae5251
md"""## Sequential Quadratic Programming

We consider the general constrained problem
```math
\begin{aligned}
\min\, & f\left(\vec{x}\right)\\
\textrm{subject to} & \begin{cases}
\vec{h}\left(\vec{x}\right)=\vec{0}\\
\vec{g}\left(\vec{x}\right)\leq\vec{0}
\end{cases}
\end{aligned}
```
where $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$, $\vec{h}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$,
$m\leq n$, and $\vec{g}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{p}$. The idea
behind the _sequential quadratic programming_ (SQP) approach
is to model the general problem at the current iterate $\vec{x}^{\left(k\right)}$
by a quadratic programming subproblem, then use the minimizer of this
subproblem to define a new iterate $\vec{x}^{\left(k+1\right)}$.
The challenge is to design the quadratic subproblem so that it yields
a good step for the general optimization problem.

We know that the extended Lagrangian function for this problem is

```math
\mathcal{L}\left(\vec{x},\vec{\lambda},\vec{\mu}\right)=f\left(\vec{x}\right)+\vec{\lambda}^\mathsf{T}\vec{h}\left(\vec{x}\right)+\vec{\mu}^\mathsf{T}\vec g\left(\vec{x}\right)\,.
```

Applying Newton's method to the Lagrangian function and linearizing
both the equality and the inequality constraints yields the following
subproblem
```math
\begin{aligned}
\min_{\vec{d}^{\left(k\right)}}\, & \frac{1}{2}\left(\vec{d}^{\left(k\right)}\right)^\mathsf{T} \mathsf{H} \mathcal{L}\left(\vec{x}^{\left(k\right)},\vec{\lambda}^{\left(k\right)},\vec{\mu}^{\left(k\right)}\right)\vec{d}^{\left(k\right)}+ \mathsf{D} \mathcal{L}\left(\vec{x}^{\left(k\right)},\vec{\lambda}^{\left(k\right)},\vec{\mu}^{\left(k\right)}\right)\vec{d}^{\left(k\right)}\\
\textrm{subject to}\, & \begin{cases}
 \mathsf{J}\vec{h}\left(\vec{x}^{\left(k\right)}\right)\vec{d}^{\left(k\right)}=-\vec{h}\left(\vec{x}^{\left(k\right)}\right)\\
 \mathsf{D} g_{j}\left(\vec{x}^{\left(k\right)}\right)\vec{d}^{\left(k\right)}=-g_{j}\left(\vec{x}^{\left(k\right)}\right)\,, & j\in W_{k}\,,
\end{cases}
\end{aligned}
```

where $\mu_{j}^{\left(k\right)}=0$, for all
$j\notin W_{k}$. We can use the active-set method for quadratic programming
to solve this subproblem. The new iterate is given by $\vec{x}^{\left(k+1\right)}$,
$\vec{\lambda}^{\left(k+1\right)}$, $\vec{\mu}^{\left(k+1\right)}$
and $W_{k+1}$.

If the SQP method is able to identify the optimal active set then
it will act like a Newton method for equality-constrained optimization
and will converge rapidly.

It is also remarkable that, far from the solution, the SQP approach
is usually able to improve the estimate of the active set and guide
the iterates towards a solution.

Non-quadratic objective functions, however, can impede progress of
the SQP algorithm, a phenomenon known as the Maratos effect. Steps
that make good progress toward a solution are rejected and the algorithm
fails to converge rapidly. These difficulties can be overcome by means
of a _second-order correction_."""

# ╔═╡ a98c8890-0190-11eb-25a4-bfd18c42b749
function myfunc(x::Vector, grad::Vector)
    if length(grad) > 0
        grad[1] = 0
        grad[2] = 0.5/sqrt(x[2])
    end
    return sqrt(x[2])
end

# ╔═╡ acfea4e0-0190-11eb-07a9-a3f38d97ff80
function myconstraint(x::Vector, grad::Vector, a, b)
    if length(grad) > 0
        grad[1] = 3a * (a*x[1] + b)^2
        grad[2] = -1
    end
    (a*x[1] + b)^3 - x[2]
end

# ╔═╡ a4824010-0190-11eb-0d97-a55353709552
let
	opt = Opt(:LD_SLSQP, 2)
	opt.lower_bounds = [-Inf, 0.]
	opt.xtol_rel = 1e-4

	opt.min_objective = myfunc
	inequality_constraint!(opt, (x,g) -> myconstraint(x,g,2,0), 1e-8)
	inequality_constraint!(opt, (x,g) -> myconstraint(x,g,-1,1), 1e-8)

	(minf, minx, ret) = optimize(opt, [1.234, 5.678])
	numevals = opt.numevals # the number of function evaluations
	minf, minx, numevals, ret
end

# ╔═╡ Cell order:
# ╟─ffed97e2-dc07-48bb-870a-737a115e69dc
# ╟─9f5f8730-0188-11eb-209f-e771579d35ce
# ╟─c30dbe90-0188-11eb-0c2b-ed5e93a86b0b
# ╟─e0be4cbe-0188-11eb-1991-5b001cf2434a
# ╠═56fd4bc0-0189-11eb-04c2-0d000ded6e97
# ╟─80edfd80-0189-11eb-1687-01a14b31ad39
# ╟─8f340330-0189-11eb-27c8-7ba7c0954ac7
# ╟─db088c90-0189-11eb-338d-81af9469d556
# ╠═6bdfb630-018a-11eb-1010-b3c332121581
# ╟─88c16630-018b-11eb-3ea1-a555d69b2e26
# ╠═4a460a20-018e-11eb-3da1-4d37893ead18
# ╠═8df4e190-018b-11eb-2eae-1d8f6c4e6248
# ╟─a8042d20-018b-11eb-0f74-3b70e587f20d
# ╠═b9b7e250-018b-11eb-1eca-51dc8250a3eb
# ╟─e14bb850-018b-11eb-0bb0-15ebdb2f7fcf
# ╠═ea6595a0-018b-11eb-07dc-f394137f59f4
# ╟─0d2eeb40-018c-11eb-0fdc-1b6668c5941f
# ╠═ee702bf0-018c-11eb-1ec9-6b481c119add
# ╠═2cb44400-018d-11eb-0138-c78f6a18b8f8
# ╟─456d9c82-018d-11eb-0e81-bf7346539a32
# ╠═84712210-018f-11eb-0bef-53cdfb93e0e2
# ╠═96ecfeb2-018e-11eb-2dda-994575ee6011
# ╟─2ddf0860-018d-11eb-2884-7dacdda18597
# ╠═c71e2f00-018e-11eb-1c30-a54d6d4d04ff
# ╟─e470e5c0-018e-11eb-229b-e521909f84e3
# ╠═596ba4f0-018f-11eb-321e-03ec8a6d18db
# ╠═91ada700-018f-11eb-37d0-757a263bae99
# ╟─abfbf760-018f-11eb-18dc-bd573cc8be06
# ╠═ba0f5450-018f-11eb-2188-4327aa9c02cc
# ╟─d8c62ef0-018f-11eb-21ce-ebe037f9d91d
# ╠═e59a7c30-018f-11eb-2b94-7b8dff2a6383
# ╠═f5f705d0-018f-11eb-2462-e7e2265b700c
# ╟─1254d540-0190-11eb-3066-b5a853eef0c0
# ╠═02787820-0190-11eb-24ed-c5296c91df6b
# ╟─aa606210-018a-11eb-29bd-5d4f9ae0966e
# ╠═f3007460-018a-11eb-3fa2-8753aceaf89a
# ╟─24c67580-018b-11eb-3eb7-47a2e2ae5251
# ╠═a98c8890-0190-11eb-25a4-bfd18c42b749
# ╠═acfea4e0-0190-11eb-07a9-a3f38d97ff80
# ╠═a4824010-0190-11eb-0d97-a55353709552


text:
## Environment and packages
cd(@__DIR__)
using Pkg; Pkg.activate("."); Pkg.instantiate()

using OrdinaryDiffEq
using ModelingToolkit
using DataDrivenDiffEq
using LinearAlgebra, Optim
using DiffEqFlux, Flux
using Plots
gr()
using JLD2, FileIO
using Statistics
# Set a random seed for reproduceable behaviour
using Random
Random.seed!(1235)

# Create a name for saving ( basically a prefix )
svname = "Scenario_3_"

## Data Generation
# Parameters
D = 0.01f0; #diffusion
r = 1.0f0; #reaction rate

# Domain
X = 1.0f0; T = 5f0;
dx = 0.04f0; dt = T/10;
x = collect(0:dx:X);
t = collect(0:dt:T);
Nx = Int64(X/dx+1);
Nt = Int64(T/dt+1);

#initial conditions
Amp = 1.0f0;
Delta = 0.5f0;
rho_ic(Δ) = Amp*(tanh.((x .- (0.5f0 - Δ/2))/(Δ/10)) - tanh.((x .- (0.5f0 + Δ/2))/(Δ/10)))/2

# Reaction term
reaction(u) = r * u .* (1f0 .- u)
# Difference Op
lap = Float32.(diagm(0 => -2.0f0 * ones(Nx), 1=> ones(Nx-1), -1 => ones(Nx-1)) ./ dx^2)
#Periodic BC
lap[1,end] = 1.0f0/dx^2
lap[end,1] = 1.0f0/dx^2

function rc_ode(rho, p, t)
    #finite difference
    D * lap * rho + reaction.(rho)
end

# Generate random measurements
probpde = ODEProblem(rc_ode, rho_ic(Delta), (0.0f0, T), saveat=dt)
solution = solve(probpde, Tsit5());

# Ideal data
X = Array(solution)
t = solution.t

# Add noise in terms of the mean
x̄ = mean(X, dims = 2)
noise_magnitude = Float32(2.5e-2)
Xₙ = X .+ (noise_magnitude*x̄) .* randn(eltype(X), size(X))

pl_solution_1 = plot(solution, color = :black, ylabel = "rho(x,t)", label = ["True Data" [nothing for i in 1:Nx-1]...])
scatter!(t, transpose(Xₙ), color = :red,label = ["Noisy Data" [nothing for i in 1:Nx-1]...], legend =:bottomright)
# Create dataset plot
pl_contour = contour(t,x,Xₙ, xlabel = "t", ylabel = "x", fill = (true, :thermal))
pl_initial_data = plot(pl_solution_1, pl_contour, layout = (1,2))
savefig(pl_initial_data, joinpath(pwd(), "plots","$(svname)full_data_$(noise_magnitude).pdf"))

## Neural Network setup
## Define the network
# Gaussian RBF as activation
rbf(x) = exp.(-(x.^2))

# Create the network to recover the reaction
rx_nn = FastChain(
    FastDense(1, 5, rbf),
    FastDense(5, 5, rbf),
    FastDense(5, 5, rbf),
    FastDense(5, 1)
)

#initialize D0 close to D/dx^2
D0 = Float32[6.5]
p1s = initial_params(rx_nn)
p2s = zeros(Float32, 4)
p = [p1s;p2s;D0]

function nn_ode(u,p,t)
    # This is the dircetional derivative
    u_cnn_1   = [p[end-4] * u[end] + p[end-3] * u[1] + p[end-2] * u[2]]
    u_cnn     = [p[end-4] * u[i-1] + p[end-3] * u[i] + p[end-2] * u[i+1] for i in 2:size(u, 1)-1]
    u_cnn_end = [p[end-4] * u[end-1] + p[end-3] * u[end] + p[end-2] * u[1]]

    [rx_nn([ui], p[1:length(p1s)])[1] for ui in u] + p[end] * vcat(u_cnn_1, u_cnn, u_cnn_end)
end

prob_nn = ODEProblem(nn_ode, Xₙ[:, 1], (t[1], t[end]), p)

## Necessary functions

function predict_pde(θ, X = Xₙ[:, 1], T = t)
  # No ReverseDiff if using Flux
  Array(solve(prob_nn, Vern7(),
        u0 = X, p = θ, tspan = (T[1], T[end]), saveat=T,
        sensealg=ForwardDiffSensitivity(convert_tspan = false)#, reltol = 1e-4, abstol = 1e-4
        ))
end

#match data and force the weights of the CNN to add up to zero
function objective_pde(p)
    pred = predict_pde(p)
    sum(abs2, pred .- Xₙ) + abs(sum(p[end-4 : end-2]))
end

# Container to track the loss
losses = Float32[]

# Callback to see the training
callback(θ,l) = begin
    push!(losses, l)
    if length(losses)%1==0
        println("Current loss after $(length(losses)) iterations: $(losses[end])")
    end
    false
end

## Train the model
# First train with ADAM
res1 = DiffEqFlux.sciml_train(objective_pde, p, ADAM(0.1f0), cb=callback, maxiters = 40)
# Train with BFGS
res2 = DiffEqFlux.sciml_train(objective_pde, res1.minimizer, BFGS(initial_stepnorm=0.1f0), cb=callback, maxiters = 10000)
p_trained = res2.minimizer

## Data evaluation -> we subsample here!
tsample = t[1]:1.0*mean(diff(t)):t[end]
X̂ = predict_pde(p_trained, Xₙ[:,1], tsample)
# Trained on noisy data vs real solution
pl_trajectory = plot(tsample, transpose(X̂), xlabel = "t", ylabel ="rho(x,t)", color = :red,
 label = ["UDE Approximation" [nothing for i in 1:25]...])
scatter!(t, transpose(Xₙ), color = :black, label = ["Measurements"  [nothing for i in 1:25]...], legend = :bottomright)
savefig(pl_trajectory, joinpath(pwd(), "plots", "$(svname)_trajectory_reconstruction.pdf"))

# Create estimates
R̂ = similar(X̂)
R̄ = similar(X̂)
for i in 1:size(X̂, 1), j in 1:size(X̂, 2)
    R̂[i,j] = rx_nn(X̂[i,j], p_trained[1:length(p1s)])[1]
    R̄[i,j] = reaction(X̂[i,j])
end

pl_reconstruction = plot(tsample, R̂', label = ["UDE Approximation" [nothing for i in 1:25]...], color = :red, xlabel = "t", ylabel = "Reaction(x,t)")
plot!(tsample, R̄', label = ["True Reaction" [nothing for i in 1:25]...], color = :black)
pl_reconstruction_error = plot(tsample, norm.(eachcol(R̂ - R̄)), color = :red, xlabel = "t", ylabel = "L2 Error")
pl_missing = plot(pl_reconstruction, pl_reconstruction_error, layout = (2,1))
savefig(pl_missing, joinpath(pwd(), "plots", "$(svname)_missingterm_reconstruction_and_error.pdf"))
pl_overall = plot(pl_trajectory, pl_missing)
savefig(pl_overall, joinpath(pwd(), "plots", "$(svname)_reconstruction.pdf"))
## Symbolic regression via sparse regression / SINDy using DataDrivenDiffEq

# Reshape the reactions
@variables u[1:1]
b = Basis(monomial_basis(u, 10), u)
# We assume (via the modeling) a common mapping of the states
# via the reaction. In other words a global function. Hence, we can
# use all measurements as one variable
Xs = Matrix(vcat(eachrow(X̂[1:end, :])...)')
Rs = Matrix(vcat(eachrow(R̂[1:end, :])...)')

# Technically this is cheating until a general DataDrivenProblem is defined properly.
# we work with the continuous form
nn_prob = ContinuousDataDrivenProblem(Xs, zeros(Float32, size(Xs, 2)), DX = Rs)

λs = Float32.(exp10.(-3:0.01:5))
opt = STLSQ(λs)

nn_res = solve(nn_prob, b, opt, maxiters = 1000, progress = true, denoise = true, normalize = false)
# Print the output
println(nn_res)
println(result(nn_res))
println(parameters(nn_res))
## Save the results

save(joinpath(pwd(), "results" ,"$(svname)recovery_$(noise_magnitude).jld2"),
    "solution", solution, "X", Xₙ, "t" , t, "neural_network" , rx_nn, "initial_parameters", p, "trained_parameters" , p_trained, # Training
    "losses", losses, "result", nn_res,  # Recovery
    ) # Estimation

ext: jl


text:
#2020 Klaus Stricker - Tests for ReadWriteDlm2
#License is MIT: http://julialang.org/license

# rwd2tests_1.jl

# 1st block modified standardtests

isequaldlm(m1, m2, t) = isequal(m1, m2) && (eltype(m1) == eltype(m2) == t)

@testset "1_     readdlm2" begin
    @test isequaldlm(readdlm2(IOBuffer("1;2\n3;4\n5;6\n"), Float64), [1. 2; 3 4; 5 6], Float64)
    @test isequaldlm(readdlm2(IOBuffer("1;2\n3;4\n5;6\n"), Int), [1 2; 3 4; 5 6], Int)
    @test isequaldlm(readdlm2(IOBuffer("1 2\n3 4\n5 6\n"),' ', Float64), [1. 2; 3 4; 5 6], Float64)
    @test isequaldlm(readdlm2(IOBuffer("1 2\n3 4\n5 6\n"), ' ', Int), [1 2; 3 4; 5 6], Int)

    @test size(readdlm2(IOBuffer("1;2;3;4\n1;2;3"))) == (2,4)
    @test size(readdlm2(IOBuffer("1; 2;3;4\n1;2;3"))) == (2,4)
    @test size(readdlm2(IOBuffer("1; 2;3;4\n1;2;3\n"))) == (2,4)
    @test size(readdlm2(IOBuffer("1;;2;3;4\n1;2;3\n"))) == (2,5)

    let result1 = reshape(Any["", "", "", "", "", "", 1.0, 1.0, "", "", "", "", "", 1.0, 2.0, "", 3.0, "", "", "", "", "", 4.0, "", "", ""], 2, 13),
        result2 = reshape(Any[1.0, 1.0, 2.0, 1.0, 3.0, "", 4.0, ""], 2, 4)

        @test isequaldlm(readdlm2(IOBuffer(";;;1;;;;2;3;;;4;\n;;;1;;;1\n")), result1, Any)
        @test isequaldlm(readdlm2(IOBuffer("   1    2 3   4 \n   1   1\n"), ' '), result1, Any)
        @test isequaldlm(readdlm2(IOBuffer("   1    2 3   4 \n   1   1\n"), ' '), result1, Any)
        @test isequaldlm(readdlm2(IOBuffer("1;2\n3;4 \n"), Float64), [[1.0, 3.0] [2.0, 4.0]], Float64)
    end

    let result1 = reshape(Any["", "", "", "", "", "", "भारत", 1.0, "", "", "", "", "", 1.0, 2.0, "", 3.0, "", "", "", "", "", 4.0, "", "", ""], 2, 13)
        @test isequaldlm(readdlm2(IOBuffer(",,,भारत,,,,2,3,,,4,\n,,,1,,,1\n"), ',', rs=()) , result1, Any)
    end

    let result1 = reshape(Any[1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 4.0, ""], 2, 4)
        @test isequaldlm(readdlm2(IOBuffer("1;2;3;4\n1;2;3")), result1, Any)
        @test isequaldlm(readdlm2(IOBuffer("1;2;3;4\n1;2;3;"),';'), result1, Any)
        @test isequaldlm(readdlm2(IOBuffer("1,2,3,4\n1,2,3"), ',', rs=()), result1, Any)
        @test isequaldlm(readdlm2(IOBuffer("1,2,3,4\r\n1,2,3\r\n"), ',', rs=()), result1, Any)
        @test isequaldlm(readdlm2(IOBuffer("1,2,3,\"4\"\r\n1,2,3\r\n"), ',', rs=()), result1, Any)
    end

    let result1 = reshape(Any["abc", "hello", "def,ghi", " \"quote\" ", "new\nline", "world"], 2, 3),
        result2 = reshape(Any["abc", "line\"", "\"hello\"", "\"def", "", "\" \"\"quote\"\" \"", "ghi\"", "", "world", "\"new", "", ""], 3, 4)

        @test isequaldlm(readdlm2(IOBuffer("abc,\"def,ghi\",\"new\nline\"\n\"hello\",\" \"\"quote\"\" \",world"), ',', rs=()), result1, Any)
        @test isequaldlm(readdlm2(IOBuffer("abc,\"def,ghi\",\"new\nline\"\n\"hello\",\" \"\"quote\"\" \",world"), ',', rs=(), quotes=false), result2, Any)
        @test isequaldlm(readdlm2(IOBuffer("abc,\"def,ghi\",\"new\nline\"\n\"hello\",\" \"\"quote\"\" \",world"), ',', rs=()), result1, Any)
        @test isequaldlm(readdlm2(IOBuffer("abc,\"def,ghi\",\"new\nline\"\n\"hello\",\" \"\"quote\"\" \",world"), ',', rs=(), quotes=false), result2, Any)

        @test isequaldlm(readdlm2(IOBuffer("abc;\"def,ghi\";\"new\nline\"\n\"hello\";\" \"\"quote\"\" \";world")), result1, Any)
        @test isequaldlm(readdlm2(IOBuffer("abc;\"def;ghi\";\"new\nline\"\n\"hello\";\" \"\"quote\"\" \";world"), quotes=false), result2, Any)
        @test isequaldlm(readdlm2(IOBuffer("abc;\"def,ghi\";\"new\nline\"\n\"hello\";\" \"\"quote\"\" \";world"), rs=()), result1, Any)
        @test isequaldlm(readdlm2(IOBuffer("abc;\"def;ghi\";\"new\nline\"\n\"hello\";\" \"\"quote\"\" \";world"), rs=(), quotes=false), result2, Any)
    end

    let result1 = reshape(Any["t", "c", "", "c"], 2, 2),
        result2 = reshape(Any["t", "\"c", "t", "c"], 2, 2)
        @test isequaldlm(readdlm2(IOBuffer("t;\n\"c\";c")), result1, Any)
        @test isequaldlm(readdlm2(IOBuffer("t;t\n\"\"\"c\";c")), result2, Any)
        @test isequaldlm(readdlm2(IOBuffer("t;\n\"c\";c"), rs=()), result1, Any)
        @test isequaldlm(readdlm2(IOBuffer("t;t\n\"\"\"c\";c"), rs=()), result2, Any)
    end

    @test isequaldlm(readdlm2(IOBuffer("\n1;2;3\n4;5;6\n\n\n"), skipblanks=false),
                    reshape(Any["",1.0,4.0,"","","",2.0,5.0,"","","",3.0,6.0,"",""], 5, 3), Any)
    @test isequaldlm(readdlm2(IOBuffer("\n1;2;3\n4;5;6\n\n\n"), Float64, skipblanks=true),
                    reshape([1.0,4.0,2.0,5.0,3.0,6.0], 2, 3), Float64)
    @test isequaldlm(readdlm2(IOBuffer("1;2\n\n4;5"), skipblanks=false),
                    reshape(Any[1.0,"",4.0,2.0,"",5.0], 3, 2), Any)
    @test isequaldlm(readdlm2(IOBuffer("1;2\n\n4;5"), Float64, skipblanks=true),
                    reshape([1.0,4.0,2.0,5.0], 2, 2), Float64)
end

@testset "1_    writedlm2" begin
    let x = bitrand(5, 10), io = IOBuffer()
        writedlm2(io, x)
        seek(io, 0)
        @test readdlm2(io, Bool) == x
    end

    let x = bitrand(5, 10)
        writedlm2("test.csv", x, ";")
        @test readdlm2("test.csv", ';', Bool, '\n') == x
        @test readdlm2("test.csv", ';', '\n') == x
    end

    let x = bitrand(5, 10)
        writedlm2("test.csv", x, ";")
        @test readdlm2("test.csv", ';', Bool) == x
        @test readdlm2("test.csv", ';') == x
    end

    let x = bitrand(5, 10)
        writedlm2("test.csv", x)
        @test readdlm2("test.csv", Bool) == x
        @test readdlm2("test.csv") == x
    end

    let x = [1,2,3], y = [4,5,6], io = IOBuffer()
        writedlm2(io, zip(x,y), ",  ")
        seek(io, 0)
        @test readdlm2(io, ',', rs=()) == [x y]
    end

    let x = [0.1 0.3 0.5], io = IOBuffer()
        writedlm2(io, x, ", ")
        seek(io, 0)
        @test read(io, String) == "0,1, 0,3, 0,5\n"
    end

    let x = [0.1 0.3 0.5], io = IOBuffer()
        writedlm2(io, x, "; ")
        seek(io, 0)
        @test read(io, String) == "0,1; 0,3; 0,5\n"
    end

    let x = [0.1 0.3 0.5], io = IOBuffer()
        writedlm2(io, x)
        seek(io, 0)
        @test readdlm2(io) == [0.1 0.3 0.5]
    end

    let x = [0.1 0.3 0.5], io = IOBuffer()
        writedlm2(io, x, decimal='€')
        seek(io, 0)
        @test readdlm2(io, rs=(r"(\d)€(\d)", s"\1.\2")) == [0.1 0.3 0.5]
    end

    let x = [0.1 0.3 0.5], io = IOBuffer()
        writedlm2(io, x, ':', decimal='€')
        seek(io, 0)
        @test readdlm2(io, ':', rs=(r"(\d)€(\d)", s"\1.\2")) == [0.1 0.3 0.5]
    end

    let x = ["abc", "def\"ghi", "jk\nl"], y = [1, ",", "\"quoted\""], io = IOBuffer()
        writedlm2(io, zip(x,y))
        seek(io, 0)
        @test readdlm2(io) == [x y]
    end

    let x = ["abc", "def\"ghi", "jk\nl"], y = [1, ",", "\"quoted\""], io = IOBuffer()
        writedlm2(io, zip(x,y))
        seek(io, 0)
        @test readdlm2(io, Any) == [x y]
    end

    let x = ["a" "b"; "d" ""], io = IOBuffer()
        writedlm2(io, x)
        seek(io, 0)
        @test readdlm2(io) == x
    end

    let x = ["\"hello\"", "world\""], io = IOBuffer()
        writedlm2(io, x, quotes=false)
        @test String(take!(io)) == "\"hello\"\nworld\"\n"
        writedlm2(io, x)
        @test String(take!(io)) == "\"\"\"hello\"\"\"\n\"world\"\"\"\n"
    end
end

@testset "1_     comments" begin
    @test isequaldlm(readdlm2(IOBuffer("#this is comment\n1;2;3\n#one more comment\n4;5;6"), Float64, comments=true), [1. 2. 3.;4. 5. 6.], Float64)
    @test isequaldlm(readdlm2(IOBuffer("#this is \n#comment\n1;2;3\n#one more \n#comment\n4;5;6"), Float64, comments=true), [1. 2. 3.;4. 5. 6.], Float64)
    @test isequaldlm(readdlm2(IOBuffer("1;2;#3\n4;5;6"), comments=true), [1. 2. "";4. 5. 6.], Any)
    @test isequaldlm(readdlm2(IOBuffer("1#;2;3\n4;5;6"), comments=true), [1. "" "";4. 5. 6.], Any)
    @test isequaldlm(readdlm2(IOBuffer("1;2;\"#3\"\n4;5;6"), comments=true), [1. 2. "#3";4. 5. 6.], Any)
    @test isequaldlm(readdlm2(IOBuffer("1;2;3\n #with leading whitespace\n4;5;6"), comments=true), [1. 2. 3.;" " "" "";4. 5. 6.], Any)
end

@testset "1_    skipstart" begin
    x = ["a" "b" "c"; "d" "e" "f"; "g" "h" "i"; "A" "B" "C"; 1 2 3; 4 5 6; 7 8 9]
    io = IOBuffer()

    writedlm2(io, x, quotes=false)
    seek(io, 0)
    (data, hdr) = readdlm2(io, header=true, skipstart=3)
    @test data == [1 2 3; 4 5 6; 7 8 9]
    @test hdr == ["A" "B" "C"]

    x = ["a" "b" "\nc"; "d" "\ne" "f"; "g" "h" "i\n"; "A" "B" "C"; 1 2 3; 4 5 6; 7 8 9]
    io = IOBuffer()

    writedlm2(io, x, quotes=true)
    seek(io, 0)
    (data, hdr) = readdlm2(io, header=true, skipstart=6)
    @test data == [1 2 3; 4 5 6; 7 8 9]
    @test hdr == ["A" "B" "C"]

    io = IOBuffer()
    writedlm2(io, x, quotes=false)
    seek(io, 0)
    (data, hdr) = readdlm2(io, header=true, skipstart=6)
    @test data == [1 2 3; 4 5 6; 7 8 9]
    @test hdr == ["A" "B" "C"]
end

@testset "1_         i18n" begin
    # source: http://www.i18nguy.com/unicode/unicode-example-utf8.zip
    let i18n_data = ["Origin (English)", "Name (English)", "Origin (Native)", "Name (Native)",
            "Australia", "Nicole Kidman", "Australia", "Nicole Kidman",
            "Austria", "Johann Strauss", "Österreich", "Johann Strauß",
            "Belgium (Flemish)", "Rene Magritte", "België", "René Magritte",
            "Belgium (French)", "Rene Magritte", "Belgique", "René Magritte",
            "Belgium (German)", "Rene Magritte", "Belgien", "René Magritte",
            "Bhutan", "Gonpo Dorji", "འབྲུག་ཡུལ།", "མགོན་པོ་རྡོ་རྗེ།",
            "Canada", "Celine Dion", "Canada", "Céline Dion",
            "Canada - Nunavut (Inuktitut)", "Susan Aglukark", "ᓄᓇᕗᒻᒥᐅᑦ", "ᓱᓴᓐ ᐊᒡᓗᒃᑲᖅ",
            "Democratic People's Rep. of Korea", "LEE Sol-Hee", "조선 민주주의 인민 공화국", "이설희",
            "Denmark", "Soren Hauch-Fausboll", "Danmark", "Søren Hauch-Fausbøll",
            "Denmark", "Soren Kierkegaard", "Danmark", "Søren Kierkegård",
            "Egypt", "Abdel Halim Hafez", "ﻣﺼﺮ", "ﻋﺑﺪﺍﻠﺣﻟﻳﻢ ﺤﺎﻓﻅ",
            "Egypt", "Om Kolthoum", "ﻣﺼﺮ", "ﺃﻡ ﻛﻟﺛﻭﻡ",
            "Eritrea", "Berhane Zeray", "ብርሃነ ዘርኣይ", "ኤርትራ",
            "Ethiopia", "Haile Gebreselassie", "ኃይሌ ገብረሥላሴ", "ኢትዮጵያ",
            "France", "Gerard Depardieu", "France", "Gérard Depardieu",
            "France", "Jean Reno", "France", "Jean Réno",
            "France", "Camille Saint-Saens", "France", "Camille Saint-Saëns",
            "France", "Mylene Demongeot", "France", "Mylène Demongeot",
            "France", "Francois Truffaut", "France", "François Truffaut",
            "France (Braille)", "Louis Braille", "⠋⠗⠁⠝⠉⠑", "⠇⠕⠥⠊⠎⠀<BR>⠃⠗⠁⠊⠇⠇⠑",
            "Georgia", "Eduard Shevardnadze", "საქართველო", "ედუარდ შევარდნაძე",
            "Germany", "Rudi Voeller", "Deutschland", "Rudi Völler",
            "Germany", "Walter Schultheiss", "Deutschland", "Walter Schultheiß",
            "Greece", "Giorgos Dalaras", "Ελλάς", "Γιώργος Νταλάρας",
            "Iceland", "Bjork Gudmundsdottir", "Ísland", "Björk Guðmundsdóttir",
            "India (Hindi)", "Madhuri Dixit", "भारत", "माधुरी दिछित",
            "Ireland", "Sinead O'Connor", "Éire", "Sinéad O'Connor",
            "Israel", "Yehoram Gaon", "ישראל", "יהורם גאון",
            "Italy", "Fabrizio DeAndre", "Italia", "Fabrizio De André",
            "Japan", "KUBOTA Toshinobu", "日本", "久保田    利伸",
            "Japan", "HAYASHIBARA Megumi", "日本", "林原 めぐみ",
            "Japan", "Mori Ogai", "日本", "森鷗外",
            "Japan", "Tex Texin", "日本", "テクス テクサン",
            "Norway", "Tor Age Bringsvaerd", "Noreg", "Tor Åge Bringsværd",
            "Pakistan (Urdu)", "Nusrat Fatah Ali Khan", "پاکستان", "نصرت فتح علی خان",
            "People's Rep. of China", "ZHANG Ziyi", "中国", "章子怡",
            "People's Rep. of China", "WONG Faye", "中国", "王菲",
            "Poland", "Lech Walesa", "Polska", "Lech Wałęsa",
            "Puerto Rico", "Olga Tanon", "Puerto Rico", "Olga Tañón",
            "Rep. of China", "Hsu Chi", "臺灣", "舒淇",
            "Rep. of China", "Ang Lee", "臺灣", "李安",
            "Rep. of Korea", "AHN Sung-Gi", "한민국", "안성기",
            "Rep. of Korea", "SHIM Eun-Ha", "한민국", "심은하",
            "Russia", "Mikhail Gorbachev", "Россия", "Михаил Горбачёв",
            "Russia", "Boris Grebenshchikov", "Россия", "Борис Гребенщиков",
            "Slovenia", "\"Frane \"\"Jezek\"\" Milcinski", "Slovenija", "Frane Milčinski - Ježek",
            "Syracuse (Sicily)", "Archimedes", "Συρακούσα", "Ἀρχιμήδης",
            "Thailand", "Thongchai McIntai", "ประเทศไทย", "ธงไชย แม็คอินไตย์",
            "U.S.A.", "Brad Pitt", "U.S.A.", "Brad Pitt",
            "Yugoslavia (Cyrillic)", "Djordje Balasevic", "Југославија", "Ђорђе Балашевић",
            "Yugoslavia (Latin)", "Djordje Balasevic", "Jugoslavija", "Đorđe Balašević"]

        i18n_arr = permutedims(reshape(i18n_data, 4, Int(floor(length(i18n_data)/4))), [2, 1])
        i18n_buff = PipeBuffer()

        writedlm2(i18n_buff, i18n_arr, ',', decimal='.')
        @test i18n_arr == readdlm2(i18n_buff, ',', rs=())

        writedlm2(i18n_buff, i18n_arr)
        @test i18n_arr == readdlm2(i18n_buff)

        hdr = i18n_arr[1:1, :]
        data = i18n_arr[2:end, :]

        writedlm2(i18n_buff, i18n_arr)
        @test (data, hdr) == readdlm2(i18n_buff, header=true)
    end
end

@testset "1_ issue #13028" begin
    for data in ["A B C", "A B C\n"]
        data,hdr = readdlm2(IOBuffer(data), ' ', header=true)
        @test hdr == AbstractString["A" "B" "C"]
        @test data == Array{Float64}(undef, 0, 3)
    end
end

@testset "1_ issue #13179" begin # fix #13179 parsing unicode lines with default delmiters
    @test isequaldlm(readdlm2(IOBuffer("# Should ignore this π\n1\tα\n2\tβ\n"), '\t', comments=true), Any[1 "α"; 2 "β"], Any)
    @test isequaldlm(readdlm2(IOBuffer("# Should ignore this π\n1;α\n2;β\n"), comments=true), Any[1 "α"; 2 "β"], Any)
    @test isequaldlm(readdlm2(IOBuffer("# Should ignore this π\n1\tα\n2\tβ\n"), '\t', rs=(), comments=true), Any[1 "α"; 2 "β"], Any)
    @test isequaldlm(readdlm2(IOBuffer("# Should ignore this π\n1;α\n2;β\n"), rs=(), comments=true), Any[1 "α"; 2 "β"], Any)
end

@testset "1_ other issues" begin
    # BigInt parser
    let data = "1;2;3"
        @test readdlm2(IOBuffer(data), BigInt) == BigInt[1 2 3]
        @test readdlm2(IOBuffer(data), BigInt, rs=()) == BigInt[1 2 3]
    end
    let data = "1 2 3"
        @test readdlm2(IOBuffer(data), ' ', BigInt) == BigInt[1 2 3]
        @test readdlm2(IOBuffer(data), ' ', BigInt, rs=()) == BigInt[1 2 3]
    end

    # Test that we can read a write protected file
    let fn = tempname()
        open(fn, "w") do f
            write(f, "Julia")
        end
        chmod(fn, 0o444)
        @test readdlm2(fn)[] == "Julia"
        rm(fn)
    end

    # issue #21180
    let data = "\"721\",\"1438\",\"1439\",\"…\",\"1\""
        @test readdlm2(IOBuffer(data), ',', rs=()) == Any[721  1438  1439  "…"  1]
    end
    let data = "\"721\";\"1438\";\"1439\";\"…\";\"1\""
        @test readdlm2(IOBuffer(data)) == Any[721  1438  1439  "…"  1]
    end

    # issue #21207
    let data = "\"1\";\"灣\"\"灣灣灣灣\";\"3\""
        @test readdlm2(IOBuffer(data)) == Any[1 "灣\"灣灣灣灣" 3]
    end

    # issue #11484: useful error message for invalid readdlm filepath arguments
    # not implemented in ReadWriteDlm2
end

@testset "1_      complex" begin
    @test readdlm2(IOBuffer("3+4im; 4+5im"), Complex{Int}) == [3+4im 4+5im]
end


text:
using LinearAlgebra, Yao.ConstGate
using Test, Random
using CuYao
using StaticArrays
using Yao.ConstGate: SWAPGate
using CuArrays

@testset "gpu instruct nbit!" begin
    Random.seed!(3)
    nbit = 6
    N = 1<<nbit
    LOC1 = SVector{2}([0, 1])
    v1 = randn(ComplexF32, N)
    vn = randn(ComplexF32, N, 333)

    for UN in [
            rand_unitary(ComplexF32, 4),
            mat(ComplexF32, CNOT),
            mat(ComplexF32, control(2,2,1=>Z)),
            mat(ComplexF32, put(2,2=>I2)),
            mat(ComplexF32, put(2,2=>P0))
            ]

        @test instruct!(v1 |> CuArray, UN, (3,1)) |> Vector ≈ instruct!(v1 |> copy, UN, (3,1))
        @test instruct!(vn |> CuArray, UN, (3,1)) |> Matrix ≈ instruct!(vn |> copy, UN, (3,1))
    end
    # sparse matrix like P0, P1 et. al. are not implemented.
end

@testset "gpu swapapply!" begin
    nbit = 6
    N = 1<<nbit
    LOC1 = SVector{2}([0, 1])
    v1 = randn(ComplexF32, N)
    vn = randn(ComplexF32, N, 333)

    @test instruct!(v1 |> CuArray, Val(:SWAP), (3,5)) |> Vector ≈ instruct!(v1 |> copy, Val(:SWAP), (3,5))
    @test instruct!(vn |> CuArray, Val(:SWAP), (3,5)) |> Matrix ≈ instruct!(vn |> copy, Val(:SWAP), (3,5))
end

@testset "gpu instruct! 1bit" begin
    nbit = 6
    N = 1<<nbit
    LOC1 = SVector{2}([0, 1])
    v1 = randn(ComplexF32, N)
    vn = randn(ComplexF32, N, 333)

    for U1 in [mat(H), mat(Y), mat(Z), mat(I2), mat(P0)]
        @test instruct!(v1 |> CuArray, U1, (3,)) |> Vector ≈ instruct!(v1 |> copy, U1, (3,))
        @test instruct!(vn |> CuArray, U1, (3,)) |> Matrix ≈ instruct!(vn |> copy, U1, (3,))
    end
    # sparse matrix like P0, P1 et. al. are not implemented.
end

@testset "gpu xyz-instruct!" begin
    nbit = 6
    N = 1<<nbit
    LOC1 = SVector{2}([0, 1])
    v1 = randn(ComplexF32, N)
    vn = randn(ComplexF32, N, 333)

    for G in [:X, :Y, :Z, :T, :H, :Tdag, :S, :Sdag]
        @test instruct!(v1 |> CuArray, Val(G), (3,)) |> Vector ≈ instruct!(v1 |> copy, Val(G), (3,))
        @test instruct!(vn |> CuArray, Val(G), (3,)) |> Matrix ≈ instruct!(vn |> copy, Val(G), (3,))
        if G != :H
            @test instruct!(v1 |> CuArray, Val(G), (1,3,4)) |> Vector ≈ instruct!(v1 |> copy, Val(G), (1,3,4))
            @test instruct!(vn |> CuArray,  Val(G),(1,3,4)) |> Matrix ≈ instruct!(vn |> copy, Val(G), (1,3,4))
        end
    end
end

@testset "gpu cn-xyz-instruct!" begin
    nbit = 6
    N = 1<<nbit
    LOC1 = SVector{2}([0, 1])
    v1 = randn(ComplexF32, N)
    vn = randn(ComplexF32, N, 333)

    for G in [:X, :Y, :Z, :T, :Tdag, :S, :Sdag]
        @test instruct!(v1 |> CuArray, Val(G), (3,), (4,5), (0, 1)) |> Vector ≈ instruct!(v1 |> copy, Val(G), (3,), (4,5), (0, 1))
        @test instruct!(vn |> CuArray, Val(G), (3,), (4,5), (0, 1)) |> Matrix ≈ instruct!(vn |> copy, Val(G), (3,), (4,5), (0, 1))
        @test instruct!(v1 |> CuArray, Val(G), (3,), (1,), (1,)) |> Vector ≈ instruct!(v1 |> copy, Val(G),(3,), (1,), (1,))
        @test instruct!(vn |> CuArray, Val(G), (3,), (1,), (1,)) |> Matrix ≈ instruct!(vn |> copy, Val(G),(3,), (1,), (1,))
    end
end

@testset "pswap" begin
    ps = put(6, (2,4)=>rot(SWAP, π/2))
    reg = rand_state(6; nbatch=10)
    @test apply!(reg |> cu, ps) |> cpu ≈ apply!(copy(reg), ps)
    @test apply!(reg |> cu, ps).state isa CuArray
end


text:
module MethodInSrc

"""
    module MethodInSrc

`MethodInSrc` provides tools to use in your test suite that verify that
a function call dispatches (or does not) to a method defined in your module.
They are meant to verify that a specialized method is called rather than a more generic one. (Or vice versa.)

This module exports `@isinsrc`, `@insrc`, `@ninsrc`,
`@isinmodule`, `@inmodule`, and  `@ninmodule`.

`@isinsrc f(x)` returns `true` if the method for `f(x)` is found under `../src`. But, it does not evaluate `f(x)`.

`@insrc f(x)` throws an error if the method for `f(x)` is not found under `../src`. Otherwise, it evaluates `f(x)`.

`@ninsrc` is the same as `@insrc` except that it throws if the method *is* under `../src`.

The first three macros are meant to be used in files in your module's "test" directory.
In this case, `@isinsrc f([x,...])` returns `true` if the method that would be called by `f([x,...])` was defined
in the module's "src" directory.

See the doc strings for `@isinmodule`, `@inmodule`, and  `@ninmodule`.

## Examples

Suppose the type `MyPackage.AMatrix` represents a square matrix whose elements are all equal to `1`.
`MyPackage` extends `Base.sum` with an efficient method for `::AMatrix`.
`MyPackage` also includes an efficient function `prod(::AMatrix)`.
But, we neglected to write `Base.prod` or `import Base: prod`.
So `MyPackage.prod` is not an extension of `Base.prod`.

### `@isinsrc`

Use `@isinsrc` to
check that both `sum` and `prod` extend `Base` functions. (`prod` does not!)
```julia
using MyPackage
using MethodInSrc
using Test

m = MyPackage.AMatrix{Int}(3)
@test @isinsrc sum(m)
@test @isinsrc prod(m)  # This will fail!
```

### `@insrc`, `@ninsrc`

Use `@insrc` if you are too lazy to write two tests,
one to verify that you have the correct method,
and another to test its correctness.

`@insrc` verifies that the methods are defined in "../src",
and then evaluate the expressions.
Here, it is assumed we know that `MyPackage.prod` is a different function and test for it.
These tests all pass.
```julia
using MyPackage
using MethodInSrc
using Test

N = 3
m = MyPackage.AMatrix{Int}(N)

# Note that `@insrc` takes the next expression as an argument.
# We need to test just the function call.
@test_throws ErrorException 1 == @insrc prod(m)   # This finds the method.
@test_throws ErrorException @insrc(prod(m)) == 1  # This does too.
@test 1 == @ninsrc prod(m)  # Do the test if the method is generic

# The following methods are found in "../src", so the expressions are evaluated
@test N^2 == @insrc sum(m)
@test (@insrc sum(m)) == N^2
@test @insrc(MyPackage.prod(m)) == 1
```
"""
MethodInSrc

using InteractiveUtils

export @isinsrc, @insrc, @ninsrc, @inmodule, @ninmodule, @isinmodule

"""
    findmethod(call::Expr)

Return the `Method` invoked by `call`.
"""
findmethod(call) = InteractiveUtils.gen_call_with_extracted_types(@__MODULE__, :which, call)

srcdir(testdir) = abspath(joinpath(testdir, "..", "src"))
srcdir() = srcdir(pwd())

"""
    methoddir(meth::Method)

Return the path to the directory containing the source for `meth`.
For methods in `Base`, return the empty string.
"""
methoddir(meth::Method) = dirname(string(meth.file))

methpath(meth::Method) = string(meth.file)

"""
    issubdir(dir::AbstractString, subdir::AbstractString)

Return `true` if `subdir` is a proper subdirectory of `dir`.
"""
function issubdir(dir::AbstractString, subdir::AbstractString)
    dirparts = splitpath(normpath(dir))
    subdirparts = splitpath(normpath(subdir))
    length(subdirparts) < length(dirparts) && return false
    subdirparts[1:length(dirparts)] == dirparts && return true
    return false
end

"""
    isdirorsubdir(dir::AbstractString, subdir::AbstractString)

Return `true` if `subdir` is equal to, or a subdirectory of, `dir`.
"""
function isdirorsubdir(dir::AbstractString, subdir::AbstractString)
    dir == subdir && return true
    return issubdir(dir, subdir)
end

"""
    @isinsrc f([x, y,...])

Return true if the method that would be called for the expression is
in the source directory ("src") of the module in which the test running.

`@isinsrc` is meant to be called from within a module's "test" directory.
The src directory is then found via "test/../src".

### Example
Verify that the method for `sum` is defined in `MyPackage`.
```julia
using MyPackage
using MethodInSrc
using Test

m = MyPackage.AMatrix(3)
@test @isinsrc sum(m)
```
"""
macro isinsrc(call)
    quote
        isdirorsubdir(srcdir(), methoddir($(findmethod(call))))
    end
end

"""
    @insrc f([x, y,...])

Evaluate the expression if the method called is defined in the source directory ("src")
of the module in which the test running. Otherwise, throw an `ErrorException`.

`@insrc` is meant to be called from within a module's "test" directory.
The src directory is then found via "test/../src".

### Example
Verify that the method for `sum` is defined in `MyPackage`.
```julia
using MyPackage
using MethodInSrc
using Test

m = MyPackage.AMatrix(3)
@test @insrc(sum(m)) == 9
```
"""
macro insrc(call)
    meth = findmethod(call)
    quote
        src_dir = srcdir()
        methdir = methoddir($meth)
        if ! isdirorsubdir(src_dir, methdir)
            throw(ErrorException("@insrc: Method defined in '$(methpath($meth))', not in '$src_dir'."))
        end
        $(esc(call))
    end
end

"""
    @ninsrc f([x, y,...])

The same as `@insrc` except that `@ninsrc` throws an error if the method
*is* in "../src", while `@insrc` throws if the method is *not* in "../src".
"""
macro ninsrc(call)
    meth = findmethod(call)
    quote
        src_dir = srcdir()
        methdir = methoddir($meth)
        if isdirorsubdir(src_dir, methdir)
            throw(ErrorException("@ninsrc: Method defined in '$(methpath($meth))', which is a subdirectory of '$src_dir'."))
        end
        $(esc(call))
    end
end

"""
    @isinmodule ModuleName f([x, y,...])

Return true if the method that would be called for the expression is
in the source directory ("src") of `ModuleName`
"""
macro isinmodule(mod, call)
    quote
        isdirorsubdir(dirname(pathof($(esc(mod)))), methoddir($(findmethod(call))))
    end
end

"""
    @inmodule ModuleName f([x, y,...])

Evaluate the the expressoin `f...` if the method called is defined in the
source directory of `ModuleName`. Otherwise throw an ErrorException.
"""
macro inmodule(mod, call)
    meth = findmethod(call)
    quote
        src_dir = dirname(pathof($(esc(mod))))
        methdir = methoddir($meth)
        if ! isdirorsubdir(src_dir, methdir)
            throw(ErrorException("@inmodule: Method defined in '$(methpath($meth))', not in '$src_dir'."))
        end
        $(esc(call))
    end
end

"""
    @ninmodule ModuleName f([x, y,...])

Evaluate the the expressoin `f...` if the method called is *not* defined in the
source directory of `ModuleName`. Otherwise throw an ErrorException.
"""
macro ninmodule(mod, call)
    meth = findmethod(call)
    quote
        src_dir = dirname(pathof($(esc(mod))))
        methdir = methoddir($meth)
        if isdirorsubdir(src_dir, methdir)
            throw(ErrorException("@ninmodule: Method defined in '$(methpath($meth))', which is a subdirectory of '$src_dir'."))
        end
        $(esc(call))
    end
end

# include code defining a type and methods for use in testing this module
include("testmethods.jl")
include("subdir/method_in_subdir.jl")

end # module


text:
using Pkg;

##Install main packages
Pkg.add([
    "MLJ",
    "Query",
    "VegaLite",
    "CategoricalArrays",
    "DataFrames",
    "Combinatorics",
    "HypothesisTests",
    "CSV",
    "PrettyPrinting",
    "Missings",
    "StatsBase",
    "Tables",
    "Hwloc",
    "ShapML",
    "PkgTemplates",
    "FreqTables",
    "Distributions",
    "IJulia"
])
Pkg.status()

##Install all MLJ relevant packages
Pkg.add([
    "Clustering",
    "DecisionTree",
    "EvoTrees",
    "GLM",
    "LightGBM",
    "LIBSVM",
    "MLJModels",
    "MLJLinearModels",
    "XGBoost",
    "ParallelKMeans"
])
Pkg.status()

#  add MLJ@0.12
#  using Pkg; Pkg.add(PackageSpec(url = "https://github.com/nredell/ShapML.jl"))
#  C:\Users\BCP\.julia\environments\v1.4


text:
struct Chain_ID
    obs_id::AbstractString
    K::Integer
    order::Integer
    replicate::Integer
    Chain_ID(obs_id,K,order,replicate) = assert_chain_id(K,order,replicate) && new(obs_id,K,order,replicate)
end

function assert_chain_id(K, order, replicate)
    K < 1 && throw(ArgumentError("Chain_ID K (# of hmm states) must be a positive integer!"))
    order < 0 && throw(ArgumentError("Chain_ID symbol order must be zero or a positive integer!"))
    replicate < 1 && throw(ArgumentError("Chain_ID replicate must be a positive integer!"))
    return true
end

struct EM_step
    iterate::Integer
    hmm::BHMM
    log_p::AbstractFloat
    delta::AbstractFloat
    converged::Bool
    EM_step(iterate,hmm,log_p,delta,converged)=assert_step(iterate,hmm,log_p) && new(iterate, hmm, log_p, delta, converged)
end

function assert_step(iterate, hmm, log_p)
    iterate < 1 && throw(ArgumentError("EM_step iterate number must be a positive integer!"))
    assert_hmm(hmm.a, hmm.A, hmm.B)
    log_p > 0.0 && throw(ArgumentError("EM_step log probability value must be 0 or negative!"))
    return true
end


text:
mutable struct Spmat{VTI<:AbstractVector, VTF<:AbstractVector}
    colptr::VTI
    rowval::VTI
    nzval::VTF

    # create 2 Spmats from complex matrix
    function Spmat{VTI, VTF}(mat::SparseMatrixCSC{Complex{Float64}, Int}) where {VTI, VTF}
        matreal = new(VTI(mat.colptr), VTI(mat.rowval), VTF(real.(mat.nzval)))
        matimag = new(VTI(mat.colptr), VTI(mat.rowval), VTF(imag.(mat.nzval)))
        return matreal, matimag
    end
    # copy constructor
    function Spmat{VTI, VTF}(mat) where {VTI, VTF}
        return new(VTI(mat.colptr), VTI(mat.rowval), VTF(mat.nzval))
    end
end


text:
module ReagentsTests

include("utils.jl")
include("test_immutablelists.jl")
include("test_internal_bags.jl")
include("test_treiberstack.jl")
include("test_msqueue.jl")
include("test_dualcontainers.jl")
include("test_channels.jl")
include("test_cancellablecontainers.jl")
include("test_anchors.jl")
include("test_doctest.jl")
include("test_aqua.jl")

end  # module ReagentsTests


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

using Markdown
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 4a816b10-c188-11ea-3935-1fc7371387ac
using Plots

# ╔═╡ 62a8f092-c198-11ea-2d26-632729a6bd7f
using StatsBase

# ╔═╡ 1c0741f0-c17a-11ea-3c6f-31f0c86dd08d
md"""
# Random Walk Models
"""

# ╔═╡ 58e6d60e-c17c-11ea-145e-ffdc97ab1f3f
md"Short function definition (one-liner):"

# ╔═╡ a93d6520-c17c-11ea-2cc7-5ff51f20df2b
jump() = rand( (-1, +1) )

# ╔═╡ bfa68c10-c17c-11ea-281c-25e48e4ab81d
md"Exactly the same"

# ╔═╡ ccef6d10-c17c-11ea-13e5-432b940ba726
function jump_()
	return rand( (-1, +1) )
end

# ╔═╡ fea8c400-c17c-11ea-35fc-e98ef849c316
jump()

# ╔═╡ 9aff5fa0-c185-11ea-145a-c3230681bd4c
md"It is better to use `tuples` in this case instead of `arrays`, tuples are not mutable and with fixed size"

# ╔═╡ 0762ce10-c17d-11ea-1f2a-172e4ed0500a
begin
	local t = (-1, +1)
	typeof(t)
end

# ╔═╡ 212df5a0-c186-11ea-0f5c-83f0dd5a1497
md"Side remark: Simple bernoulli function"

# ╔═╡ 8ce720ae-c185-11ea-380c-2df4848fb56a
bernoulli(p) = rand() < p

# ╔═╡ f52fc0a0-c185-11ea-2126-b10b11defe13
bernoulli(0.25)

# ╔═╡ 891b2f20-c186-11ea-024f-372f0eef600c
md"## Random walks"

# ╔═╡ 15990630-c186-11ea-2059-e1e301040dae
[jump() for i ∈ 1:10]

# ╔═╡ 62d228a2-c186-11ea-3ced-1db44003e269
function walk(n)
	x = 0

	for i ∈ 1:n
		x += jump() #x = x + jump()
	end

	return x
end

# ╔═╡ 0b66e5a2-c187-11ea-36de-279c4d6aa180
walk(10)

# ╔═╡ 432008a0-c187-11ea-3606-159ae4ed8460
function trajectory(n)
	x = 0

	xs = Int64[x] # inititalize with an initial value and the type of the array

	for i ∈ 1:n
		x += jump() # x = x + jump()

		push!(xs, x)
	end

	return xs
end

# ╔═╡ 131fcc00-c18a-11ea-09dc-0906c430d399
md"""steps = $(@bind steps html"<input type='range' min=5 max=100>" ) """

# ╔═╡ 1afa1280-c18c-11ea-1d1b-7f277d76e68c
md"**$(steps)** steps"

# ╔═╡ 8cfcb9a2-c187-11ea-1129-ff3d55e25de4
trajectories = trajectory(steps)

# ╔═╡ 3b677390-c188-11ea-1079-01f262b1d76f
md"## Plotting"

# ╔═╡ 643c5510-c188-11ea-2d15-fbf4939ae53c
begin
	plot(trajectories, m=:o, size=(300, 200))
	hline!([0], ls=:dash)
end

# ╔═╡ b417eaf0-c18c-11ea-3ae7-19f88eab6b9f
begin
	local num_walkers = 1000
	local num_steps = 100

	local p = plot(size=(500, 400), leg=false)

	for i ∈ 1:num_walkers
		traj = trajectory(num_steps)

		plot!(traj)
	end

	p
end

# ╔═╡ eb35c2e0-c18d-11ea-3209-bd013b156820
md"## Precompute the data"

# ╔═╡ 1ec51ce0-c18f-11ea-2166-6923cbcef091
n = 20

# ╔═╡ 2c106ee0-c18f-11ea-21dc-db7d6706f5cb
traj = trajectory(n)

# ╔═╡ 3ee28c60-c18f-11ea-3297-3d26c8d37627
@bind traj_steps html"<input type='range' min=1 max=20>"

# ╔═╡ 1ac71ed0-c190-11ea-2cfe-4d7cfd514800
md"`traj_steps: ` $(traj_steps)"

# ╔═╡ 727357e0-c18e-11ea-24f4-81a80c9e5f36
begin
	plot(traj[1:traj_steps])
	xlims!(0, length(traj))
	ylims!(minimum(traj), maximum(traj))
end

# ╔═╡ 55e2ff60-c191-11ea-1b6b-b9a74d8fab5f
md"## Run N walkers simultaneously"

# ╔═╡ 60381bd0-c191-11ea-36a6-69333299f23b
md"Rule: NEVER copy and paste (excpet occasionally)

Instead: Automate "

# ╔═╡ e1c193c0-c191-11ea-0c16-67d6d7611c26
begin
	local n = 20

	traj_ = [trajectory(n)]

	push!(traj_, trajectory(n))
end

# ╔═╡ 6fad5d40-c192-11ea-2573-01f0eb6c17ed
md"`append!` extracts the elements of a vector"

# ╔═╡ f4e1b4c0-c192-11ea-3782-79f7c856869f
begin
	local v = [1, 2, 3 ,4]

	append!(v, [5, 6])
end

# ╔═╡ 0a1011c2-c193-11ea-04a2-e71ff37bf002
md"a vector of vectors using an array comprehension"

# ╔═╡ d8b930b0-c193-11ea-0d53-63e3111f83a3
n_ = 20

# ╔═╡ c8dad7c0-c193-11ea-04a6-cf21953820d6
walkers = [trajectory(n_) for i ∈ 1:100]

# ╔═╡ 0f2b2db0-c194-11ea-06a1-6b2192a8e677
md"Use `end` to get the last element in a vector"

# ╔═╡ 61d6a940-c194-11ea-1d9f-4d6550479b45
[ traj_[end] for traj_ in walkers ] # iterating through walkers, a Vector of Vectors

# ╔═╡ 02bf22b0-c195-11ea-2d1b-ddebd9d7a069
md"Convert Vector of Vectors into a matrix:"

# ╔═╡ 13a7bfb0-c195-11ea-13d2-41437b78aae1
hcat(walkers...) # "splat" or like the "spread" operator in js

# ╔═╡ 188fbdb0-c196-11ea-16c2-0bbd3bd3c5f2
md"a more efficient way:"

# ╔═╡ b3074f30-c195-11ea-0025-af543bf613be
walkers_matrix = reduce(hcat, walkers)

# ╔═╡ 2a00013e-c196-11ea-1d45-3d8ff26e1e53
md"las element of the first row:"

# ╔═╡ f2221830-c195-11ea-09ff-2fb3e85bde34
walkers_matrix[1, end]

# ╔═╡ 5ed8c0a0-c196-11ea-2e0c-370732384c6b
md"last element of the first and second rows:"

# ╔═╡ 43f2c420-c196-11ea-3bfa-3320a1f4f92e
walkers_matrix[1:2, end]

# ╔═╡ 6e5598a0-c196-11ea-11d2-8f08d13fb07f
md"get the last element of all of the rows:"

# ╔═╡ 8b5d1592-c196-11ea-2a9f-17fbc031a2bd
walkers_matrix[1:end, end]

# ╔═╡ 9cbb09a0-c196-11ea-0975-c59ab691f596
md"also:"

# ╔═╡ a19bf82e-c196-11ea-024c-c7de6f0fbdb2
final_positions = walkers_matrix[:, end]

# ╔═╡ aed262a0-c196-11ea-310c-976af5163b2a
md"## Where are the walkers at last time step?"

# ╔═╡ e3dc4420-c196-11ea-37be-1388e05cdfc5
scatter(final_positions)

# ╔═╡ 8e8a1d20-c197-11ea-2241-7fed80778896
begin
	local n = 100

	local num_walkers = 10_000

	local walkers = [trajectory(n) for i ∈ 1:num_walkers]

	final_positions_ = [traj_[end] for traj_ ∈ walkers]
end

# ╔═╡ a46bcf22-c198-11ea-3496-cfdbdfcac559
scatter(final_positions_)

# ╔═╡ 527b1810-c198-11ea-3634-5f82562c0394
import Pkg; Pkg.add("StatsBase")

# ╔═╡ 76ce6910-c198-11ea-1a99-e33a29b2860f
counts = countmap(final_positions_)

# ╔═╡ c6170a90-c198-11ea-19c6-a377288360bd
scatter(counts)

# ╔═╡ Cell order:
# ╟─1c0741f0-c17a-11ea-3c6f-31f0c86dd08d
# ╟─58e6d60e-c17c-11ea-145e-ffdc97ab1f3f
# ╠═a93d6520-c17c-11ea-2cc7-5ff51f20df2b
# ╟─bfa68c10-c17c-11ea-281c-25e48e4ab81d
# ╠═ccef6d10-c17c-11ea-13e5-432b940ba726
# ╠═fea8c400-c17c-11ea-35fc-e98ef849c316
# ╟─9aff5fa0-c185-11ea-145a-c3230681bd4c
# ╠═0762ce10-c17d-11ea-1f2a-172e4ed0500a
# ╟─212df5a0-c186-11ea-0f5c-83f0dd5a1497
# ╠═8ce720ae-c185-11ea-380c-2df4848fb56a
# ╠═f52fc0a0-c185-11ea-2126-b10b11defe13
# ╟─891b2f20-c186-11ea-024f-372f0eef600c
# ╠═15990630-c186-11ea-2059-e1e301040dae
# ╠═62d228a2-c186-11ea-3ced-1db44003e269
# ╠═0b66e5a2-c187-11ea-36de-279c4d6aa180
# ╠═432008a0-c187-11ea-3606-159ae4ed8460
# ╟─131fcc00-c18a-11ea-09dc-0906c430d399
# ╟─1afa1280-c18c-11ea-1d1b-7f277d76e68c
# ╠═8cfcb9a2-c187-11ea-1129-ff3d55e25de4
# ╟─3b677390-c188-11ea-1079-01f262b1d76f
# ╠═4a816b10-c188-11ea-3935-1fc7371387ac
# ╠═643c5510-c188-11ea-2d15-fbf4939ae53c
# ╠═b417eaf0-c18c-11ea-3ae7-19f88eab6b9f
# ╟─eb35c2e0-c18d-11ea-3209-bd013b156820
# ╠═1ec51ce0-c18f-11ea-2166-6923cbcef091
# ╠═2c106ee0-c18f-11ea-21dc-db7d6706f5cb
# ╟─3ee28c60-c18f-11ea-3297-3d26c8d37627
# ╟─1ac71ed0-c190-11ea-2cfe-4d7cfd514800
# ╠═727357e0-c18e-11ea-24f4-81a80c9e5f36
# ╟─55e2ff60-c191-11ea-1b6b-b9a74d8fab5f
# ╟─60381bd0-c191-11ea-36a6-69333299f23b
# ╠═e1c193c0-c191-11ea-0c16-67d6d7611c26
# ╟─6fad5d40-c192-11ea-2573-01f0eb6c17ed
# ╠═f4e1b4c0-c192-11ea-3782-79f7c856869f
# ╟─0a1011c2-c193-11ea-04a2-e71ff37bf002
# ╠═d8b930b0-c193-11ea-0d53-63e3111f83a3
# ╠═c8dad7c0-c193-11ea-04a6-cf21953820d6
# ╟─0f2b2db0-c194-11ea-06a1-6b2192a8e677
# ╠═61d6a940-c194-11ea-1d9f-4d6550479b45
# ╟─02bf22b0-c195-11ea-2d1b-ddebd9d7a069
# ╠═13a7bfb0-c195-11ea-13d2-41437b78aae1
# ╟─188fbdb0-c196-11ea-16c2-0bbd3bd3c5f2
# ╠═b3074f30-c195-11ea-0025-af543bf613be
# ╟─2a00013e-c196-11ea-1d45-3d8ff26e1e53
# ╠═f2221830-c195-11ea-09ff-2fb3e85bde34
# ╟─5ed8c0a0-c196-11ea-2e0c-370732384c6b
# ╠═43f2c420-c196-11ea-3bfa-3320a1f4f92e
# ╟─6e5598a0-c196-11ea-11d2-8f08d13fb07f
# ╠═8b5d1592-c196-11ea-2a9f-17fbc031a2bd
# ╟─9cbb09a0-c196-11ea-0975-c59ab691f596
# ╠═a19bf82e-c196-11ea-024c-c7de6f0fbdb2
# ╟─aed262a0-c196-11ea-310c-976af5163b2a
# ╠═e3dc4420-c196-11ea-37be-1388e05cdfc5
# ╠═8e8a1d20-c197-11ea-2241-7fed80778896
# ╠═a46bcf22-c198-11ea-3496-cfdbdfcac559
# ╠═527b1810-c198-11ea-3634-5f82562c0394
# ╠═62a8f092-c198-11ea-2d26-632729a6bd7f
# ╠═76ce6910-c198-11ea-1a99-e33a29b2860f
# ╠═c6170a90-c198-11ea-19c6-a377288360bd

max_stars_repo_name: elihugarret/my-notebooks


text:
## Time discretization and Solution Post Processing
#
# Currently implements fourth order Runge Kutta for CG and DG methods.
#
# Additional features include:
#   - MUSCL Limiter for (currently) 1D, p=1 problems
#
#
import NumericExtensions: add!, multiply!

function globalSolve{T}(p::DG,MMatrix,Q::Array{T},t,cache)
    k = Array(T,1,numcells(p),nbf(p))
    for c in p.mesh
        k[1,cid(c),:] = MMatrix[cid(c)]\RHS(p,c,Q,t,cache)
    end
    k
end

LHS(p::Galerkin,basis) = [factorize(Mlocal(p,c,basis)) for c in mesh(p)]

# RK4 Solver
function _solve(p::Galerkin, ℚ; cache = generateMatrices(p))
    basis = phi(p)
    # ℚ[1] is the inital condition
    M = LHS(p,basis)
    for i in 2:p.nt
        k1 = globalSolve(p,M,ℚ[i-1,:,:],i*p.Δt,cache)
        k2 = globalSolve(p,M,ℚ[i-1,:,:] + 0.5*p.Δt*k1,i*p.Δt,cache)
        k3 = globalSolve(p,M,ℚ[i-1,:,:] + 0.5*p.Δt*k2,i*p.Δt,cache)
        k4 = globalSolve(p,M,ℚ[i-1,:,:] + p.Δt*k3,i*p.Δt,cache)
        for j = 1:size(ℚ,2)
            for k = 1:size(ℚ,3)
                ℚ[i,j,k] = ℚ[i-1,j,k] + p.Δt/6 * (k1[1,j,k]+2*k2[1,j,k]+2*k3[1,j,k]+k4[1,j,k])
            end
        end
        postProcess!(p, ℚ, i, cache)
    end
    ℚ
end

postProcess(args...) = nothing

function postProcess!(p::DG, ℚ, i, cache)
    # post-processing
    if isa(p,DG1D) && p.useMUSCL
        @assert p.p == 1 && isa(p,DG1D) # For now
        for c in p.mesh
            Qk = ℚ[i,c.coord,:]
            ql = evaluate(Qk,basis,coord(:l))
            qr = evaluate(Qk,basis,coord(:r))

            h = (c.right-c.left)

            slope = (qr-ql)/h
            qkbar = (ql+qr)/2

            s = sign(slope)
            minabsa = abs(slope)

            for face in (:l,:r)
                if has_neighbor(p.mesh,c, face)
                    N = neighbor(p.mesh, c, face)
                    Qk = ℚ[i,N.coord,:]
                    ql = evaluate(Qk,basis,coord(:l))
                    qr = evaluate(Qk,basis,coord(:r))
                    qNbar = (ql+qr)/2

                    eeslope = n⁻(face)*(qNbar-qkbar)/h

                    minabsa = min(minabsa,abs(eeslope))

                    s += sign(eeslope)
                end
            end

            s/=3
            # Hack in knowledge about WaveState - will fix after
            # assignment is due

            local q1c, q2c
            if abs(s).q1 != 1
                q1c = [qkbar.q1 for i=1:(p.p+1)]
            else
                slope = s*minabsa
                q1c = [(qkbar+(node)*(h/2)*slope).q1 for node in cache.nodes]
            end

            if abs(s).q2 != 1
                q2c = [qkbar.q2 for i=1:(p.p+1)]
            else
                slope = s*minabsa
                q2c = [(qkbar+(node)*(h/2)*slope).q2 for node in cache.nodes]
            end

            #@show (s, q1c, q2c)

            cs = typeof(s)[typeof(s)(q1c[i],q2c[i]) for i = 1:length(q1c)]
            ℚ[i,c.coord,:] = cs
        end
    end
end

# Nodal interpolation
interpolate(f, nodes) = map(f,nodes)

function solve(p::DG)
    if p.mesh === nothing
        if isa(p,DG1D)
            p.mesh = generateMesh(Mesh1D,p.a,p.b,p.K; periodic=false)
        else
            error("Can only generate 1D uniform meshes. Please specify a mesh!")
        end
    end

    cache = generateMatrices(p)

    # Save coefficients for all timesteps (nt)
    # for all K cells
    if isa(p.q₀,Vector)
        T = eltype(p.q₀)
        ℚ = Array(T,p.nt,numcells(p),nbf(p))
        ℚ[1,:] = p.q₀ #[zero(Coeffs{p.p+1}) for i = 1:p.K]
    elseif isa(p.q₀,Function)
        T = typeof(p.q₀(zerop(p)))
        ℚ = Array(T,p.nt,numcells(p),nbf(p))
        # Use nodal interpolation
        ns = nodes(p)
        for c in p.mesh
            ℚ[1,cid(c),:] = interpolate(p.q₀,map(n->𝜒⁻¹(p.mesh,c,n),ns))
        end
    elseif isa(p.q₀, Number)
        # Uniform initial conditions
        T = typeof(p.q₀)
        ℚ = Array(T,p.nt,numcells(p),nbf(p))

        ℚ[1,:,:] = p.q₀*ones(1,numcells(p),nbf(p))
    else
        @assert "Unknown option"
    end
    @assert eltype(ℚ) != Any

    _solve(p,ℚ; cache=cache)
end


text:
module LazyStrings

# Write your package code here.
export StringWrapper
include("StringWrappers.jl")

export ReversedStrings
include("ReversedStrings.jl")

export CharMappedString
"""
    CharMappedString(f::Function,x) <: AbstractString

String implementation lazily transforming characters.
Used for parsing with [`MappedSequenceParser`](@ref).
"""
struct CharMappedString{S<:AbstractString,M<:Function} <: StringWrapper
    x::S
    f::M
    function CharMappedString(f::Function,x::AbstractString)
        new{typeof(x),typeof(f)}(x,f)
    end
end

export lmap
lmap(f::Function,x::AbstractString) =
    CharMappedString(f::Function,x::AbstractString)

tuple_pos(i) = i[1]
tuple_state(i) = i[2]

@inline Base.@propagate_inbounds Base.getindex(x::CharMappedString,i::Integer) =
    x.f(getindex(x.x,i))
@inline Base.@propagate_inbounds Base.iterate(x::CharMappedString{<:AbstractString}) =
    let i=iterate(x.x)
        i===nothing && return nothing
        x.f(tuple_pos(i)), tuple_state(i)
    end
@inline Base.@propagate_inbounds Base.iterate(x::CharMappedString{<:AbstractString},i::Integer) =
    let j=iterate(x.x,i)
        j===nothing && return nothing
        x.f(tuple_pos(j)), tuple_state(j)
    end
@inline Base.@propagate_inbounds Base.SubString(x::CharMappedString,start::Int,stop::Int) =
    lmap(x.f, SubString(x.x,start,stop))

end


text:
function gradient(
    models::Array{<:UQModel}, inputs::Array{<:UQInput}, x::DataFrame, output::Symbol
)
    samples = copy(x)

    random_names = names(filter(i -> isa(i, RandomUQInput), inputs))

    function f(x)
        samples = sample(inputs, size(x, 1))
        samples[:, random_names] .= x

        evaluate!(models, samples)

        return samples[:, output][1]
    end

    reference = Array{Float64,2}(samples[:, random_names])

    g = grad(forward_fdm(2, 1), f, reference)[1]

    return (; zip(random_names, g)...)
end

function gradient_in_standard_normal_space(
    models::Array{<:UQModel}, inputs::Array{<:UQInput}, x::DataFrame, output::Symbol
)
    samples = copy(x)

    random_names = names(filter(i -> isa(i, RandomUQInput), inputs))
    to_standard_normal_space!(inputs, samples)

    function f(x)
        samples[:, random_names] .= x
        to_physical_space!(inputs, samples)

        evaluate!(models, samples)

        return samples[:, output][1]
    end

    reference = Array{Float64,2}(samples[:, random_names])

    g = grad(forward_fdm(2, 1), f, reference)[1]

    return (; zip(random_names, g)...)
end

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
using SeisIO, SeisNoise, JLD2, PlotlyJS, StatsBase, Sockets, ORCA, Statistics

include("../pairing.jl")
include("../reference.jl")
include("../stacking.jl")

function plot_seismograms(finame::String, stn::String; norm_factor=nothing, sparse::Int64=1, foname::String="", show::Bool=true, timeslice::Array{Int64,1}=[1, -1])
    f = jldopen(finame)
    timestamplist = f["info/DLtimestamplist"]

    if timeslice[2]==-1
        timeslice = [1, length(timestamplist)]
    end

    # set up plot
    layout = Layout(width=1200, height=2000,
                    xaxis_title="Time [hours]",
                    yaxis_title="Days since $(timestamplist[1])",
                    font=attr(size=18),
                    showlegend=false,
                    title="$stn")
    p = PlotlyJS.plot([NaN], layout)

    for (titer, time) in enumerate(timestamplist[timeslice[1]:timeslice[2]])
        ts = f["$time/$stn"]
        if ts.misc["dlerror"] == 1 continue end
        # extract SeisChannel if necessary
        if typeof(ts)==SeisData ts=ts[1] end

        taxis = (collect(0:length(ts.x)-1) ./ ts.fs) ./ 3600

        if norm_factor == nothing
            norm = maximum(ts.x)
        else
            norm = norm_factor
        end

        trace = scatter(;x=taxis[1:sparse:end], y=((ts.x ./ (norm)) .+ titer .- 1)[1:sparse:end], mode="lines", line_color="rgb(0,0,0)", name=time)
        addtraces!(p, trace)
    end
    deletetraces!(p, 1)
    if foname !== ""
        PlotlyJS.savefig(p, foname)
    end
    if show
        display(p)
        readline()
    end

    close(f)
end

function plot_seismograms(finame::String; basefoname::String="", show::Bool=false, sparse::Int64=1)
    f=jldopen(finame)
    stlist = f["info/stationlist"]
    for stn in stlist
        println("Plotting $stn")
        plot_seismograms(finame, stn, foname="$basefoname.$stn.pdf", show=show, sparse=sparse)
    end
end

function plot_xcorrs(basefiname::String, stn1::String, stn2::String; reference::Union{Bool, String}=false, type::String="heatmap", foname::String="", show::Bool=true, filter::Union{Bool, Array{Float64,1}}=false, phase_smoothing::Float64=0., timeslice::Array{Int64,1}=[1,-1], stack::String="selective", threshold::Float64=0.0, metric::String="cc", slice::Union{Bool, Float64, Array{Float64,1}}=false)
    f = jldopen(basefiname*".jld2")
    timestamplist = f["info/timestamplist"]

    # set up plot
    layout = Layout(width=1200, height=800,
                    xaxis_title="Lags [s]",
                    yaxis_title="Days since $(timestamplist[1])",
                    font=attr(size=18),
                    showlegend=false,
                    title="$stn1.$stn2")
    p = PlotlyJS.plot([NaN], layout)

    stnpair = stn1*"."*stn2
    stnpairrev = stn2*"."*stn1

    if type=="heatmap" xcorr_heat = []; tscounter=0 end

    if stack=="selective"
        rmList = Array{Float64,1}()
    end

    if timeslice[2]==-1
        timeslice=[1,length(timestamplist)]
    end

    for (titer, time) in enumerate(timestamplist[timeslice[1]:timeslice[2]])
        f_cur = jldopen(basefiname*".$time.jld2")

        # load cross-correlation or its reverse
        if stnpair ∈ keys(f_cur[time])
            xcorr = f_cur["$time/$stnpair"]
            if filter!=false
                xcorr.corr = bandpass(xcorr.corr, filter[1], filter[2], xcorr.fs, corners=4, zerophase=false)
            end
            nWins = length(xcorr.corr[1,:])
            if stack=="selective"
                if reference!=false
                    f_ref=jldopen(reference)
                    ref = f_ref[stnpair]

                    if filter!=false
                        ref.corr = bandpass(ref.corr, filter[1], filter[2], xcorr.fs, corners=4, zerophase=false)
                    end
                    xcorr, ccList = selective_stacking(xcorr, ref, threshold=threshold, metric=metric, filter=filter)
                    close(f_ref)
                else
                    xcorr, ccList = selective_stacking(xcorr, threshold=threshold, metric=metric, filter=filter)
                end
                nRem = length(findall(x->(x<threshold), ccList))
                push!(rmList, nRem / nWins)
            else
                stack!(xcorr, allstack=true, phase_smoothing=phase_smoothing)
            end

        elseif stnpairrev ∈ keys(f_cur[time])
            xcorr = f_cur["$time/$stnpairrev"]
            if filter!=false
                xcorr.corr = bandpass(xcorr.corr, filter[1], filter[2], xcorr.fs, corners=4, zerophase=false)
            end
            nWins = length(xcorr.corr[1,:])
            if stack=="selective"
                if reference!=false
                    f_ref=jldopen(reference)
                    ref = f_ref[stnpairrev]
                    if filter!=false
                        ref.corr = bandpass(ref.corr, filter[1], filter[2], xcorr.fs, corners=4, zerophase=false)
                    end
                    xcorr, ccList = selective_stacking(xcorr, ref, threshold=threshold, metric=metric, filter=filter)
                    close(f_ref)
                else
                    xcorr, ccList = selective_stacking(xcorr, threshold=threshold, metric=metric, filter=filter)
                end
                nRem = length(findall(x->(x<threshold), ccList))
                push!(rmList, nRem / nWins)
            else
                stack!(xcorr, allstack=true, phase_smoothing=phase_smoothing)
            end
            xcorr.corr = reverse(xcorr.corr, dims=1)

        else
            xcorr = CorrData()
            xcorr.corr = zeros(4001,1)
        end

        # make lags visible outside of loop
        global lags = -xcorr.maxlag:1/xcorr.fs:xcorr.maxlag

        norm_factor = maximum(abs.(xcorr.corr[:]))

        try
            if type=="wiggles"
                # add trace for wiggle plot
                trace = PlotlyJS.scatter(;x=lags, y=xcorr.corr[:, 1] ./ (2 * norm_factor) .+ titer .- 1, mode="lines", line_color="black", name=time)
                addtraces!(p, trace)

            elseif type=="heatmap"
                # append to array for heat map. We cannot plot before it is fully filled
                append!(xcorr_heat, (xcorr.corr[:, 1]./ (norm_factor)))
                tscounter+=1
            end
        catch y
            println(y)
            println("Plotting failed")
        end
        close(f_cur)
    end

    if stack=="selective"
        p1=PlotlyJS.plot(rmList.*100)
        display(p1)
        readline()
    end

    if type=="heatmap"
        # dont plot if no cross-correlations were found. This would give div by 0
        if tscounter==0 println("No cross-correlations to plot. Exiting."); exit() end
        xcorr_heat = reshape(xcorr_heat, convert(Int64, length(xcorr_heat)/tscounter), tscounter)
        trace = PlotlyJS.heatmap(x=lags, z=xcorr_heat)
        addtraces!(p, trace)
    end
    deletetraces!(p, 1)

    if foname !== ""
        PlotlyJS.savefig(p, foname)
    end
    if show
        display(p)
        readline()
    end

    close(f)
end

function plot_xcorrs(basefiname::String, basefoname::String; show::Bool=false, type::String="heatmap")
    f=jldopen(basefiname*".jld2")
    stlist = f["info/stationlist"]
    stniter=1
    for stn1 in stlist
        for stn2 in stlist[stniter:end]
            println("Processing $stn1--$stn2")
            plot_xcorrs(basefiname, stn1, stn2; type=type, foname=basefoname*".$stn1.$stn2.png", show=show)
        end
        stniter+=1
    end
end

function plot_corrcoeff(basefiname::String, reference::String)
    # read base file for time and station information
    f=jldopen(basefiname*".jld2")
    stlist = f["info/stationlist"]
    tslist = f["info/timestamplist"]
    close(f)
    # reference cross-correlations
    f_ref = jldopen(reference)

    # empty dictionary containing the correlation coefficients for each station pair through time
    ccDict = Dict{String, Array{Float64,1}}()

    # iterate over station pairs present at the current time step
    for t in tslist
        println("Processing $t")
        # open file containing cross-correlations for the current time stamp
        f_cur = jldopen(basefiname*".$t.jld2")
        for stnpair in keys(f_cur[t])
            # read cross-correlations for current station pair
            xcorr = f_cur["$t/$stnpair"]
            # read corresponding reference
            ref = f_ref[stnpair]
            # compute correlation-coefficient
            ccList = get_cc(xcorr, ref)
            # append to the dictionary
            if haskey(ccDict, stnpair)
                append!(ccDict[stnpair], ccList)
            else
                ccDict[stnpair] = ccList
            end
        end
        close(f_cur) # current timestamp cross-correlations
    end
    close(f_ref) # reference cross-correlations
    return ccDict
end

function plot_corrcoeff(basefiname::String, reference::String, stnpair::String)
    # read base file for time and station information
    f=jldopen(basefiname*".jld2")
    tslist = f["info/timestamplist"]
    close(f)
    # reference cross-correlations
    f_ref = jldopen(reference)

    # empty dictionary containing the correlation coefficients for each station pair through time
    ccDict = Dict{String, Array{Float64,1}}()

    # empty dictionary containing the maximum amplitude of each cross-correlation window
    maxAmp = Dict{String, Array{Float32,1}}()

    # iterate over station pairs present at the current time step
    for t in tslist
        # open file containing cross-correlations for the current time stamp
        f_cur = jldopen(basefiname*".$t.jld2")
        # read cross-correlations for current station pair
        xcorr = try f_cur["$t/$stnpair"] catch; continue end
        dist = xcorr.misc["dist"]
        # read corresponding reference
        ref = f_ref[stnpair]
        # compute correlation-coefficient
        ccList = get_cc(xcorr, ref)
        # extract maximum amplitude for each windowed cross-correlation
        maxWinAmp = maximum(xcorr.corr, dims=1)[:]

        # append to the dictionary
        if haskey(ccDict, stnpair)
            append!(ccDict[stnpair], ccList)
        else
            ccDict[stnpair] = ccList
        end
        if haskey(maxAmp, stnpair)
            append!(maxAmp[stnpair], maxWinAmp)
        else
            maxAmp[stnpair] = maxWinAmp
        end
        close(f_cur) # current timestamp cross-correlations
    end
    close(f_ref) # reference cross-correlations
    return ccDict, maxAmp, dist
end

function plot_reference(finame::String, stn1::String, stn2::String)
    f = jldopen(finame)
    stname = stn1*"."*stn2
    data = f[stname]

    lags = -data.maxlag:1/data.fs:data.maxlag

    # set up plot
    layout = Layout(width=1200, height=800,
                    xaxis_title="Lags [s]",
                    font =attr(size=18),
                    showlegend=true,
                    title="$stname",
                    paper_bgcolor="rgba(0,0,0,0)",
                    plot_bgcolor="rgba(0,0,0,0)")

    p = PlotlyJS.plot([NaN], layout)
    trace1 = PlotlyJS.scatter(;x=lags, y=data.corr[:, 1], mode="lines", line_color="rgb(0,0,0)")
    addtraces!(p, trace1)
    deletetraces!(p, 1)
    display(p)
    readline()
end

function plot_convergence(finame::String, stn1::String, stn2::String; maxlag::Float64=100.0, fs::Float64=20.0)
    f = jldopen(finame)
    stname = stn1*"."*stn2
    data = f[stname]

    lags = -maxlag:1/fs:maxlag

    # set up plot
    layout = Layout(width=1200, height=800,
                    xaxis_title="# of daily stacked cross-correlations",
                    font=attr(size=18),
                    showlegend=true,
                    title="$stname RMS convergence")

    p = PlotlyJS.plot([NaN], layout)

    trace1 = PlotlyJS.scatter(;x=collect(1:length(data))./46, y=data, mode="lines", line_color="rgb(0,0,0)")
    addtraces!(p, trace1)
    deletetraces!(p, 1)
    display(p)
    readline()
end

function plot_convergence(finame::String; sparse::Int64=1)
    f=jldopen(finame)
    lags = -100.0:0.05:100.0

    # set up plot
    layout = Layout(width=1200, height=800,
                    xaxis_title="# of daily stacked cross-correlations",
                    font=attr(size=18),
                    showlegend=true,
                    title="RMS Convergence")
    p = PlotlyJS.plot([NaN], layout)
    for stpair in keys(f)[1:200]
        (net1, stn1, comp1, net2, stn2, comp2) = filter!(x->x≠"", split(stpair, "."))
        if comp1==comp2 && stn1!=stn2 && comp1=="BP3"
            color="rgba(150,0,0,1.0)"
        elseif stn1==stn2
            continue
        else
            color="rgba(0,0,0,0.1)"
        end
        data = f[stpair]
        trace = PlotlyJS.scatter(;x=(collect(1:length(data))[1:sparse:end])./46, y=data[1:sparse:end], mode="lines", line_color=color, name=stpair)
        addtraces!(p, trace)
    end
    deletetraces!(p, 1)
    display(p)
    readline()
end

function count_zeros(ref::Union{Array{Float32,1}, Array{Float64,1}}, data::Union{Array{Float32,1}, Array{Float64,1}}; percent::Bool=true)
    # get number of zeros in two given arrays, and compute the difference.
    # in the expected use case, this finds the amount of data zeroed by earthquake removal
    nZeros = count(x->x==0, data)
    nZerosRef = count(x->x==0, ref)
    rem = nZeros-nZerosRef
    # convert number of removed zeros to percent of total data length
    if percent
        rem /= length(ref)
        rem *= 100
    end
    return rem
end

function plot_removal(orig_finame::String, finame::String; kurtlist::Array{Float64,1}, staltalist::Array{Float64,1})
    # open original input file
    f_orig = jldopen(orig_finame*".jld2")
    # time stamps and stations
    tslist = f_orig["info/DLtimestamplist"]
    stlist = f_orig["info/stationlist"]

    # initialize 2D array to store number of removed data points
    removed = zeros(length(staltalist), length(kurtlist), length(stlist))
    # initialize 3D array to store peak amplitude
    amps = zeros(length(staltalist), length(kurtlist), length(tslist), length(stlist))
    # initialize 3D array to store average amplitude
    avgs = zeros(length(staltalist), length(kurtlist), length(tslist), length(stlist))

    for (i, stalta) in enumerate(staltalist)
        for (j, kurt) in enumerate(kurtlist)
            println("Processing stalta: $stalta, kurtosis: $kurt")
            # open file for each stalta threshold and kurtosis threshold
            f=jldopen(finame*"_$(stalta)_$(kurt).jld2")
            # count number of zeros for all time stamps
            for (t,ts) in enumerate(tslist)
                for (s, stn) in enumerate(stlist)
                    orig_data = f_orig["$ts/$stn"].x
                    data = f["$ts/$stn"].x
                    removed[i,j,s] += count_zeros(orig_data, data)
                    amps[i,j,t,s] = maximum(abs.(data))
                    avgs[i,j,t,s] = mean(data)
                end
            end
            close(f)

            # average percent removal over days
            removed[i,j,:] ./= length(tslist)
        end
    end
    close(f_orig)

    for (s, stn) in enumerate(stlist)
        # set up plot
        layout = Layout(width=1200, height=800,
                        xaxis_title="STA/LTA Threshold",
                        yaxis_title="Kurtosis Threshold",
                        font=attr(size=18),
                        showlegend=true,
                        title="$stn % removal")

        p = PlotlyJS.plot([NaN], layout)
        trace = PlotlyJS.heatmap(x=staltalist, y=kurtlist, z=removed[:,:,s])
        addtraces!(p, trace)
        deletetraces!(p, 1)
        display(p)
        readline()
    end
end

function plot_dvv(finame::String, stn1::String, stn2::String)
    f=jldopen(finame)
    dvv = f["$stn1.$stn2"]
    tslist = f["info/timestamplist"]

    layout = Layout(width=1200, height=800,
                    xaxis_title="Days since $(tslist[1])",
                    yaxis_title="dv/v (%)",
                    font=attr(size=18),
                    showlegend=false,
                    title="$stn1--$stn2")

    p = PlotlyJS.plot([NaN], layout)
    trace = PlotlyJS.scatter(;x=collect(1:length(tslist)), y=dvv, mode="lines+markers", line_color="rgba(0,0,0,1)")
    addtraces!(p, trace)
    deletetraces!(p, 1)
    display(p)
    readline()
end

function plot_dvv(finame::String, stack=false)
    f=jldopen(finame)
    tslist = f["info/timestamplist"]

    layout = Layout(width=1200, height=800,
                    xaxis_title="Days since $(tslist[1])",
                    yaxis_title="dv/v (%)",
                    font=attr(size=18),
                    showlegend=false)
    p = PlotlyJS.plot([NaN], layout)

    if stack
        dvvfull = zeros(31)
        count = 0
    end

    for stpair in keys(f)[1:780]
        (net1, stn1, comp1, net2, stn2, comp2) = filter!(x->x≠"", split(stpair, "."))
        if stpair == "BP.CCRB..BP1.BP.SMNB..BP1"
            color="rgba(150,0,0,1)"
        elseif stn1==stn2
            continue
        elseif stn1!=stn2 && comp1!=comp2
            continue
        else
            color = "rgba(0,0,0,0.2)"
        end

        if stack
            dvv = f[stpair]
            dvvfull .+= dvv
            count += 1
        else
            trace = PlotlyJS.scatter(;x=collect(1:length(tslist)), y=dvv, mode="lines+markers", line_color=color)
            addtraces!(p,trace)
        end
    end

    if stack
        trace = PlotlyJS.scatter(;x=collect(1:length(tslist)), y=dvvfull./count, mode="lines+markers", line_color="rgba(0,0,0,1)")
        addtraces!(p, trace)
    end

    close(f)
    deletetraces!(p,1)
    display(p)
    readline()
end

function plot_dvv(basefiname::String, methods::Array{String,1}, percent::Array{Bool,1})
    f1=jldopen(basefiname*"_$(methods[1]).jld2")
    tslist = f1["info/timestamplist"]
    close(f1)
    layout = Layout(width=1200, height=800,
                    xaxis_title="Days since $(tslist[1])",
                    yaxis_title="dv/v (%)",
                    yaxis_range=[-0.1, 0.1],
                    font=attr(size=18),
                    showlegend=true)
    p = PlotlyJS.plot([NaN], layout)

    for (i,method) in enumerate(methods)
        f=jldopen(basefiname*"_$(methods[i]).jld2")

        dvvfull = zeros(31)
        count=0

        for stpair in keys(f)[1:780]
            (net1, stn1, comp1, net2, stn2, comp2) = filter!(x->x≠"", split(stpair, "."))
            if stpair == "BP.CCRB..BP1.BP.SMNB..BP1"
                color="rgba(150,0,0,1)"
            elseif stn1==stn2
                continue
            elseif stn1!=stn2 && comp1!=comp2
                continue
            else
                color = "rgba(0,0,0,0.2)"
            end

            dvv = f[stpair]
            dvvfull .+= dvv
            count += 1
        end

        if percent[i]!=true dvvfull .*= 100 end
        trace = PlotlyJS.scatter(;x=collect(1:length(tslist)), y=dvvfull./count, mode="lines+markers", line_color="rgba(0,0,0,1)", name="$method")
        addtraces!(p, trace)
        close(f)
    end
    deletetraces!(p, 1)
    display(p)
    readline()
end

max_stars_repo_name: kura-okubo/SeisXcorrelation

max_issues_repo_name: kura-okubo/SeisXcorrelation

max_forks_repo_name: kura-okubo/SeisXcorrelation