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

using Markdown
using InteractiveUtils

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

# ╔═╡ 0c45e9aa-e0ba-11eb-0f56-af8192849f70
begin
  using Pkg
  Pkg.add([
    "PlutoUI",
    "Plots",
  ])
  using PlutoUI
  using Plots
  #plotly()
  #gr()
  using Random
end

# ╔═╡ e840f8ca-4ed6-4839-baeb-312cc457e5be
PlutoUI.TableOfContents(aside=true)

# ╔═╡ 293031b4-e0bb-11eb-1f5b-9da7ddf56b6c
md"""
# Motivations
When we study undergraduate maths, in particular, linear algebra, the symbol ``\mathscr{L}(E, F)`` is sometimes employed to mean the space of all linear transformations from ``E`` to ``F``. But then, when we study a little bit higher, our books and professors start to use ``\mathscr{L}(E, F)`` to mean the space of all linear **_continuous_** transformations from ``E`` to ``F``. The topic we want to discuss here is
> Why sometimes there is continuity and somtimes not in the notation of ``\mathscr{L}(E, F)``?

Let's first read the following statement.
> To define a norm on the space ``%\mathscr{L}(E, F)`` of linear transformations from
> ``E`` to ``F`` (both vector sapces on ``\mathbb{R}`` or both on ``\mathbb{C}``),
> one way is to equipe the ``E`` and ``F`` with norms, and the linear transformations with continuity.
>
> In other words, the space of linear continuous transformations
> ``%\mathscr{L}(E, F)`` admits a norm.

We will understand why this is true later. For the moment being, I would like us to
think together about what drives us to define a norm for such entities as linear transformations.

According to what I have learned, one such motivation comes from calculus:$(HTML("
"))
When we study calculus, we all have to learn the second derivative.
Let us recall how the second derivative is introduced in an abstract setting.

Let ``E, F`` be two normed vector spaces (We need norms to be able to speak of
differentiability), and let ``f: U \subset E \to F`` be a function, where ``U`` is open in ``E\,.`` (The term "_open_" simply means that, for every point ``a`` in ``U``, we are assured that every neighboring point of ``a``, provided close enough, lies also in ``U``.) Then, by definition, ``f`` is differentiable at a point ``a \in U`` if there exists a linear transformation ``f'(a): E \to F`` s.t.
```math
\newcommand{\norm}[1]{\lVert{#1}\rVert}
f(a + h) = f(a) + f'(a)(h) + o(\norm{h})\,,
```
where the `` o(\norm{h})`` simply means
```math
\lim_{\norm{h}\to 0} \frac{\norm{f(a + h) - f(a) - f'(a)(h)}}{\norm{h}} = 0\,.
```

> **(?)** The derivative ``f'(a)``, does it **_suffice_** to be **_merely linear_**?
> Or it needs to be **_continuous_** as well?

Note that, if we are to further discuss whether ``f`` is twice differentiable at ``a``, we
need that

- ``f'`` be defined on an open set containing ``a`` denoted ``N_{a}\,``.
  - ``f': N_{a} \to G``, where ``G`` denotes the space of linear transformations from ``E`` to ``F``.
- ``G`` **_have a norm_** in order to discuss the differentiability of ``f'``.

So, we see one necessity -- We need a norm for linear transfomrations from ``E`` to ``F``.
"""

# ╔═╡ b65e4364-e0cf-11eb-0d0e-2f5418a1286a
md"""
# A Norm on Linear Transformations
How do we define a norm on linear transformations ``E \to F\,``? I guess we can find
one by trial and error.

Let ``T: E \to F`` be a linear transformation. Then
for all coefficient ``a`` and all vector ``x \in E``, we have
```math
\newcommand{\norm}[1]{\lVert{#1}\rVert}
\newcommand{\abs}[1]{\lvert{#1}\rvert}
\norm{T(ax)} = \abs{a}\;\norm{T(x)}\,,
```

In particular, if there exists some ``x_{0} \in E`` s.t. ``T(x_{0}) \ne 0``, then we see that
```math
\begin{align}
  \lim_{a \to \infty} \norm{aT(x_{0})} &= \infty\,, \\
  \lim_{a \to -\infty} \norm{aT(x_{0})} &= \infty\,, \\
  \lim_{a \to 0} \norm{aT(x_{0})} &= 0\,.
\end{align}
```

If no such ``x_{0}`` exists, then, uninterestingly, ``T = 0``, the zero function.

Coined in English, the above just tells us that
> In every direction ``x`` in which ``T`` does not zero out, the norm ``\norm{T(ax)}`` tends to ``\infty`` as it goes to the extremities.

So one idea could be to consider a representative for each and every direction, and then try to produce a norm from these representatives, like a progressive election:

- either ``\norm{T} := \sup_{\norm{x}=1} \norm{T(x)}``
- or ``\quad\;\,\norm{T} := \inf_{\norm{x}=1} \,\norm{T(x)}``

The definition via ``\inf``  has some drawbacks. For example,

- if ``\norm{T(x_{0})} = 0`` for some ``x_{0} \ne 0`` then ``\norm{T(\frac{1}{\norm{x_{0}}} x_{0})} = \frac{1}{\norm{x_{0}}} \norm{T(x_{0})} = 0\,.``
  - Since, ``\frac{1}{\norm{x_{0}}} x_{0}`` is a unit vector, this implies that
    ``\inf_{\norm{x}=1} \norm{T(x)} = 0\,.``
  - Had the ``\inf`` definition been a valid norm, we would have
    ``\norm{T} = 0``, which implies ``T = 0``, the zero function.
    This is the drawback I mentioned: when ``T`` zeros out in some direction,
    ``T`` has to zero out in all directions in order for the ``\inf`` definition
    to become a valid norm. This is a strict requirement. Let alone saying
    requiring so much has not readily made the ``\inf`` definition a real norm.
    It was just a necessary condition.

Let's see if ``\inf``'s counterpart, the ``\sup`` definition, can make it until the end. And then we will come back to see if it is still possible to have a norm defined via ``\inf\,.``
"""

# ╔═╡ 6e593c67-1d05-41ae-bee6-6b63605b28ef
md"""
## ``\norm{T} := \sup_{\norm{x}=1} \norm{T(x)}``
First of all, to be able to define the norm like this, we must require that the set
```math
\left\{\,T(x) \,\middle\vert\; \norm{x}=1\,\right\}
```
be bounded. This is not automatic. Indeed, as we shall see, there exist linear functions that do not satisfy this.

We show that the norm ``\norm{\cdot}`` defined on the set

```math
V := \left\{T \,\middle|\; T: E \to F \;\text{linear} \quad\text{and}\quad \sup_{\norm{x}=1} \norm{T(x)} < \infty  \right\}
```
is a norm.

We'd better first show that ``V`` is a vector space. Since the demonstration should be quite straightforward, we only pick a few points to prove.

- ``T, S \in V \implies T + S \in V``
  - ```math
    \sup_{\norm{x}=1} \norm{(T+S)(x)} \le \sup_{\norm{x}=1} \left(\norm{T(x)} +
    \norm{S(x)}\right)
    \le \sup_{\norm{x}=1} \norm{T(x)} + \sup_{\norm{x}=1} \norm{S(x)}
    ```
- ``a \;\text{coefficient}, T \in S \implies aT \in S``
  - ```math
    \sup_{\norm{x}=1} \norm{(aT)(x)} = \sup_{\norm{x}=1} \abs{a}\norm{T(x)}
    = \abs{a}\sup_{\norm{x}=1} \norm{T(x)}
    ```
  - Note that the above is stronger than an inequality. We prove it stronger because later on we can take this result to establish one condition that a norm must satisfy.

It is also quite straightforward to prove that ``\norm{\cdot}`` is a norm on ``V``.

- ``\norm{T} \ge 0 \quad\forall\; T \in V``
- ``\norm{T} = 0 \iff T = 0``
- ``\norm{aT} = \abs{a} \norm{T} \quad\forall\; a, T``
- ``\norm{T+S} \le \norm{T} + \norm{S} \quad\forall\; T, S \in V``

In conclusion,
> we've found a norm defined on ``V``, a subspace of the space of all linear functions from ``E`` to ``F``.

Great! $(HTML("
"))
Remember our example on derivatives? If we ask the derivative ``f'(a)`` to lie in ``V``, then there is a norm on every element of ``V`` and it is valid for us
to discuss ``f''``.

> But what does this have to do with continuity?

It turns out that the demand
```math
\sup_{\norm{x}=1} \norm{T(x)} < \infty
```
coupled with the fact that ``T`` is linear constitutes continuity.

Indeed, assuming that ``T`` is linear, then it is well-known that
```math
\sup_{\norm{x}=1} \norm{T(x)} < \infty \iff
T \;\text{continuous at}\; 0 \iff
T \;\text{continuous on}\; E\,.
```

Those who don't believe may benefit from proving these implications by themselves.

Note that, in general, ``V`` is a proper subset of the set of all linear functions from ``E`` to ``F``. That is to say, there exist discontinuous linear transformations. Cf. [Appendix B](## B. Examples of Discontinuous Linear Transformations)
"""

# ╔═╡ 1d65c197-3da7-49b8-a97b-95d69607e269
md"""
## Visualization
Let's use our computer to help us visualize the situation better. More precisely, let's denote ``\Gamma = \left\{\,x \in \mathbb{R}^{n} \,\middle\vert\; \norm{x}=1\,\right\}``, where ``n = 2 \;\text{or}\; 3``, and try to plot the set ``T(\Gamma)`` for the following linear transformations:

- ``T: \mathbb{R}^{2} \to \mathbb{R}^{2}``
- ``T: \mathbb{R}^{3} \to \mathbb{R}^{3}``
- ``T: \mathbb{R}^{3} \to \mathbb{R}^{2}``

**Rmk.** One may adjust the values inside the next cell, and the following plots will change accordingly. (If one changes `Tu, Tv, Tw`, one'd better change `xlim, ylim` in the plot functions as well, so that one does not lose sight by focusing on the wrong region.)
"""

# ╔═╡ e3bee7b9-bb79-4473-8665-30066aeb1c53
begin
  Tu = [3. 4.]'
  Tv = [-1. 0.]'
  Tw = [1. -2.]'
  n_angles = 30
  md"""
  `seed` = $(@bind seed Slider(0:42, show_value=true, default=15))
  """
end

# ╔═╡ 50c79c60-f332-4c40-a32f-2de732842224
function T(coeff, basis)
  return basis * coeff
end

# ╔═╡ c6828156-9a1f-4dc5-8150-0d99972db825
begin
  ϵ = 1e-9
  xy = zeros((2, n_angles))
  for (col, θ) in enumerate(range(0+ϵ, 2π-ϵ; length=n_angles))
    xy[1:end, col] = [cos(θ), sin(θ)]
  end
  pts22 = T(xy, [Tu Tv])
end

# ╔═╡ 3836840d-9a62-4b67-9d56-5f4347b55d75
begin
  xyz = zeros((3, n_angles^2))
  counter = 1
  for ϕ in range(0+ϵ, π-ϵ; length=n_angles), θ in range(0+ϵ, 2π-ϵ; length=n_angles)
    xyz[1:end, counter] = [sin(ϕ)*cos(θ), sin(ϕ)*sin(θ), cos(ϕ)]
    global counter += 1
  end
  Random.seed!(seed)
  pts33 = T(xyz, float(rand(-10:10, (3,3))))
end

# ╔═╡ 6595c9d5-1005-498c-8aad-0ad07e90b4e0
pts32 = T(xyz, [Tu Tv Tw])

# ╔═╡ c5439555-3cd8-44a7-a4df-78810b381efd
let
  # stop
  sct = Plots.scatter(
    1,
    xlim=(-7, 7),
    ylim=(-7, 7),
    aspect_ratio=:equal,
    background_color=:black,
    title="T: R² to R²",
    #marker=2,
    label=false,
  )
  @gif for col = 1:size(pts22, 2)
    pt = pts22[1:end, col]
    push!(sct, pt[1], pt[2])
  end every 1
end

# ╔═╡ 01c4b8f1-709b-4a7f-ab5d-00e51949fd0b
let
  plt = Plots.plot3d(
    1,
    xlim=(-15, 15),
    ylim=(-20, 20),
    zlim=(-10, 10),
    aspect_ratio=:equal,
    background_color=:black,
    title="T: R³ to R³",
    label=false,
  )
  @gif for col in 1:size(pts33, 2)
    push!(plt, pts33[1, col], pts33[2, col], pts33[3, col])
  end every 5
end

# ╔═╡ 38dc587f-b10b-4c72-9d29-5487d2bfddba
let
  sct = Plots.scatter(
    1,
    xlim=(-7, 7),
    ylim=(-7, 7),
    aspect_ratio=:equal,
    background_color=:black,
    title="T: R³ to R²",
    marker=3,
    label=false,
  )
  @gif for col = 1:size(pts32, 2)
    pt = pts32[1:end, col]
    push!(sct, pt[1], pt[2])
  end every 3
end

# ╔═╡ e3531807-65f0-4d0d-ab97-c6f21897c3a0
md"""
As we can see, ``T(\Gamma)`` is

- **_ellipse_** in the case ``\,T: \mathbb{R}^{2} \to \mathbb{R}^{2}``
- **_ellipsoid_** in the case ``\,T: \mathbb{R}^{3} \to \mathbb{R}^{3}``
- **_something bounded_** in the case ``\,T: \mathbb{R}^{3} \to \mathbb{R}^{2}`` (Note that it is a **2D plot**, not a 3D one.)

We observed a general elliptic pattern here. In fact, there exist techniques with which we can foresee that the result is ellipsoid in the case of ``T: \mathbb{R}^{n} \to \mathbb{R}^{n}``, e.g. Singular Value Decomposition (SVD). Uninformed but interested readers should refer to textbooks or the Internet, because SVD is out of scope of this article.
"""

# ╔═╡ 5fa9dfdf-7ae3-42ac-a02c-b89afbdcad28
md"""
## ``\norm{T} := \inf_{\norm{x}=1} \norm{T(x)}``
Unlike for ``\sup``, since any norm is always bounded below by ``0``, ``\inf_{\norm{x}=1} \norm{T(x)}`` always exist for all ``T: E \to F``.
So here we don't have to restrict to a subset ``V`` of the set of all linear
transformations from ``E`` to ``F``.

In general, this ``\inf`` definition does not give us a norm,
as the following example shows. Consider
``E = F = \mathbb{R}^{2}`` and the two linear functions ``T, S: E \to F``
```math
\begin{align}
  T: e_{1} \mapsto e_{1},\; e_{2} \mapsto 0 \\
  S: e_{2} \mapsto e_{2},\; e_{1} \mapsto 0 \\
\end{align}
```
where ``\{e_{1}, e_{2}\}`` denotes the canonical basis for ``\mathbb{R}^{2}``.

Then ``T+S = I``, the identity function ``\mathbb{R}^{2} \to \mathbb{R}^{2}``. Indeed,
for all ``x \in \mathbb{R}^{2}``, we see that
```math
\begin{align}
(T+S)(x) &= (T+S)(x_{1}e_{1} + x_{2}e_{2}) \\
         &= T(x_{1}e_{1} + x_{2}e_{2}) + S(x_{1}e_{1} + x_{2}e_{2}) \\
         &= x_{1}e_{1} + x_{2}e_{2} \\
         &= x\,.
\end{align}
```

As a consequence,
```math
\inf_{\norm{x}=1} \norm{(T+S)(x)} = 1 > 0 + 0 =
\inf_{\norm{x}=1} \norm{T(x)} + \inf_{\norm{x}=1} \norm{S(x)},
```
which says that the triangular inequality simply cannot be satisfied.
"""

# ╔═╡ 1e673c5f-bd3d-4e5b-a094-e0dbee6c4f8d
md"""
## A Comment on Linear Algebra
This comment is targeted at a small portion of readers of the current article, including myself 10 years ago. More specifically, it is targeted at some undergraduates like I was a decade ago.

When we enter into college, mathematics major, the first year undergraduate mathematics conventionally teach us calculus and linear algebra.
Then we go on to study ODE, PDE, differential geometry, complex analysis, and so on, and so forth. But when asked, "what good is linear algebra?", a portion among us will not be able to give an interesting answer at the age of a college student.

One viewpoint that this article tries to convey is that
> Actually, linearity is the simplest and almost always the first thing that
> we could think of to do. Although there is _algebra_ in the name of
> _linear algebra_, linear algebra is not constrained to be used in algebra.
> Even in calculus, which is classified as an _analysis_ course,
> such fundamental entities as the derivatives are linearities, like we just saw.
>
> Linear algebra is used in a lot of places, both in various math branches
> and in practical engineering.
>
> We simplify things and we study linear algebra so that
> we know (at least) how to cope with simplified, linear things.
"""

# ╔═╡ 04895de2-1f06-4a55-9989-9fa26505a744
md"""
# References
The materials presented here are generally considered quite elementary and could probably be found in many books and on the Internet. In particular, the author first read these in the following book

- [Cours de calcul différentiel, Henri Cartan](https://www.amazon.com/Cours-calcul-diff%C3%A9rentiel-Henri-Cartan/dp/2705667024)
"""

# ╔═╡ 83ed0c76-b9cc-4a9d-8583-dd82f7447775
md"""
# Appendices
At first I thought I was going to write sth in these appendices. It turned out that

01. I am running out of time before the submission deadline of [`some1`](https://www.3blue1brown.com/blog/some1)
02. There are, unsurprisingly, excellent resources online on the topics that these appendices try to cover, and I am not confident that my own writing will render the ideas clearer than they do.

As a consequence, I will temporarily refer the interested readers to these resources.

## A. Continuity Is Automatic when ``\dim(E)`` Is Finite
That is,
> If ``T: E \to F`` is linear and ``\dim(E)`` is finite, then, ``T`` is a continuous function.

**Proof.** $(HTML("
"))
Let ``\{e_{1}, e_{2}, \ldots, e_{n}\}`` be a basis for ``E``. Then for all ``x \in E`` with ``\norm{x}_{E} = 1``, we have
```math
  \norm{T(x)}_{F} =
  %\norm{T(\sum_{i=1}^{n} x_{i}e_{i})}_{F} =
  \norm{\sum_{i=1}^{n} x_{i}T(e_{i})}_{F} \le
  \sum_{i=1}^{n} \abs{x_{i}}\, \norm{T(e_{i})}_{F} \le
  \left(\max_{i=1,2,\ldots, n} \norm{T(e_{i})}_{F}\right) \sum_{i=1}^{n} \abs{x_{i}}
```
It can be shown that the mapping ``x \mapsto \sum_{i=1}^{n} \abs{x_{i}}`` defines a norm on ``E``. Additionally, we know that all norms on a finite-dimensional vector space are **_equivalent_**, whence, continuing the foregoing inequality gives us
```math
  \norm{T(x)}_{F} \le
  M \sum_{i=1}^{n} \abs{x_{i}} \le
  M C \norm{x}_{E} =
  M C
```
for some constant ``C > 0``. (We have denoted ``M = \max_{i=1,2,\ldots, n} \norm{T(e_{i})}_{F}`` and used the fact that ``\norm{x}_{E} = 1``.)

This completes the proof that ``T`` is continuous.

## B. Examples of Discontinuous Linear Transformations
- This Wikipedea entry consists of much more than a mere example: <https://en.wikipedia.org/wiki/Discontinuous_linear_map>
"""

# ╔═╡ Cell order:
# ╟─0c45e9aa-e0ba-11eb-0f56-af8192849f70
# ╟─e840f8ca-4ed6-4839-baeb-312cc457e5be
# ╟─293031b4-e0bb-11eb-1f5b-9da7ddf56b6c
# ╟─b65e4364-e0cf-11eb-0d0e-2f5418a1286a
# ╟─6e593c67-1d05-41ae-bee6-6b63605b28ef
# ╟─1d65c197-3da7-49b8-a97b-95d69607e269
# ╠═e3bee7b9-bb79-4473-8665-30066aeb1c53
# ╠═50c79c60-f332-4c40-a32f-2de732842224
# ╟─c6828156-9a1f-4dc5-8150-0d99972db825
# ╟─6595c9d5-1005-498c-8aad-0ad07e90b4e0
# ╟─3836840d-9a62-4b67-9d56-5f4347b55d75
# ╟─c5439555-3cd8-44a7-a4df-78810b381efd
# ╟─01c4b8f1-709b-4a7f-ab5d-00e51949fd0b
# ╟─38dc587f-b10b-4c72-9d29-5487d2bfddba
# ╟─e3531807-65f0-4d0d-ab97-c6f21897c3a0
# ╟─5fa9dfdf-7ae3-42ac-a02c-b89afbdcad28
# ╟─1e673c5f-bd3d-4e5b-a094-e0dbee6c4f8d
# ╟─04895de2-1f06-4a55-9989-9fa26505a744
# ╟─83ed0c76-b9cc-4a9d-8583-dd82f7447775


text:
function F_yeast(x, parms)
  R, L, RL, G, Ga, Gbg, Gd = x
  k1, k2, k3, k4, k5, k6, k7, k8 = parms

  receptor_uregulate = k1
  receptor_dregulate = k2 * R
  ligand_binding = k3 * L * R
  ligand_degrade = k4 * RL
  G_activation = k5 * RL * G
  G_dephosphorylate = k6 * Ga
  G_rebinding = k7 * Gd * Gbg
  bound_receptor_uregulate = k8

  rates = [
    receptor_uregulate,
    receptor_dregulate,
    ligand_binding,
    ligand_degrade,
    G_activation,
    G_dephosphorylate,
    G_rebinding,
    bound_receptor_uregulate
  ]
  return rates
end

function gillespie_yeast(;
  R=50, L=2, RL=0, G=50, Ga=0, Gbg=0, Gd=0,
  k1=0.0038, k2=4e-4, k3=0.042, k4=0.01, k5=0.011, k6=0.1, k7=1.05e-3, k8=3.21)
  # initial conditions
  x0 = [R, L, RL, G, Ga, Gbg, Gd]
  # stoichiometries
  nu = [
    [1 0 0 0 0 0 0];
    [-1 0 0 0 0 0 0];
    [-1 0 1 0 0 0 0];
    [1 0 -1 0 0 0 0];
    [0 0 0 -1 1 1 0];
    [0 0 0 0 -1 1 0];
    [0 0 0 1 0 -1 -1];
    [0 0 1 0 0 0 0]
  ]
  # parameter vector
  parms = [k1, k2, k3, k4, k5, k6, k7, k8]

  return x0, nu, parms
end

x0, nu, parms = gillespie_yeast()
F = F_yeast

ext: jl


text:
using Phylodendron2
using Test

@testset "Linking and unlinking" begin
	a = Node("A")
	b = Node("B")
	c = Node("C")
	d = Node("D")
	e = Node("E")

	br_ab = Branch()

	link!(a, b, br_ab)
	link!(a, c)
	link!(c, d)
	link!(c, e)

	@test length(a.links) == 2
	@test a.links[1].to == b
	@test a.links[2].to == c
	@test length(b.links) == 1
	@test b.links[1].to == a
	@test length(c.links) == 3
	@test c.links[1].to == a
	@test c.links[2].to == d
	@test c.links[3].to == e
	@test length(d.links) == 1
	@test d.links[1].to == c
	@test length(e.links) == 1
	@test e.links[1].to == c

	@test unlink!(a, b) == br_ab
	@test isempty(b.links)
	@test length(a.links) == 1
	@test a.links[1].to == c
end # testset "Linking and unlinking"


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.7.0-beta3
#     language: julia
#     name: julia-1.7
# ---

# %%
using Statistics
using Plots

function mcpi1(n)
    c = 0
    for _ in 1:n
        c += rand()^2 + rand()^2 ≤ 1
    end
    4c/n
end

function mcpi2(n)
    c = 0
    for _ in 1:n
        x, y = rand(), rand()
        c += (x^2 + y^2 ≤ 1 && (1 - x)^2 + (1 - y)^2 ≤ 1)
    end
    2c/n + 2
end

function mcpi3(n)
    c = 0
    for _ in 1:n
        x, y = rand(), (1 - √3/2)*rand() + √3/2
        c += (x^2 + y^2 > 1 && (1 - x)^2 + y^2 > 1)
    end
    6(1 - √3/4 - (1 - √3/2)*c/n)
end

function sim(mcpi, n, L)
    X = Vector{Float64}(undef, L)
    Threads.@threads for i in 1:L
        X[i] = mcpi(n)
    end
    X
end

@show Threads.nthreads()
@time pi1 = sim(mcpi1, 10^6, 10^4)
@time pi2 = sim(mcpi2, 10^6, 10^4)
@time pi3 = sim(mcpi3, 10^6, 10^4)

@show mean(pi1), std(pi1)
@show mean(pi2), std(pi2)
@show mean(pi3), std(pi3)

plot()
stephist!(pi1; norm=true, label="mcpi1", lw=1.5, ls=:solid)
stephist!(pi2; norm=true, label="mcpi2", lw=1.5, ls=:dash)
stephist!(pi3; norm=true, label="mcpi3", lw=1.5, ls=:dashdot)

# %%
using Base64
showimg(mime, fn; tag="img") = open(fn) do f
    base64 = base64encode(f)
    display("text/html", """<$tag src="data:$mime;base64,$base64" />""")
end
showimg("image/jpeg", "1BD08057-D77D-44FF-9D3A-7620DA4CF740.jpeg"; tag="img width=80%")

# %%

max_stars_repo_path: 0013/Monte Carlo calculation of pi.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0013/Monte Carlo calculation of pi.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0013/Monte Carlo calculation of pi.jl

max_forks_repo_name: genkuroki/public


text:
module ScatteringTransform

using ChainRules
using CUDA
using Flux: mse, update!
using FFTW
using Wavelets, ContinuousWavelets
using Zygote, Flux, Shearlab, LinearAlgebra, AbstractFFTs
using Flux
using FourierFilterFlux
using Adapt
using RecipesBase
using Base: tail
using ChainRulesCore
using Plots # who cares about weight really?
using Statistics, LineSearches, Optim
using Dates

import Adapt: adapt
import ChainRules: rrule
import Zygote: has_chain_rrule, rrule
import Wavelets: eltypes

include("shared.jl")
export scatteringTransform, stFlux, stParallel, eltypes, depth
include("pathLocs.jl")
export pathLocs
include("scattered.jl")
export Scattered, ScatteredFull, ScatteredOut, pathLocs, nonZeroPaths
include("parallel/basicTypes.jl")
include("Flux/flPool.jl")
export RationPool, nPoolDims, outdims, poolSize
include("Flux/flTransform.jl")
export normalize
export cu

function scatteringTransform(inputSize, m, backend::UnionAll; kwargs...)
    backend(inputSize, m; kwargs...)
end

scatteringTransform(inputSize, m; kwargs...) =
    scatteringTransform(inputSize, m, stFlux; kwargs...)

include("Flux/flUtilities.jl")
export getWavelets, flatten, roll, importantCoords, batchOff, getParameters, getMeanFreq
include("parallel/parallelCore.jl") # there's enough weird stuff going on in
# here that I'm isolating it in a module
using .parallel
function (stPara::stParallel)(x; nonlinearity = abs, thin = false,
    outputSubsample = (-1, -1), subsam = true,
    totalScales = [-1 for i = 1:depth(stPara)+1],
    percentage = 0.9, fftPlans = -1)
    st(x, stPara, nonlinearity; thin = thin, outputSubsample = outputSubsample,
        subsam = subsam, totalScales = totalScales, percentage = percentage,
        fftPlans = fftPlans)
end

function ScatteredFull(layers::scatteringTransform{S,1}, X::Array{T,N};
    totalScales = [-1 for i = 1:depth(layers)+1],
    outputSubsample = (-1, -1)) where {T<:Real,N,S}
    if N == 1
        X = reshape(X, (size(X)..., 1))
    end

    n, q, dataSizes, outputSizes, resultingSize =
        parallel.calculateSizes(layers, outputSubsample, size(X),
            totalScales = totalScales)
    numInLayer = getQ(layers, n, totalScales; product = false)
    addedLayers = [numInLayer[min(i, 2):i] for i = 1:depth(layers)+1]
    if 1 == N
        zerr = [zeros(T, n[i], addedLayers[i]...) for i = 1:depth(layers)+1]
        output = [zeros(T, resultingSize[i], addedLayers[i]...) for i = 1:depth(layers)+1]
    else
        zerr = [zeros(T, n[i], q[i], size(X)[2:end]...) for i = 1:depth(layers)+1]
        output = [zeros(T, resultingSize[i], addedLayers[i]..., size(X)[2:end]...)
                  for i = 1:depth(layers)+1]
        @info "" [size(x) for x in output]
    end
    zerr[1][:, 1, axes(X)[2:end]...] = copy(X)
    return ScatteredFull{T,N + 1}(depth(layers), 1, zerr, output)
end

function ScatteredOut(layers::ST, X::Array{T,N};
    totalScales = [-1 for i = 1:depth(layers)+1],
    outputSubsample = (-1, -1)) where {ST<:scatteringTransform,T<:Real,N,S}
    if N == 1
        X = reshape(X, (size(X)..., 1))
    end

    n, q, dataSizes, outputSizes, resultingSize =
        parallel.calculateSizes(layers, outputSubsample, size(X), totalScales = totalScales)
    addedLayers = parallel.getListOfScaleDims(layers, n, totalScales)
    @info addedLayers
    if 1 == N
        output = [zeros(T, resultingSize[i], addedLayers[i]...) for i = 1:depth(layers)+1]
    else
        output = [zeros(T, resultingSize[i], addedLayers[i]..., size(X)[2:end]...)
                  for i = 1:depth(layers)+1]
    end
    @info [size(x) for x in output]
    return ScatteredOut{T,N + 1}(depth(layers), 1, output)
end

roll(toRoll, stP::stParallel, originalInput; varargs...) = parallel.roll(stP, toRoll, originalInput; varargs...)
export roll, wrap, sizes, calculateThinStSizes, createFFTPlans,
    createRemoteFFTPlan, computeAllWavelets, plotAllWavelets
export spInverse, aTanh, Tanh, ReLU, piecewiseLinear, plInverse
export st
export pathToThinIndex, MatrixAggregator, plotCoordinate, reshapeFlattened,
    numberSkipped, maxPooling, numScales, incrementKeeper, numInLayer
export loadSyntheticMatFile, transformFolder, flatten
include("Flux/adjoints.jl")
include("Flux/interpretationTools.jl")
export ∇st, plotFirstLayer1D, gifFirstLayer, plotSecondLayer, addNextPath, jointPlot
end # end Module


text:
"""
    add_taper!(structure, length)

Adds a tapered segment of `length`. The variables for the parameters
(dispersion, nonlinearity, etc.) in the tapered segment are linearly
interpolated between the parameters for the modes immediately preceding
and succeeding the tapered region. The medium must be the same both
before and after the taper

Not implemented
"""
function add_taper!(structure, length)
    error("add_taper!: Not implemented")
end


text:
using Plots, LaTeXStrings, ProgressBars
cd(dirname(@__FILE__))
function LogisticMap(x₀::Vector{Float64} , r::Float64, n::Real)
    x = x₀
    for i in 1:n
    x = 4 .* r .* x .* (1 .- x)
    end
    return x
end

####Part A####
r_s=collect(range(0.0, 1.0, step=10^(-4)))
x₀_s = rand(100)
x_s=zeros(length(r_s), length(x₀_s))

for (i,r) in enumerate(r_s)
    x_s[i,:]=LogisticMap(x₀_s, r, 10^4 )
end

scatter(r_s, x_s, legend=false, ms=0.5, color=:black, dpi=400, xlabel=L"r", ylabel=L"x", title=L"Logistic\ Map\ :\ x_n_+_1\ =\ 4rx_n(1-x_n)")
savefig("../../computational_physics/PSet9/Figs/Bifurcation.pdf")

####Part B####
####Zoomed Plot####
r_zoom= collect(range(0.70, 1.0, step=10^(-4)))
x₀_zoom = rand(range(0.70, 1.0, length=100), 100)
x_zoom=zeros(length(r_zoom), length(x₀_zoom))
for (i,r) in enumerate(r_zoom)
    x_zoom[i,:]=LogisticMap(x₀_zoom, r, 10^4 )
end
scatter(r_zoom, x_zoom, legend=false, ms=0.5, color=:black, dpi=400, xlabel=L"r", ylabel=L"x",
    title=L"Logistic\ Map\ :\ x_n_+_1\ =\ 4rx_n(1-x_n),\ 0.7<r<1.0")
savefig("../../computational_physics/PSet9/Figs/Zoomed_Bifurcation.pdf")
####Part C####
####Feigenbaum constants#####

#### σ ####
rb= collect(range(0.85, 0.99, step=10^(-4)))
x₀=rand(range(0.88, 0.99, length=256), 256)
x=zeros(length(rb), length(x₀))
for (i,r) in ProgressBar(enumerate(rb))
    x[i,:]=LogisticMap(x₀, r, 10^4 )
end
x
dₙ₋₂_index = findfirst(i->length(unique(round.(x[i, :], digits=9))) == 4, 1:1401)
dₙ₋₁_index = findfirst(i->length(unique(round.(x[i, :], digits=9))) == 8, 1:1401)
dₙ_index = findfirst(i->length(unique(round.(x[i, :], digits=9))) == 16, 1:1401)
δ=(rb[dₙ₋₁_index]-rb[dₙ₋₂_index])/(rb[dₙ_index]-rb[dₙ₋₁_index])
# σ ∼ 5.0999
#### α ####
dₙ₊₁_index = findfirst(i->length(unique(round.(x[i, :], digits=9))) == 32, 1:1401)
rₙ₋₂_index = findfirst(i-> 1/2 ∈ unique(ceil.(x[i, :], digits=3)), dₙ₋₂_index:dₙ₋₁_index) + dₙ₋₂_index
rₙ₋₁_index = findfirst(i-> 1/2 ∈ unique(ceil.(x[i, :], digits=3)),dₙ₋₁_index:dₙ_index) + dₙ₋₁_index
rₙ_index = findfirst(i-> 1/2 ∈ unique(ceil.(x[i, :], digits=3)), dₙ_index:dₙ₊₁_index) + dₙ_index
xsₙ₋₂ = sort(unique(round.(x[rₙ₋₂_index, :], digits=3)))
xsₙ₋₁ = sort(unique(round.(x[rₙ₋₁_index, :], digits=3)))
xsₙ = sort(unique(round.(x[rₙ_index, :], digits=3)))
α = (xsₙ₋₁[2] - xsₙ₋₁[1] - (xsₙ₋₂[4]- xsₙ₋₂[3] ))/( xsₙ[2] - xsₙ[1] - (xsₙ₋₁[2] - xsₙ₋₁[1]))


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sKubeAggregatorPkgApisApiregistrationV1APIServiceSpec <: SwaggerModel
    caBundle::Any # spec type: Union{ Nothing, Vector{UInt8} } # spec name: caBundle
    group::Any # spec type: Union{ Nothing, String } # spec name: group
    groupPriorityMinimum::Any # spec type: Union{ Nothing, Int32 } # spec name: groupPriorityMinimum
    insecureSkipTLSVerify::Any # spec type: Union{ Nothing, Bool } # spec name: insecureSkipTLSVerify
    service::Any # spec type: Union{ Nothing, IoK8sKubeAggregatorPkgApisApiregistrationV1ServiceReference } # spec name: service
    version::Any # spec type: Union{ Nothing, String } # spec name: version
    versionPriority::Any # spec type: Union{ Nothing, Int32 } # spec name: versionPriority

    function IoK8sKubeAggregatorPkgApisApiregistrationV1APIServiceSpec(;caBundle=nothing, group=nothing, groupPriorityMinimum=nothing, insecureSkipTLSVerify=nothing, service=nothing, version=nothing, versionPriority=nothing)
        o = new()
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1APIServiceSpec, Symbol("caBundle"), caBundle)
        setfield!(o, Symbol("caBundle"), caBundle)
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1APIServiceSpec, Symbol("group"), group)
        setfield!(o, Symbol("group"), group)
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1APIServiceSpec, Symbol("groupPriorityMinimum"), groupPriorityMinimum)
        setfield!(o, Symbol("groupPriorityMinimum"), groupPriorityMinimum)
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1APIServiceSpec, Symbol("insecureSkipTLSVerify"), insecureSkipTLSVerify)
        setfield!(o, Symbol("insecureSkipTLSVerify"), insecureSkipTLSVerify)
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1APIServiceSpec, Symbol("service"), service)
        setfield!(o, Symbol("service"), service)
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1APIServiceSpec, Symbol("version"), version)
        setfield!(o, Symbol("version"), version)
        validate_property(IoK8sKubeAggregatorPkgApisApiregistrationV1APIServiceSpec, Symbol("versionPriority"), versionPriority)
        setfield!(o, Symbol("versionPriority"), versionPriority)
        o
    end
end # type IoK8sKubeAggregatorPkgApisApiregistrationV1APIServiceSpec

const _property_map_IoK8sKubeAggregatorPkgApisApiregistrationV1APIServiceSpec = Dict{Symbol,Symbol}(Symbol("caBundle")=>Symbol("caBundle"), Symbol("group")=>Symbol("group"), Symbol("groupPriorityMinimum")=>Symbol("groupPriorityMinimum"), Symbol("insecureSkipTLSVerify")=>Symbol("insecureSkipTLSVerify"), Symbol("service")=>Symbol("service"), Symbol("version")=>Symbol("version"), Symbol("versionPriority")=>Symbol("versionPriority"))
const _property_types_IoK8sKubeAggregatorPkgApisApiregistrationV1APIServiceSpec = Dict{Symbol,String}(Symbol("caBundle")=>"Vector{UInt8}", Symbol("group")=>"String", Symbol("groupPriorityMinimum")=>"Int32", Symbol("insecureSkipTLSVerify")=>"Bool", Symbol("service")=>"IoK8sKubeAggregatorPkgApisApiregistrationV1ServiceReference", Symbol("version")=>"String", Symbol("versionPriority")=>"Int32")
Base.propertynames(::Type{ IoK8sKubeAggregatorPkgApisApiregistrationV1APIServiceSpec }) = collect(keys(_property_map_IoK8sKubeAggregatorPkgApisApiregistrationV1APIServiceSpec))
Swagger.property_type(::Type{ IoK8sKubeAggregatorPkgApisApiregistrationV1APIServiceSpec }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sKubeAggregatorPkgApisApiregistrationV1APIServiceSpec[name]))}
Swagger.field_name(::Type{ IoK8sKubeAggregatorPkgApisApiregistrationV1APIServiceSpec }, property_name::Symbol) =  _property_map_IoK8sKubeAggregatorPkgApisApiregistrationV1APIServiceSpec[property_name]

function check_required(o::IoK8sKubeAggregatorPkgApisApiregistrationV1APIServiceSpec)
    (getproperty(o, Symbol("groupPriorityMinimum")) === nothing) && (return false)
    (getproperty(o, Symbol("service")) === nothing) && (return false)
    (getproperty(o, Symbol("versionPriority")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sKubeAggregatorPkgApisApiregistrationV1APIServiceSpec }, name::Symbol, val)
    if name === Symbol("caBundle")
    end
end


text:
function print_colored(in::String, color1)
    if (VERSION.major == 0 && VERSION.minor <= 6)
        return print_with_color(color1, in)
    else
        return printstyled(in, color = color1)
    end
end

test_list = ["testRealNumbers.jl"; "testComplexNumbers.jl"; "testForwardDiff.jl"; "testDual_.jl"; "testHyperDualNumbers.jl"; "testDates.jl"]

println("Running tests:\n")
for (current_test, i) in zip(test_list, 1:length(test_list))
    println("------------------------------------------------------------")
    println("  * $(current_test) *")
    include(current_test)
    println("------------------------------------------------------------")
    if (i < length(test_list))
        println("")
    end
end


text:
function berechneallgemeineSchnittpunkte(PX1,PX2,theta,ParameterdarstellungX1,ParameterdarstellungX2,a,b,X1,X2,N,JJ,VX1,VX2,HX1,HX2,counter)

pl = plot(leg=false,xticks = -3.5:1:3.5,yticks = -3.5:1:3.5,ratio = :equal);

SPX1 = [];
SPX2 = [];
alleSchnittpunkteX = [];
alleSchnittpunkteY = [];
t = a:0.001:b;

Xrows,Xcols = size(X1);
Yrows,Ycols = size(X2);
Strahl = 1;
x = 7.3916;
Z = 2;
S = 4;
k = 2;
r = 2.741;
eps = 1e-7;

            if almostequal(VX1(k), X1[k])
                if almostequal(VX2(k,Z,S), tand(theta[S])*X1[k]+PX2[Z,S]-tand(theta[S])*PX1[Z,S])
                    if almostequal(HX1(k,Z,S), (X2[k]-PX2[Z,S]+tand(theta[S])*PX1[Z,S])/tand(theta[S]))
                        if almostequal(HX2(k), X2[k])
                            if counter == 1
                                println("Ihr habt die Schnittpunkte mit den Pixelkanten korrekt berechnet!");
                            end
                            for S = 1:size(PX1,2)
                                for Z = 1:size(PX1,1)
                                    p = zeros(size(t,2),1);
                                    y = zeros(size(t,2),1);
                                    HorzSchnitt = zeros(size(X1,2),2);
                                    VertSchnitt = zeros(size(X1,2),2);
                                    Schnittpunktesortiert = [];
                                    counter1 = 1;
                                    if theta[S] == 90
                                    plot!([PX1[Z,S]; PX1[Z,S]],[a; b],linewidth=1.5,color="black",reuse=true);
                                    else
                                        for k = 1:size(t,2)
                                            if tand(theta[S])*t[k]+PX2[Z,S]-tand(theta[S])*PX1[Z,S] >= a && tand(theta[S])*t[k]+PX2[Z,S]-tand(theta[S])*PX1[Z,S] <= b
                                                y[counter1] = tand(theta[S])*t[k]+PX2[Z,S]-tand(theta[S])*PX1[Z,S];
                                                p[counter1] = t[k];
                                                counter1 = counter1+1;
                                            end
                                        end
                                    end
                                    plot!(p,y,linewidth=1.5,color="black",reuse=true);
                                    for k = 1:size(X1,2)
                                        if theta[S] == 90
                                            HorzSchnitt[k,1] = PX1[Z,S];
                                            HorzSchnitt[k,2] = X2[k];
                                        else
                                            if tand(theta[S])*X1[k]+PX2[Z,S]-tand(theta[S])*PX1[Z,S] >= a && tand(theta[S])*X1[k]+PX2[Z,S]-tand(theta[S])*PX1[Z,S] <= b
                                                VertSchnitt[k,1] = X1[k];
                                                VertSchnitt[k,2] = tand(theta[S])*X1[k]+PX2[Z,S]-tand(theta[S])*PX1[Z,S];
                                            end
                                            if (X2[k]-PX2[Z,S]+tand(theta[S])*PX1[Z,S])/tand(theta[S]) >= a && (X2[k]-PX2[Z,S]+tand(theta[S])*PX1[Z,S])/tand(theta[S]) <= b
                                                HorzSchnitt[k,1] = (X2[k]-PX2[Z,S]+tand(theta[S])*PX1[Z,S])/tand(theta[S]);
                                                HorzSchnitt[k,2] = X2[k];
                                            end
                                        end
                                    end

                                    Schnittpunkte = vcat(HorzSchnitt, VertSchnitt);
                                    Schnittpunkteneu1 = [];
                                    Schnittpunkteneu2 = [];
                                    Zeile = 1
                                    for q = 1:size(Schnittpunkte,1)
                                        if Schnittpunkte[q,1] != 0 && Schnittpunkte[q,2] != 0
                                            Schnittpunkteneu1 = resize(Schnittpunkteneu1,Zeile);
                                            Schnittpunkteneu1[Zeile] = Schnittpunkte[q,1];
                                            Schnittpunkteneu2 = resize(Schnittpunkteneu2,Zeile);
                                            Schnittpunkteneu2[Zeile] = Schnittpunkte[q,2];
                                            Zeile = Zeile+1;
                                        end
                                    end
                                    Schnittpunkteecht = hcat(Schnittpunkteneu1, Schnittpunkteneu2);
                                    if theta[S] < 90
                                        Schnittpunktesortiert = sortieren(Schnittpunkteecht,1,0);
                                    elseif theta[S] == 90
                                        Schnittpunktesortiert = sortieren(Schnittpunkteecht,2,0);
                                    else
                                        Schnittpunktesortiert = sortieren(Schnittpunkteecht,1,1);
                                    end
                                    X1Teil = Schnittpunktesortiert[:,1];
                                    X2Teil = Schnittpunktesortiert[:,2];
                                    plot!(X1Teil,X2Teil,linewidth=1.5,color="black",markersize=18,reuse=true);
                                    scatter!(X1Teil,X2Teil,markersize=3,markercolor=:orange)
                                    X1Teil = resize(X1Teil,size(Schnittpunkte,1));
                                    X2Teil = resize(X2Teil,size(Schnittpunkte,1));
                                    Schnittvektorsortiert = hcat(X1Teil, X2Teil);
                                        #display(Schnittvektorsortiert)
                                    Schnittrows,Schnittcols = size(Schnittpunkte);
                                    diff = Ycols+Xcols-Schnittrows; #damit alle Anzahln von Schnittpunkten von der Länge der Zeilen übereinstimmen
                                    if diff > 0
                                        Nullen = zeros(diff,2);
                                        Schnittnull = vcat(Schnittpunktesortiert, Nullen);
                                        if Strahl == 1
                                            alleSchnittpunkteX = round.(Schnittnull[:1];digits = 4);
                                            alleSchnittpunkteY = round.(Schnittnull[:,2];digits = 4);
                                            Strahl = Strahl+1;
                                        else
                                            alleSchnittpunkteX = round.(hcat(alleSchnittpunkteX, Schnittnull[:,1]);digits = 4);
                                            alleSchnittpunkteY = round.(hcat(alleSchnittpunkteY, Schnittnull[:,2]);digits = 4);
                                        end
                                    else
                                        if Strahl == 1
                                            Schnittpunktesortiert1 = resize(Schnittpunktesortiert[:,1],size(Schnittpunkte,1));
                                            Schnittpunktesortiert2 = resize(Schnittpunktesortiert[:,2],size(Schnittpunkte,1));
                                            alleSchnittpunkteX = round.(Schnittpunktesortiert1;digits = 4);
                                            alleSchnittpunkteY = round.(Schnittpunktesortiert2;digits = 4);
                                            Strahl = Strahl+1;
                                        else
                                            SpalteX1 = Schnittpunktesortiert[:,1];
                                            SpalteX1 = resize(SpalteX1,size(Schnittpunkte,1));
                                            SpalteY2 = Schnittpunktesortiert[:,2];
                                            SpalteY2 = resize(SpalteY2,size(Schnittpunkte,1));
                                            alleSchnittpunkteX = round.(hcat(alleSchnittpunkteX, SpalteX1);digits=4);
                                            alleSchnittpunkteY = round.(hcat(alleSchnittpunkteY, SpalteY2);digits=4);
                                        end
                                        #display(alleSchnittpunkteX)
                                    end
                                end
                            end

                            plot!([-1.5; 1.5],[0.5; 0.5],linewidth=1.5,color="blue",reuse=true);
                            plot!([-1.5; 1.5],[-0.5; -0.5],linewidth=1.5,color="blue",reuse=true);
                            plot!([-1.5; 1.5],[1.5; 1.5],linewidth=1.5,color="blue",reuse=true);
                            plot!([-1.5; 1.5],[-1.5; -1.5],linewidth=1.5,color="blue",reuse=true);
                            plot!([-1.5; -1.5],[-1.5; 1.5],linewidth=1.5,color="blue",reuse=true);
                            plot!([1.5; 1.5],[-1.5; 1.5],linewidth=1.5,color="blue",reuse=true);
                            plot!([0.5; 0.5],[-1.5; 1.5],linewidth=1.5,color="blue",reuse=true);
                            plot!([-0.5; -0.5],[-1.5; 1.5],linewidth=1.5,color="blue",reuse=true);
                            if counter == 1
                                plot!([-3.5; 3.5],[0.25; 0.25],linewidth=1.5,color="black",reuse=true);
                                plot!([-3.5; 3.5],[-0.75; -0.75],linewidth=1.5,color="black",reuse=true);
                                plot!([-3.5; 3.5],[1.25; 1.25],linewidth=1.5,color="black",reuse=true);
                            end
                            #axis equal
                            #grid on
                            #
                            #ax.XGrid = 'on';
                            #
                            #ay.XGrid = 'on';
                            rowsx = size(alleSchnittpunkteX,1);
                            for Spalte1 = 1:size(alleSchnittpunkteX,2)
                                doppelt = 1;
                                for Zeile1 = 1:rowsx-1
                                    if almostequal(alleSchnittpunkteX[Zeile1,Spalte1], alleSchnittpunkteX[Zeile1+1,Spalte1]) && almostequal(alleSchnittpunkteY[Zeile1,Spalte1], alleSchnittpunkteY[Zeile1+1,Spalte1]) && alleSchnittpunkteX[Zeile1+1,Spalte1] != 0 && alleSchnittpunkteY[Zeile1+1,Spalte1] != 0
                                        for hoch = 1:(size(alleSchnittpunkteX,1)-1-Zeile1)
                                            #if alleSchnittpunkteX[Zeile1+hoch,Spalte1] != 0 || alleSchnittpunkteY[Zeile1+hoch,Spalte1] != 0
                                                alleSchnittpunkteX[Zeile1+hoch,Spalte1] = alleSchnittpunkteX[Zeile1+hoch+1,Spalte1];
                                                alleSchnittpunkteY[Zeile1+hoch,Spalte1] = alleSchnittpunkteY[Zeile1+hoch+1,Spalte1];
                                                                        #
                                            #end
                                        end
                                        alleSchnittpunkteX[rowsx-doppelt+1,Spalte1] = 0;
                                        alleSchnittpunkteY[rowsx-doppelt+1,Spalte1] = 0;
                                        doppelt = doppelt+1;
                                    end
                                end
                                #if doppelt != 1
                                #    for l = rowsx-doppelt-1:rowsx
                                #        alleSchnittpunkteX[l,Spalte1] = 0;
                                #        alleSchnittpunkteY[l,Spalte1] = 0;
                                #    end
                                #end
                            end

                            SPX1 = alleSchnittpunkteX[1:rowsx,:];
                            SPX2 = alleSchnittpunkteY[1:rowsx,:];

                        else
                            println("Überprüft erneut eure Eingabe für die X2 Koordinate der Schnittpunkte mit den horizontalen Pixelkanten.");
                            SPX1 = [];
                            SPX2 = [];
                        end
                    else
                        if almostequal(HX2(k), X2[k])
                            println("Überprüft erneut eure Eingabe für die X1 Koordinate der Schnittpunkte mit den horizontalen Pixelkanten.");
                            SPX1 = [];
                            SPX2 = [];
                        else
                            println("Überprüft erneut eure Eingaben für die X1 und X2 Koordinaten der Schnittpunkte mit den horizontalen Pixelkanten.");
                            SPX1 = [];
                            SPX2 = [];
                        end
                    end
                else
                    if almostequal(HX1(k,Z,S), (X2[k]-PX2[Z,S]+tand(theta[S])*PX1[Z,S])/tand(theta[S]))
                        if almostequal(HX2(k), X2[k])
                            println("Überprüft erneut eure Eingabe für die X2 Koordinate der Schnittpunkte mit den vertikalem Pixelkanten.");
                            SPX1 = [];
                            SPX2 = [];
                        else
                            println("Überprüft erneut eure Eingabe für die X2 Koordinaten der Schnittpunkte mit den horizontalen und vertikalen Pixelkanten.");
                            SPX1 = [];
                            SPX2 = [];
                        end
                    else
                        if almostequal(HX2(k), X2[k])
                            println("Überprüft erneut eure Eingabe für die X1 Koordinate der Schnittpunkte mit den horizontalen und die X2 Koordinate der vertikalen Pixelkanten.");
                            SPX1 = [];
                            SPX2 = [];
                        else
                            println("Überprüft erneut eure Eingaben für die X1 und X2 Koordinaten der Schnittpunkte mit den horizontalen und die X2 Koordinate der vertikalen Pixelkanten.");
                            SPX1 = [];
                            SPX2 = [];
                        end
                    end
                end
            else
                if almostequal(VX2(k,Z,S), tand(theta[S])*X1[k]+PX2[Z,S]-tand(theta[S])*PX1[Z,S])
                    if almostequal(HX1(k,Z,S), (X2[k]-PX2[Z,S]+tand(theta[S])*PX1[Z,S])/tand(theta[S]))
                        if almostequal(HX2(k), X2[k])
                            println("Überprüft erneut eure Eingabe für die X1 Koordinate der Schnittpunkte mit den vertikalen Pixelkanten.");
                            SPX1 = [];
                            SPX2 = [];
                        else
                            println("Überprüft erneut eure Eingabe für die X2 Koordinate der Schnittpunkte mit den horizontalen und die X1 Koordinate für die vertikalen Pixelkanten.");
                            SPX1 = [];
                            SPX2 = [];
                        end
                    else
                        if almostequal(HX2(k),X2[k])
                            println("Überprüft erneut eure Eingabe für die X1 Koordinate der Schnittpunkte mit den horizontalen und vertikalen Pixelkanten.");
                            SPX1 = [];
                            SPX2 = [];
                        else
                            println("Überprüft erneut eure Eingaben für die X1 und X2 Koordinaten der Schnittpunkte mit den horizontalen und die X1 Koordinate für die vertikalen Pixelkanten.");
                            SPX1 = [];
                            SPX2 = [];
                        end
                    end
                else
                    if almostequal(HX1(k,Z,S),(X2[k]-PX2[Z,S]+tand(theta[S])*PX1[Z,S])/tand(theta[S]))
                        if almostequal(HX2(k), X2[k])
                            println("Überprüft erneut eure Eingabe für die X2 Koordinate der Schnittpunkte mit den vertikalem und die X1 Koordinate für die vertikalen Pixelkanten.");
                            SPX1 = [];
                            SPX2 = [];
                        else
                            println("Überprüft erneut eure Eingabe für die X1 und X2 Koordinaten der Schnittpunkte mit den vertikalen und die X2 Koordinate für die horizontalen Pixelkanten.");
                            SPX1 = [];
                            SPX2 = [];
                        end
                    else
                        if almostequal(HX2(k), X2[k])
                            println("Überprüft erneut eure Eingabe für die X1 Koordinate der Schnittpunkte mit den horizontalen und vertikalen Pixelkanten und der X2 Koordinate der vertikalen Pixelkanten.");
                            SPX1 = [];
                            SPX2 = [];
                        else
                            println("All eure Formeln sind leider unkorrekt.");
                            SPX1 = [];
                            SPX2 = [];
                        end
                    end
                end
            end

return SPX1, SPX2,pl;

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

abstract type AbstractCmd end

# libuv process option flags
const UV_PROCESS_WINDOWS_VERBATIM_ARGUMENTS = UInt8(1 << 2)
const UV_PROCESS_DETACHED = UInt8(1 << 3)
const UV_PROCESS_WINDOWS_HIDE = UInt8(1 << 4)

struct Cmd <: AbstractCmd
    exec::Vector{String}
    ignorestatus::Bool
    flags::UInt32 # libuv process flags
    env::Union{Array{String},Nothing}
    dir::String
    Cmd(exec::Vector{String}) =
        new(exec, false, 0x00, nothing, "")
    Cmd(cmd::Cmd, ignorestatus, flags, env, dir) =
        new(cmd.exec, ignorestatus, flags, env,
            dir === cmd.dir ? dir : cstr(dir))
    function Cmd(cmd::Cmd; ignorestatus::Bool=cmd.ignorestatus, env=cmd.env, dir::AbstractString=cmd.dir,
                 detach::Bool = 0 != cmd.flags & UV_PROCESS_DETACHED,
                 windows_verbatim::Bool = 0 != cmd.flags & UV_PROCESS_WINDOWS_VERBATIM_ARGUMENTS,
                 windows_hide::Bool = 0 != cmd.flags & UV_PROCESS_WINDOWS_HIDE)
        flags = detach * UV_PROCESS_DETACHED |
                windows_verbatim * UV_PROCESS_WINDOWS_VERBATIM_ARGUMENTS |
                windows_hide * UV_PROCESS_WINDOWS_HIDE
        new(cmd.exec, ignorestatus, flags, byteenv(env),
            dir === cmd.dir ? dir : cstr(dir))
    end
end

"""
    Cmd(cmd::Cmd; ignorestatus, detach, windows_verbatim, windows_hide, env, dir)

Construct a new `Cmd` object, representing an external program and arguments, from `cmd`,
while changing the settings of the optional keyword arguments:

* `ignorestatus::Bool`: If `true` (defaults to `false`), then the `Cmd` will not throw an
  error if the return code is nonzero.
* `detach::Bool`: If `true` (defaults to `false`), then the `Cmd` will be run in a new
  process group, allowing it to outlive the `julia` process and not have Ctrl-C passed to
  it.
* `windows_verbatim::Bool`: If `true` (defaults to `false`), then on Windows the `Cmd` will
  send a command-line string to the process with no quoting or escaping of arguments, even
  arguments containing spaces. (On Windows, arguments are sent to a program as a single
  "command-line" string, and programs are responsible for parsing it into arguments. By
  default, empty arguments and arguments with spaces or tabs are quoted with double quotes
  `"` in the command line, and `\\` or `"` are preceded by backslashes.
  `windows_verbatim=true` is useful for launching programs that parse their command line in
  nonstandard ways.) Has no effect on non-Windows systems.
* `windows_hide::Bool`: If `true` (defaults to `false`), then on Windows no new console
  window is displayed when the `Cmd` is executed. This has no effect if a console is
  already open or on non-Windows systems.
* `env`: Set environment variables to use when running the `Cmd`. `env` is either a
  dictionary mapping strings to strings, an array of strings of the form `"var=val"`, an
  array or tuple of `"var"=>val` pairs, or `nothing`. In order to modify (rather than
  replace) the existing environment, create `env` by `copy(ENV)` and then set
  `env["var"]=val` as desired.
* `dir::AbstractString`: Specify a working directory for the command (instead
  of the current directory).

For any keywords that are not specified, the current settings from `cmd` are used. Normally,
to create a `Cmd` object in the first place, one uses backticks, e.g.

    Cmd(`echo "Hello world"`, ignorestatus=true, detach=false)
"""
Cmd

hash(x::Cmd, h::UInt) = hash(x.exec, hash(x.env, hash(x.ignorestatus, hash(x.dir, hash(x.flags, h)))))
==(x::Cmd, y::Cmd) = x.exec == y.exec && x.env == y.env && x.ignorestatus == y.ignorestatus &&
                     x.dir == y.dir && isequal(x.flags, y.flags)

struct OrCmds <: AbstractCmd
    a::AbstractCmd
    b::AbstractCmd
    OrCmds(a::AbstractCmd, b::AbstractCmd) = new(a, b)
end

struct ErrOrCmds <: AbstractCmd
    a::AbstractCmd
    b::AbstractCmd
    ErrOrCmds(a::AbstractCmd, b::AbstractCmd) = new(a, b)
end

struct AndCmds <: AbstractCmd
    a::AbstractCmd
    b::AbstractCmd
    AndCmds(a::AbstractCmd, b::AbstractCmd) = new(a, b)
end

hash(x::AndCmds, h::UInt) = hash(x.a, hash(x.b, h))
==(x::AndCmds, y::AndCmds) = x.a == y.a && x.b == y.b

shell_escape(cmd::Cmd; special::AbstractString="") =
    shell_escape(cmd.exec..., special=special)
shell_escape_posixly(cmd::Cmd) =
    shell_escape_posixly(cmd.exec...)

function show(io::IO, cmd::Cmd)
    print_env = cmd.env !== nothing
    print_dir = !isempty(cmd.dir)
    (print_env || print_dir) && print(io, "setenv(")
    print(io, '`')
    join(io, map(cmd.exec) do arg
        replace(sprint(context=io) do io
            with_output_color(:underline, io) do io
                print_shell_word(io, arg, shell_special)
            end
        end, '`' => "\\`")
    end, ' ')
    print(io, '`')
    print_env && (print(io, ","); show(io, cmd.env))
    print_dir && (print(io, "; dir="); show(io, cmd.dir))
    (print_dir || print_env) && print(io, ")")
end

function show(io::IO, cmds::Union{OrCmds,ErrOrCmds})
    print(io, "pipeline(")
    show(io, cmds.a)
    print(io, ", ")
    print(io, isa(cmds, ErrOrCmds) ? "stderr=" : "stdout=")
    show(io, cmds.b)
    print(io, ")")
end

function show(io::IO, cmds::AndCmds)
    show(io, cmds.a)
    print(io, " & ")
    show(io, cmds.b)
end

const STDIN_NO  = 0
const STDOUT_NO = 1
const STDERR_NO = 2

struct FileRedirect
    filename::AbstractString
    append::Bool
    function FileRedirect(filename, append)
        if lowercase(filename) == (@static Sys.iswindows() ? "nul" : "/dev/null")
            @warn "For portability use devnull instead of a file redirect" maxlog=1
        end
        return new(filename, append)
    end
end

rawhandle(::DevNull) = C_NULL
rawhandle(x::OS_HANDLE) = x
if OS_HANDLE !== RawFD
    rawhandle(x::RawFD) = Libc._get_osfhandle(x)
end

const Redirectable = Union{IO, FileRedirect, RawFD, OS_HANDLE}
const StdIOSet = NTuple{3, Redirectable}

struct CmdRedirect <: AbstractCmd
    cmd::AbstractCmd
    handle::Redirectable
    stream_no::Int
end

function show(io::IO, cr::CmdRedirect)
    print(io, "pipeline(")
    show(io, cr.cmd)
    print(io, ", ")
    if cr.stream_no == STDOUT_NO
        print(io, "stdout=")
    elseif cr.stream_no == STDERR_NO
        print(io, "stderr=")
    elseif cr.stream_no == STDIN_NO
        print(io, "stdin=")
    end
    show(io, cr.handle)
    print(io, ")")
end

"""
    ignorestatus(command)

Mark a command object so that running it will not throw an error if the result code is non-zero.
"""
ignorestatus(cmd::Cmd) = Cmd(cmd, ignorestatus=true)
ignorestatus(cmd::Union{OrCmds,AndCmds}) =
    typeof(cmd)(ignorestatus(cmd.a), ignorestatus(cmd.b))

"""
    detach(command)

Mark a command object so that it will be run in a new process group, allowing it to outlive the julia process, and not have Ctrl-C interrupts passed to it.
"""
detach(cmd::Cmd) = Cmd(cmd; detach=true)

# like String(s), but throw an error if s contains NUL, since
# libuv requires NUL-terminated strings
function cstr(s)
    if Base.containsnul(s)
        throw(ArgumentError("strings containing NUL cannot be passed to spawned processes"))
    end
    return String(s)
end

# convert various env representations into an array of "key=val" strings
byteenv(env::AbstractArray{<:AbstractString}) =
    String[cstr(x) for x in env]
byteenv(env::AbstractDict) =
    String[cstr(string(k)*"="*string(v)) for (k,v) in env]
byteenv(env::Nothing) = nothing
byteenv(env::Union{AbstractVector{Pair{T}}, Tuple{Vararg{Pair{T}}}}) where {T<:AbstractString} =
    String[cstr(k*"="*string(v)) for (k,v) in env]

"""
    setenv(command::Cmd, env; dir="")

Set environment variables to use when running the given `command`. `env` is either a
dictionary mapping strings to strings, an array of strings of the form `"var=val"`, or zero
or more `"var"=>val` pair arguments. In order to modify (rather than replace) the existing
environment, create `env` by `copy(ENV)` and then setting `env["var"]=val` as desired, or
use `withenv`.

The `dir` keyword argument can be used to specify a working directory for the command.
"""
setenv(cmd::Cmd, env; dir="") = Cmd(cmd; env=byteenv(env), dir=dir)
setenv(cmd::Cmd, env::Pair{<:AbstractString}...; dir="") =
    setenv(cmd, env; dir=dir)
setenv(cmd::Cmd; dir="") = Cmd(cmd; dir=dir)

(&)(left::AbstractCmd, right::AbstractCmd) = AndCmds(left, right)
redir_out(src::AbstractCmd, dest::AbstractCmd) = OrCmds(src, dest)
redir_err(src::AbstractCmd, dest::AbstractCmd) = ErrOrCmds(src, dest)

# Stream Redirects
redir_out(dest::Redirectable, src::AbstractCmd) = CmdRedirect(src, dest, STDIN_NO)
redir_out(src::AbstractCmd, dest::Redirectable) = CmdRedirect(src, dest, STDOUT_NO)
redir_err(src::AbstractCmd, dest::Redirectable) = CmdRedirect(src, dest, STDERR_NO)

# File redirects
redir_out(src::AbstractCmd, dest::AbstractString) = CmdRedirect(src, FileRedirect(dest, false), STDOUT_NO)
redir_out(src::AbstractString, dest::AbstractCmd) = CmdRedirect(dest, FileRedirect(src, false), STDIN_NO)
redir_err(src::AbstractCmd, dest::AbstractString) = CmdRedirect(src, FileRedirect(dest, false), STDERR_NO)
redir_out_append(src::AbstractCmd, dest::AbstractString) = CmdRedirect(src, FileRedirect(dest, true), STDOUT_NO)
redir_err_append(src::AbstractCmd, dest::AbstractString) = CmdRedirect(src, FileRedirect(dest, true), STDERR_NO)

"""
    pipeline(command; stdin, stdout, stderr, append=false)

Redirect I/O to or from the given `command`. Keyword arguments specify which of the
command's streams should be redirected. `append` controls whether file output appends to the
file. This is a more general version of the 2-argument `pipeline` function.
`pipeline(from, to)` is equivalent to `pipeline(from, stdout=to)` when `from` is a command,
and to `pipeline(to, stdin=from)` when `from` is another kind of data source.

**Examples**:

```julia
run(pipeline(`dothings`, stdout="out.txt", stderr="errs.txt"))
run(pipeline(`update`, stdout="log.txt", append=true))
```
"""
function pipeline(cmd::AbstractCmd; stdin=nothing, stdout=nothing, stderr=nothing, append::Bool=false)
    if append && stdout === nothing && stderr === nothing
        error("append set to true, but no output redirections specified")
    end
    if stdin !== nothing
        cmd = redir_out(stdin, cmd)
    end
    if stdout !== nothing
        cmd = append ? redir_out_append(cmd, stdout) : redir_out(cmd, stdout)
    end
    if stderr !== nothing
        cmd = append ? redir_err_append(cmd, stderr) : redir_err(cmd, stderr)
    end
    return cmd
end

pipeline(cmd::AbstractCmd, dest) = pipeline(cmd, stdout=dest)
pipeline(src::Union{Redirectable,AbstractString}, cmd::AbstractCmd) = pipeline(cmd, stdin=src)

"""
    pipeline(from, to, ...)

Create a pipeline from a data source to a destination. The source and destination can be
commands, I/O streams, strings, or results of other `pipeline` calls. At least one argument
must be a command. Strings refer to filenames. When called with more than two arguments,
they are chained together from left to right. For example, `pipeline(a,b,c)` is equivalent to
`pipeline(pipeline(a,b),c)`. This provides a more concise way to specify multi-stage
pipelines.

**Examples**:

```julia
run(pipeline(`ls`, `grep xyz`))
run(pipeline(`ls`, "out.txt"))
run(pipeline("out.txt", `grep xyz`))
```
"""
pipeline(a, b, c, d...) = pipeline(pipeline(a,b), c, d...)

mutable struct Process <: AbstractPipe
    cmd::Cmd
    handle::Ptr{Cvoid}
    in::IO
    out::IO
    err::IO
    exitcode::Int64
    termsignal::Int32
    exitnotify::Condition
    closenotify::Condition
    function Process(cmd::Cmd, handle::Ptr{Cvoid},
                     in::Union{Redirectable, Ptr{Cvoid}},
                     out::Union{Redirectable, Ptr{Cvoid}},
                     err::Union{Redirectable, Ptr{Cvoid}})
        if !isa(in, IO)
            in = devnull
        end
        if !isa(out, IO)
            out = devnull
        end
        if !isa(err, IO)
            err = devnull
        end
        this = new(cmd, handle, in, out, err,
                   typemin(fieldtype(Process, :exitcode)),
                   typemin(fieldtype(Process, :termsignal)),
                   Condition(), Condition())
        finalizer(uvfinalize, this)
        return this
    end
end
pipe_reader(p::Process) = p.out
pipe_writer(p::Process) = p.in

struct ProcessChain <: AbstractPipe
    processes::Vector{Process}
    in::Redirectable
    out::Redirectable
    err::Redirectable
    ProcessChain(stdios::StdIOSet) = new(Process[], stdios[1], stdios[2], stdios[3])
end
pipe_reader(p::ProcessChain) = p.out
pipe_writer(p::ProcessChain) = p.in

function _jl_spawn(file, argv, cmd::Cmd, stdio)
    loop = eventloop()
    handles = Tuple{Cint, UInt}[ # assuming little-endian layout
        let h = rawhandle(io)
            h === C_NULL     ? (0x00, UInt(0)) :
            h isa OS_HANDLE  ? (0x02, UInt(cconvert(@static(Sys.iswindows() ? Ptr{Cvoid} : Cint), h))) :
            h isa Ptr{Cvoid} ? (0x04, UInt(h)) :
            error("invalid spawn handle $h from $io")
        end
        for io in stdio]
    proc = Libc.malloc(_sizeof_uv_process)
    disassociate_julia_struct(proc)
    error = ccall(:jl_spawn, Int32,
              (Cstring, Ptr{Cstring}, Ptr{Cvoid}, Ptr{Cvoid},
               Ptr{Tuple{Cint, UInt}}, Int,
               UInt32, Ptr{Cstring}, Cstring, Ptr{Cvoid}),
        file, argv, loop, proc,
        handles, length(handles),
        cmd.flags,
        cmd.env === nothing ? C_NULL : cmd.env,
        isempty(cmd.dir) ? C_NULL : cmd.dir,
        uv_jl_return_spawn::Ptr{Cvoid})
    if error != 0
        ccall(:jl_forceclose_uv, Cvoid, (Ptr{Cvoid},), proc)
        throw(_UVError("could not spawn " * string(cmd), error))
    end
    return proc
end

function uvfinalize(proc::Process)
    if proc.handle != C_NULL
        disassociate_julia_struct(proc.handle)
        ccall(:jl_close_uv, Cvoid, (Ptr{Cvoid},), proc.handle)
        proc.handle = C_NULL
    end
    nothing
end

function uv_return_spawn(p::Ptr{Cvoid}, exit_status::Int64, termsignal::Int32)
    data = ccall(:jl_uv_process_data, Ptr{Cvoid}, (Ptr{Cvoid},), p)
    data == C_NULL && return
    proc = unsafe_pointer_to_objref(data)::Process
    proc.exitcode = exit_status
    proc.termsignal = termsignal
    ccall(:jl_close_uv, Cvoid, (Ptr{Cvoid},), proc.handle)
    notify(proc.exitnotify)
    nothing
end

function _uv_hook_close(proc::Process)
    proc.handle = C_NULL
    notify(proc.closenotify)
end

function _spawn(redirect::CmdRedirect, stdios::StdIOSet; chain::Union{ProcessChain, Nothing}=nothing)
    _spawn(redirect.cmd,
           (redirect.stream_no == STDIN_NO  ? redirect.handle : stdios[1],
            redirect.stream_no == STDOUT_NO ? redirect.handle : stdios[2],
            redirect.stream_no == STDERR_NO ? redirect.handle : stdios[3]),
           chain=chain)
end

function _spawn(cmds::OrCmds, stdios::StdIOSet; chain::Union{ProcessChain, Nothing}=nothing)
    if chain === nothing
        chain = ProcessChain(stdios)
    end
    in_pipe, out_pipe = link_pipe(false, false)
    try
        _spawn(cmds.a, (stdios[1], out_pipe, stdios[3]), chain=chain)
        _spawn(cmds.b, (in_pipe, stdios[2], stdios[3]), chain=chain)
    finally
        close_pipe_sync(out_pipe)
        close_pipe_sync(in_pipe)
    end
    return chain
end

function _spawn(cmds::ErrOrCmds, stdios::StdIOSet; chain::Union{ProcessChain, Nothing}=nothing)
    if chain === nothing
        chain = ProcessChain(stdios)
    end
    in_pipe, out_pipe = link_pipe(false, false)
    try
        _spawn(cmds.a, (stdios[1], stdios[2], out_pipe), chain=chain)
        _spawn(cmds.b, (in_pipe, stdios[2], stdios[3]), chain=chain)
    finally
        close_pipe_sync(out_pipe)
        close_pipe_sync(in_pipe)
    end
    return chain
end

function setup_stdio(stdio::PipeEndpoint, child_readable::Bool)
    if stdio.status == StatusInit
        rd, wr = link_pipe(!child_readable, child_readable)
        open_pipe!(stdio, child_readable ? wr : rd, !child_readable, child_readable)
        return (child_readable ? rd : wr, true)
    end
    return (stdio, false)
end

function setup_stdio(stdio::Pipe, child_readable::Bool)
    if stdio.in.status == StatusInit && stdio.out.status == StatusInit
        link_pipe!(stdio)
    end
    io = child_readable ? stdio.out : stdio.in
    return (io, false)
end

setup_stdio(stdio::AbstractPipe, readable::Bool) =
    setup_stdio(readable ? pipe_reader(stdio) : pipe_writer(stdio), readable)

function setup_stdio(stdio::IOStream, child_readable::Bool)
    io = RawFD(fd(stdio))
    return (io, false)
end

function setup_stdio(stdio::FileRedirect, child_readable::Bool)
    if child_readable
        attr = JL_O_RDONLY
        perm = zero(S_IRUSR)
    else
        attr = JL_O_WRONLY | JL_O_CREAT
        attr |= stdio.append ? JL_O_APPEND : JL_O_TRUNC
        perm = S_IRUSR | S_IWUSR | S_IRGRP | S_IROTH
    end
    io = Filesystem.open(stdio.filename, attr, perm)
    return (io, true)
end

function setup_stdio(io, child_readable::Bool)
    # if there is no specialization,
    # assume that rawhandle is defined for it
    return (io, false)
end

close_stdio(stdio::OS_HANDLE) = close_pipe_sync(stdio)
close_stdio(stdio::Nothing) = nothing
close_stdio(stdio) = close(stdio)

function setup_stdio(anon::Function, stdio::StdIOSet)
    in, close_in = setup_stdio(stdio[1], true)
    try
        out, close_out = setup_stdio(stdio[2], false)
        try
            err, close_err = setup_stdio(stdio[3], false)
            try
                anon((in, out, err))
            finally
                close_err && close_stdio(err)
            end
        finally
            close_out && close_stdio(out)
        end
    finally
        close_in && close_stdio(in)
    end
    nothing
end

function _spawn(cmd::Cmd, stdios::StdIOSet; chain::Union{ProcessChain, Nothing}=nothing)
    if isempty(cmd.exec)
        throw(ArgumentError("cannot spawn empty command"))
    end
    pp = Process(cmd, C_NULL, stdios[1], stdios[2], stdios[3])
    setup_stdio(stdios) do stdios
        handle = _jl_spawn(cmd.exec[1], cmd.exec, cmd, stdios)
        associate_julia_struct(handle, pp)
        pp.handle = handle
    end
    if chain !== nothing
        push!(chain.processes, pp)
    end
    return pp
end

function _spawn(cmds::AndCmds, stdios::StdIOSet; chain::Union{ProcessChain, Nothing}=nothing)
    if chain === nothing
        chain = ProcessChain(stdios)
    end
    setup_stdio(stdios) do stdios
        _spawn(cmds.a, stdios, chain=chain)
        _spawn(cmds.b, stdios, chain=chain)
    end
    return chain
end

# INTERNAL
# returns stdios:
# A set of up to 256 stdio instructions, where each entry can be either:
#   | - An IO to be passed to the child
#   | - devnull to pass /dev/null
#   | - An Filesystem.File object to redirect the output to
#   \ - A string specifying a filename to be opened

spawn_opts_swallow(stdios::StdIOSet) = (stdios,)
spawn_opts_swallow(in::Redirectable=devnull, out::Redirectable=devnull, err::Redirectable=devnull, args...) =
    ((in, out, err), args...)
spawn_opts_inherit(stdios::StdIOSet) = (stdios,)
# pass original descriptors to child processes by default, because we might
# have already exhausted and closed the libuv object for our standard streams.
# this caused issue #8529.
spawn_opts_inherit(in::Redirectable=RawFD(0), out::Redirectable=RawFD(1), err::Redirectable=RawFD(2), args...) =
    ((in, out, err), args...)

_spawn(cmds::AbstractCmd, args...; chain::Union{ProcessChain, Nothing}=nothing) =
    _spawn(cmds, spawn_opts_swallow(args...)...; chain=chain)

function eachline(cmd::AbstractCmd; keep::Bool=false)
    _stdout = Pipe()
    processes = _spawn(cmd, (devnull, _stdout, stderr))
    close(_stdout.in)
    out = _stdout.out
    # implicitly close after reading lines, since we opened
    return EachLine(out, keep=keep,
        ondone=()->(close(out); success(processes) || pipeline_error(processes)))::EachLine
end

function open(cmds::AbstractCmd, mode::AbstractString, other::Redirectable=devnull)
    if mode == "r+" || mode == "w+"
        return open(cmds, other, read = true, write = true)
    elseif mode == "r"
        return open(cmds, other)
    elseif mode == "w"
        return open(cmds, other, write = true)
    else
        throw(ArgumentError("mode must be \"r\" or \"w\", not \"$mode\""))
    end
end

# return a Process object to read-to/write-from the pipeline
"""
    open(command, stdio=devnull; write::Bool = false, read::Bool = !write)

Start running `command` asynchronously, and return a tuple `(stream,process)`.  If `read` is
true, then `stream` reads from the process's standard output and `stdio` optionally
specifies the process's standard input stream.  If `write` is true, then `stream` writes to
the process's standard input and `stdio` optionally specifies the process's standard output
stream.
"""
function open(cmds::AbstractCmd, other::Redirectable=devnull; write::Bool = false, read::Bool = !write)
    if read && write
        other === devnull || throw(ArgumentError("no other stream can be specified in read-write mode"))
        in = Pipe()
        out = Pipe()
        processes = _spawn(cmds, (in,out,stderr))
        close(in.out)
        close(out.in)
    elseif read
        in = other
        out = Pipe()
        processes = _spawn(cmds, (in,out,stderr))
        close(out.in)
    elseif write
        in = Pipe()
        out = other
        processes = _spawn(cmds, (in,out,stderr))
        close(in.out)
    else
        processes = _spawn(cmds)
    end
    return processes
end

"""
    open(f::Function, command, mode::AbstractString="r", stdio=devnull)

Similar to `open(command, mode, stdio)`, but calls `f(stream)` on the resulting process
stream, then closes the input stream and waits for the process to complete.
Returns the value returned by `f`.
"""
function open(f::Function, cmds::AbstractCmd, args...)
    P = open(cmds, args...)
    ret = try
        f(P)
    catch e
        kill(P)
        rethrow(e)
    finally
        close(P.in)
    end
    success(P) || pipeline_error(P)
    return ret
end

function read(cmd::AbstractCmd)
    procs = open(cmd, "r", devnull)
    bytes = read(procs.out)
    success(procs) || pipeline_error(procs)
    return bytes
end

read(cmd::AbstractCmd, ::Type{String}) = String(read(cmd))

"""
    run(command, args...; wait::Bool = true)

Run a command object, constructed with backticks. Throws an error if anything goes wrong,
including the process exiting with a non-zero status (when `wait` is true).

If `wait` is false, the process runs asynchronously. You can later wait for it and check
its exit status by calling `success` on the returned process object.

When `wait` is false, the process' I/O streams are directed to `devnull`.
When `wait` is true, I/O streams are shared with the parent process.
Use [`pipeline`](@ref) to control I/O redirection.
"""
function run(cmds::AbstractCmd, args...; wait::Bool = true)
    if wait
        ps = _spawn(cmds, spawn_opts_inherit(args...)...)
        success(ps) || pipeline_error(ps)
    else
        ps = _spawn(cmds, spawn_opts_swallow(args...)...)
    end
    return ps
end

# some common signal numbers that are usually available on all platforms
# and might be useful as arguments to `kill` or testing against `Process.termsignal`
const SIGHUP  = 1
const SIGINT  = 2
const SIGQUIT = 3 # !windows
const SIGKILL = 9
const SIGPIPE = 13 # !windows
const SIGTERM = 15

function test_success(proc::Process)
    @assert process_exited(proc)
    if proc.exitcode < 0
        #TODO: this codepath is not currently tested
        throw(_UVError("could not start process $(string(proc.cmd))", proc.exitcode))
    end
    proc.exitcode == 0 && (proc.termsignal == 0 || proc.termsignal == SIGPIPE)
end

function success(x::Process)
    wait(x)
    return test_success(x)
end
success(procs::Vector{Process}) = mapreduce(success, &, procs)
success(procs::ProcessChain) = success(procs.processes)

"""
    success(command)

Run a command object, constructed with backticks, and tell whether it was successful (exited
with a code of 0). An exception is raised if the process cannot be started.
"""
success(cmd::AbstractCmd) = success(_spawn(cmd))

function pipeline_error(proc::Process)
    if !proc.cmd.ignorestatus
        error("failed process: ", proc, " [", proc.exitcode, "]")
    end
    nothing
end

function pipeline_error(procs::ProcessChain)
    failed = Process[]
    for p = procs.processes
        if !test_success(p) && !p.cmd.ignorestatus
            push!(failed, p)
        end
    end
    isempty(failed) && return nothing
    length(failed) == 1 && pipeline_error(failed[1])
    msg = "failed processes:"
    for proc in failed
        msg = string(msg, "\n  ", proc, " [", proc.exitcode, "]")
    end
    error(msg)
end

"""
    kill(p::Process, signum=SIGTERM)

Send a signal to a process. The default is to terminate the process.
Returns successfully if the process has already exited, but throws an
error if killing the process failed for other reasons (e.g. insufficient
permissions).
"""
function kill(p::Process, signum::Integer)
    if process_running(p)
        @assert p.handle != C_NULL
        err = ccall(:uv_process_kill, Int32, (Ptr{Cvoid}, Int32), p.handle, signum)
        if err != 0 && err != UV_ESRCH
            throw(_UVError("kill", err))
        end
    end
end
kill(ps::Vector{Process}) = foreach(kill, ps)
kill(ps::ProcessChain) = foreach(kill, ps.processes)
kill(p::Process) = kill(p, SIGTERM)

function _contains_newline(bufptr::Ptr{Cvoid}, len::Int32)
    return (ccall(:memchr, Ptr{Cvoid}, (Ptr{Cvoid},Int32,Csize_t), bufptr, '\n', len) != C_NULL)
end

## process status ##

"""
    process_running(p::Process)

Determine whether a process is currently running.
"""
process_running(s::Process) = s.exitcode == typemin(fieldtype(Process, :exitcode))
process_running(s::Vector{Process}) = any(process_running, s)
process_running(s::ProcessChain) = process_running(s.processes)

"""
    process_exited(p::Process)

Determine whether a process has exited.
"""
process_exited(s::Process) = !process_running(s)
process_exited(s::Vector{Process}) = all(process_exited, s)
process_exited(s::ProcessChain) = process_exited(s.processes)

process_signaled(s::Process) = (s.termsignal > 0)

#process_stopped (s::Process) = false #not supported by libuv. Do we need this?
#process_stop_signal(s::Process) = false #not supported by libuv. Do we need this?

function process_status(s::Process)
    process_running(s) ? "ProcessRunning" :
    process_signaled(s) ? "ProcessSignaled("*string(s.termsignal)*")" :
    #process_stopped(s) ? "ProcessStopped("*string(process_stop_signal(s))*")" :
    process_exited(s) ? "ProcessExited("*string(s.exitcode)*")" :
    error("process status error")
end

## implementation of `cmd` syntax ##

arg_gen()          = String[]
arg_gen(x::AbstractString) = String[cstr(x)]
arg_gen(cmd::Cmd)  = cmd.exec

function arg_gen(head)
    if isiterable(typeof(head))
        vals = String[]
        for x in head
            push!(vals, cstr(string(x)))
        end
        return vals
    else
        return String[cstr(string(head))]
    end
end

function arg_gen(head, tail...)
    head = arg_gen(head)
    tail = arg_gen(tail...)
    vals = String[]
    for h = head, t = tail
        push!(vals, cstr(string(h,t)))
    end
    return vals
end

function cmd_gen(parsed)
    args = String[]
    for arg in parsed
        append!(args, arg_gen(arg...))
    end
    return Cmd(args)
end

"""
    @cmd str

Similar to `cmd`, generate a `Cmd` from the `str` string which represents the shell command(s) to be executed.
The [`Cmd`](@ref) object can be run as a process and can outlive the spawning julia process (see `Cmd` for more).

# Examples
```jldoctest
julia> cm = @cmd " echo 1 "
`echo 1`

julia> run(cm)
1
Process(`echo 1`, ProcessExited(0))
```
"""
macro cmd(str)
    return :(cmd_gen($(esc(shell_parse(str, special=shell_special)[1]))))
end

wait(x::Process)      = if !process_exited(x); stream_wait(x, x.exitnotify); end
wait(x::ProcessChain) = for p in x.processes; wait(p); end

show(io::IO, p::Process) = print(io, "Process(", p.cmd, ", ", process_status(p), ")")

# allow the elements of the Cmd to be accessed as an array or iterator
for f in (:length, :firstindex, :lastindex, :keys, :first, :last, :iterate)
    @eval $f(cmd::Cmd) = $f(cmd.exec)
end
eltype(::Type{Cmd}) = eltype(fieldtype(Cmd, :exec))
for f in (:iterate, :getindex)
    @eval $f(cmd::Cmd, i) = $f(cmd.exec, i)
end

max_stars_repo_name: pochoi/julia

max_issues_repo_name: amilsted/julia

max_forks_repo_name: amilsted/julia


text:
#================================================================

Common typedefs and functions used throughout the LSH module.

================================================================#

#========================
Global variables and constants
========================#

# available_similarities: a set of all of the similarity functions that have
# been associated with hash functions via the register_similarity! macro.
const available_similarities = Set()

# Defaults to use for common arguments
const DEFAULT_N_HASHES = 1
const DEFAULT_DTYPE = Float32
const DEFAULT_RESIZE_POW2 = false

#========================
Abstract typedefs
========================#

abstract type LSHFunction end

@doc """
    abstract type SymmetricLSHFunction <: LSHFunction end

A symmetric locality-sensitive hashing function. A `SymmetricLSHFunction` uses the same hash function to insert items into a hash table as well as query the table for collisions. If `hashfn` is a `SymmetricLSHFunction`, you can compute the hash for an input `x` as `hashfn(x)`.

See also: [`AsymmetricLSHFunction`](@ref)
"""
abstract type SymmetricLSHFunction <: LSHFunction end

@doc """
    abstract type AsymmetricLSHFunction <: LSHFunction end

An asymmetric locality-sensitive hashing function. An `AsymmetricLSHFunction` uses one hash function to insert items into a hash table, and a different hash function to query the table for collisions. If `hashfn` is an `AsymmetricLSHFunction`, you can compute the indexing hash for an input `x` with `index_hash(hashfn,x)`, and the querying hash with `query_hash(hashfn,x)`.

See also: [`SymmetricLSHFunction`](@ref)
"""
abstract type AsymmetricLSHFunction <: LSHFunction end

function Base.show(io::IO, hashfn::LSHFunction)
    args = "n_hashes = $(n_hashes(hashfn))"
    args = args * ", similarity = $(similarity(hashfn))"
    args = args * ", hashtype = $(hashtype(hashfn))"
    msg = "$(typeof(hashfn))($args)"
    print(io, msg)
end

#========================
Similarity function API
========================#

# Value type to encode different similarity functions
struct SimilarityFunction{F} end
SimilarityFunction(sim) = SimilarityFunction{sim}()

#========================
LSHFunction API
========================#

macro register_similarity! end
function LSHFunction end
function lsh_family end

@doc """
    collision_probability(
        hashfn::H,
        sim;
        n_hashes::Union{Symbol,Integer}=:auto
    ) where {H <: LSHFunction}

Compute the probability of hash collision between two inputs with similarity `sim` for an [`LSHFunction`](@ref) of type `H`. This function returns the probability that `n_hashes` hashes simultaneously collide.

# Arguments
- `hashfn::LSHFunction`: the `LSHFunction` for which we want to compute the probability of collision.
- `sim`: a similarity (or vector of similarities), computed using the similarity function returned by `similarity(hashfn)`.

# Keyword arguments
- `n_hashes::Union{Symbol,Integer}` (default: `:auto`): the number of hash functions to use to compute the probability of collision. If the probability that a single hash collides is `p`, then the probability that `n_hashes` hashes simultaneously collide is `p^n_hashes`. As a result,

  ```
  collision_probability(hashfn, sim; n_hashes=N)
  ```

  is the same as

  ```
  collision_probability(hashfn, sim; n_hashes=1).^N
  ```

  If `n_hashes = :auto` then this function will select the number of hashes to be `n_hashes(hashfn)` (using the [`n_hashes`](@ref) function from the [`LSHFunction`](@ref) API).

# Examples
The probability that a single MinHash hash function causes a hash collision between inputs `A` and `B` is equal to `jaccard(A,B)`:

```jldoctest; setup = :(using LSHFunctions)
julia> hashfn = MinHash();

julia> A = Set(["a", "b", "c"]);

julia> B = Set(["b", "c", "d"]);

julia> jaccard(A,B)
0.5

julia> collision_probability(hashfn, jaccard(A,B); n_hashes=1)
0.5
```

If our [`MinHash`](@ref) struct keeps track of `N` hash functions simultaneously, then the probability of collision is `jaccard(A,B)^N`:

```jldoctest; setup = :(using LSHFunctions)
julia> hashfn = MinHash(10);

julia> A = Set(["a", "b", "c"]);

julia> B = Set(["b", "c", "d"]);

julia> collision_probability(hashfn, jaccard(A,B)) ==
       collision_probability(hashfn, jaccard(A,B); n_hashes=10) ==
       collision_probability(hashfn, jaccard(A,B); n_hashes=1)^10
true
```

See also: [`n_hashes`](@ref), [`similarity`](@ref)
"""
@generated function collision_probability(hashfn::LSHFunction, sim;
                                          n_hashes::Union{Symbol,Integer} = :auto)

    error_msg = :("n_hashes must be :auto or a positive Integer" |>
                  ErrorException |>
                  throw)

    n_hashes = begin
        if n_hashes <: Symbol
            quote
                if n_hashes != :auto
                    $error_msg
                end

                n_hashes = _n_hashes(hashfn)
            end
        else
            quote
                if n_hashes ≤ 0
                    $error_msg
                end
                nh = n_hashes
            end
        end
    end

    quote
        $n_hashes
        single_hash_collision_probability(hashfn, sim).^n_hashes
    end
end

@doc """
    collision_probability(hashfn::LSHFunction, x, y;
                          n_hashes::Union{Symbol,Integer} = :auto)

Computes the probability of a hash collision between two inputs `x` and `y` for a given hash function `hashfn`. This is the same as calling

    collision_probability(hashfn, similarity(hashfn)(x,y); n_hashes=n_hashes)

# Examples
The following snippet computes the probability of collision between two sets `A` and `B` for a single MinHash. For MinHash, this probability is just equal to the Jaccard similarity between `A` and `B`.

```jldoctest; setup = :(using LSHFunctions)
julia> hashfn = MinHash();

julia> A = Set(["a", "b", "c"]);

julia> B = Set(["a", "b", "c"]);

julia> similarity(hashfn) == jaccard
true

julia> collision_probability(hashfn, A, B) ==
       collision_probability(hashfn, jaccard(A,B)) ==
       jaccard(A,B)
true
```

We can use the `n_hashes` argument to specify the probability that `n_hashes` MinHash hash functions simultaneously collide. If left unspecified, then we'll simply use `n_hashes(hashfn)` as the number of hash functions:

```jldoctest; setup = :(using LSHFunctions)
julia> hashfn = MinHash(10);

julia> A = Set(["a", "b", "c"]);

julia> B = Set(["a", "b", "c"]);

julia> collision_probability(hashfn, A, B) ==
       collision_probability(hashfn, A, B; n_hashes=10) ==
       collision_probability(hashfn, A, B; n_hashes=1)^10
true
```
"""
collision_probability(hashfn::LSHFunction, A, B; kws...) =
    collision_probability(hashfn, similarity(hashfn)(A,B); kws...)

#=
The following functions must be defined for all LSHFunction subtypes
=#

@doc """
    similarity(hashfn::LSHFunction)

Returns the similarity function that `hashfn` hashes on.

# Arguments
- `hashfn::AbstractLSHFunction`: the hash function whose similarity we would like to retrieve.

# Examples
```jldoctest; setup = :(using LSHFunctions)
julia> hashfn = LSHFunction(cossim);

julia> similarity(hashfn) == cossim
true

julia> hashfn = LSHFunction(ℓ1);

julia> similarity(hashfn) == ℓ1
true
```
"""
function similarity end

@doc """
    hashtype(hashfn::LSHFunction)

Returns the type of hash generated by a hash function.

# Examples
```jldoctest; setup = :(using LSHFunctions)
julia> hashfn = LSHFunction(cossim);

julia> hashtype(hashfn)
Bool

julia> hashfn = LSHFunction(ℓ1);

julia> hashtype(hashfn)
Int32
```
"""
function hashtype end

@doc """
    n_hashes(hashfn::LSHFunction)

Return the number of hashes computed by `hashfn`.

# Examples
```jldoctest; setup = :(using LSHFunctions)
julia> hashfn = SimHash();

julia> n_hashes(hashfn)
$(DEFAULT_N_HASHES)

julia> hashfn = SimHash(12);

julia> n_hashes(hashfn)
12

julia> hashes = hashfn(rand(25));

julia> length(hashes)
12
```
"""
function n_hashes end

# Alias for n_hashes that's occasionally useful when we need to process
# variables that are named n_hashes
const _n_hashes = n_hashes

# The function
#
#       single_hash_collision_probability(hashfn::H, sim)
#
# must be implemented for every subtype H of LSHFunction. Note that users don't
# access this function directly; instead, they use the collision_probability
# function exported by the LSH API.
function single_hash_collision_probability end

#========================
SymmetricLSHFunction API
========================#

@doc """
    index_hash(hashfn::SymmetricLSHFunction, x)

Identical to calling `hashfn(x)`.

See also: [`query_hash`](@ref), [`SymmetricLSHFunction`](@ref)
"""
index_hash(hashfn::SymmetricLSHFunction, x) = hashfn(x)

@doc """
    query_hash(hashfn::SymmetricLSHFunction, x)

Identical to calling `hashfn(x)`.

See also: [`index_hash`](@ref), [`SymmetricLSHFunction`](@ref)
"""
query_hash(hashfn::SymmetricLSHFunction, x) = hashfn(x)

#=
The following functions must be defined for all SymmetricLSHFunction subtypes
=#
function (::SymmetricLSHFunction) end

#========================
Abstract typedefs
========================#

#=
The following functions must be defined for all AsymmetricLSHFunction subtypes
=#

@doc """
    index_hash(hashfn::AsymmetricLSHFunction, x)

Computes the indexing hash (the hash used to insert items into the hash table) for an `AsymmetricLSHFunction` with input `x`.

See also: [`query_hash`](@ref), [`AsymmetricLSHFunction`](@ref)
"""
function index_hash end

@doc """
    query_hash(hashfn::AsymmetricLSHFunction, x)

Computes the querying hash (the hash used to query for items in the hash table) for an `AsymmetricLSHFunction` with input `x`.

See also: [`index_hash`](@ref), [`AsymmetricLSHFunction`](@ref)
"""
function query_hash end

#========================
Documentation utilities
========================#

available_similarities_as_strings() = available_similarities .|> string |> sort

### Docstring generators for common keyword arguments
N_HASHES_DOCSTR(; default = DEFAULT_N_HASHES) = """
`n_hashes::Integer` (default: `$(default)`): the number of hash functions to generate."""

DTYPE_DOCSTR(hashfn; default = DEFAULT_DTYPE) = """
`dtype::Type` (default: `$(default)`): the data type to use in the $(hashfn) internals. For performance reasons you should pick `dtype` to match the type of the data you're hashing."""

RESIZE_POW2_DOCSTR(hashfn; default = DEFAULT_RESIZE_POW2) = """
`resize_pow2::Bool` (default: `$(default)`): affects the way in which the returned `$(hashfn)` resizes to hash inputs of different sizes. If you think you'll be hashing inputs of many different sizes, it's more efficient to set `resize_pow2 = true`."""


text:
module FragProcessing

include("gp_fm_reg.jl")
include("processing.jl")
end # module


text:
function get_method_count(cur_package::Module)

  method_count = 0

  for cur_name in cur_package.names(cur_package, true)

    isdefined(cur_package, cur_name) || continue

    cur_obj = getfield(cur_package, cur_name)

    isa(cur_obj, Base.Callable) || continue

    method_count += length(methods(cur_obj))

  end

  method_count

end


text:
#---------------------------------------------------------
# # Radiocarbon with OCIM1
#---------------------------------------------------------

#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`radiocarbon_OCIM.ipynb`](@__NBVIEWER_ROOT_URL__/examples/generated/radiocarbon_OCIM.ipynb)

#---------------------------------------------------------
# ## Setup the Radiocarbon model
#---------------------------------------------------------

# We will use the same model of the radiocarbon that is in the trcaer-transport-operator example.

# Radiocarbon, ¹⁴C, is produced by cosmic rays in the lower stratosphere and upper troposphere.
# It quickly reacts with oxygen to produce ¹⁴CO₂, which is then mixed throughout the troposphere and enters the ocean through air–sea gas exchange.
# Because the halflife of radiocarbon is only 5730 years a significant amount of decay can occur before the dissolved inorganic radiocarbon (DI¹⁴C) can mix uniformally throughout the ocean.
# As such the ¹⁴C serves as a tracer label for water that was recently in contact with the atmosphere.

# ### Tracer Equation

# Mathematically, the ¹⁴C tracer concentration, denoted $R$ (for Radiocarbon), satisfies the following tracer equation:
#
# $$\frac{\partial R}{\partial t} + \nabla \cdot \left[ \boldsymbol{u} - \mathbf{K} \cdot \nabla \right] R = \Lambda(R_\mathsf{atm} - R) - R / \tau,$$
#
# where $\Lambda(R_\mathsf{atm} - R)$ represents the air–sea exchanges and $R / \tau$ the radioactive decay rate.
# ($\tau$ is the radioactive decay timescale.)

# The discretized tracer is thus given by
#
# $$\frac{\partial \boldsymbol{R}}{\partial t} + \mathbf{T} \, \boldsymbol{R} = \mathbf{\Lambda}(R_\mathsf{atm} - \boldsymbol{R}) - \boldsymbol{R} / \tau.$$

# ### Translation to AIBECS Code

# We will perform an idealized radiocarbon simulation in our model and use the ocean circulation defined earlier using AIBECS.
# In this model we prescribe the atmospheric concentration, $R_\mathsf{atm}$, to be simply equal to 1.
# (We do not specify its unit or its specific value because it is not important for determining the age of a water parcel — only the decay rate does.)

# To use AIBECS, one must put the equations into the generic form of
#
# $$\frac{\partial \boldsymbol{x}}{\partial t} + \mathbf{T}(\boldsymbol{p}) \, \boldsymbol{x} = \boldsymbol{G}(\boldsymbol{x}, \boldsymbol{p}),$$
#
# where $\boldsymbol{x}$ is the state vector, $\boldsymbol{p}$ is the vector of model parameters, $\mathbf{T}(\boldsymbol{p})$ is the transport operator, and $\boldsymbol{G}(\boldsymbol{x}, \boldsymbol{p})$ is the local sources minus sinks.

# In our radiocarbon-model context, with $\boldsymbol{x} = \boldsymbol{R}$, we have that
#
# $$\boldsymbol{G}(\boldsymbol{x}, \boldsymbol{p}) = \mathbf{\Lambda}(R_\mathsf{atm} - \boldsymbol{x}) - \boldsymbol{x} / \tau.$$
#
# Hence, we must create `T(p)` and `G(x,p)` to give AIBECS the means to simulate the tracer distribution and/or its evolution in time.

# Like for any models using AIBECS, we start by telling Julia just that:

using AIBECS

# We then load the OCIM1 model (the grid and the transport) via

wet3D, grd, T = OCIM1.load() ;

# #### Sources and Sinks

# The local sources and sinks are thus simply given by

function G(x, p)
    τ, Ratm = p.τ, p.Ratm
    return Λ(Ratm .- x, p) - x / τ
end

# where `τ` is the decay rate timescale and `Ratm` is the atmospheric concentration of radiocarbon.

# We must define the air–sea exchange rate, `Λ(x,p)`, which requires us to define which boxes are located at the surface first.
# This is done, e.g., via

iwet = findall(vec(wet3D))
surface_boxes = grd.depth_3D[iwet] .== grd.depth[1]

# The air–sea exchange rate is then given by

function Λ(x, p)
    λ, h = p.λ, p.h
    return λ / h * surface_boxes .* x
end

# where `λ` is the piston velocity and `h` is the height of the top layer of the model grid.

# #### Parameters

# For the air–sea gas exchange, we use a constant piston velocity $\lambda$ of 50m / 10years, which will happen in the top layer, of height given by, well, the height of the top layer.
# And for the radioactive decay we use a timescale $\tau$ of 5730/log(2) years.
# We define these as parameters using the dedicated API from the AIBECS:

t = empty_parameter_table()                   # initialize an empty table of parameters
add_parameter!(t, :τ, 5730u"yr"/log(2)) # radioactive decay e-folding timescale
add_parameter!(t, :λ, 50u"m" / 10u"yr") # piston velocity
add_parameter!(t, :h, grd.δdepth[1])    # height of top layer
add_parameter!(t, :Ratm, 1.0u"mol/m^3") # atmospheric concentration
t

# shows the parameters that you just created.

# We now generate a new object to contain all these parameters via

initialize_Parameters_type(t, "C14_OCIM_parameters") # creates the type for parameters
p = C14_OCIM_parameters()                            # creates the parameters object

# #### Generate the state function and its Jacobian

# The last step for the setup is for AIBECS to create $\boldsymbol{F}(\boldsymbol{x}, \boldsymbol{p}) = \boldsymbol{G}(\boldsymbol{x}, \boldsymbol{p}) - \mathbf{T}(\boldsymbol{p}) \, \boldsymbol{x}$, which defines the rate of change of the state, $\boldsymbol{x}$.
# This is done via

F, ∇ₓF = state_function_and_Jacobian(p -> T, G) # generates the state function (and its Jacobian!)
nb = length(iwet)
x = zeros(nb)
F(x,p)

#md # !!! note
#md #     Here, AIBECS has automatically created `∇ₓF`, i.e., $\nabla_{\boldsymbol{x}}\boldsymbol{F}(\boldsymbol{x}, \boldsymbol{p})$, which is the Jacobian of the system.
#md #     This Jacobian will be useful in the simulations below.

#nb # > **Note**
#nb # > Here, AIBECS has automatically created `∇ₓF`, i.e., $\nabla_{\boldsymbol{x}}\boldsymbol{F}(\boldsymbol{x}, \boldsymbol{p})$, which is the Jacobian of the system.
#nb # > This Jacobian will be useful in the simulations below.

# That's it!
# Your model is entirely setup and ready to be used for simulations!

#---------------------------------------------------------
# ## Run the Radiocarbon model
#---------------------------------------------------------

# ### Compute the steady-state

# Because this model, embedded in the OCIM1 circulation, is much bigger than when embedded in the 5-box shoebox model, it would take quite some computational resournces (in time) to implicitly time-step.
# Instead, we will directly solve for the steady-state.
# For that, we simply define the steady-state problem and solve it via

prob = SteadyStateProblem(F, ∇ₓF, x, p) # define the problem
R = solve(prob, CTKAlg()).u             # solve the problem

# This should take a few seconds on a laptop.
# Once the radiocarbon concentration is computed, we can convert it into the corresponding age, via

C14age = -log.(R) * p.τ * u"s" .|> u"yr"

# ### Plot the radiocarbon age

# First, we must rearrange the age from its column vector shape into the corresponding 3D field.
# We start by filling a 3D array of the same size as the grid with some `NaN`s.
# (It serves as a blank canvas that we will "color" in with our computed age.)

C14age_3D = fill(NaN, size(grd))     # creates a 3D array of NaNs

# We use the indices of wet boxes, `iwet`, to fill in the 3D array, via

C14age_3D[iwet] .= ustrip.(C14age)   # Fills the wet grid boxes with the age values

# where we have removed the unit (for plotting).
# We then pick a layer to plot, e.g., the first layer with a depth larger than 700m.
# For that, we create `iz`, which is the index of that layer:

iz = findfirst(grd.depth .> 700u"m") # aim for a depth of ~ 700 m

# We then take a horizontal 2D slice of the 3D age at index `iz`

C14age_3D_1000m_yr = C14age_3D[:,:,iz]

# Finally, we plot the radiocarbon age using Cartopy, via:

ENV["MPLBACKEND"]="qt5agg"
using PyPlot, PyCall
clf()
ccrs = pyimport("cartopy.crs")
ax = subplot(projection=ccrs.EqualEarth(central_longitude=-155.0))
ax.coastlines()
lon_cyc = ustrip.([grd.lon; grd.lon[1] + 360u"°"]) # making it cyclic for Cartopy
age_cyc = hcat(C14age_3D_1000m_yr, C14age_3D_1000m_yr[:,1])
p = contourf(lon_cyc, ustrip.(grd.lat), age_cyc, levels=0:100:2000, transform=ccrs.PlateCarree(), zorder=-1)
colorbar(p, orientation="horizontal")
title("¹⁴C age at $(string(round(typeof(1u"m"),grd.depth[iz]))) depth using the OCIM1 circulation")
gcf() # gets the current figure to display