text:
# Parts of this code were taken / derived from NetworkX. See LICENSE for
# licensing details.

"""Creates a complete graph with `n` vertices. A complete graph has edges
connecting each pair of vertices.
"""
function CompleteGraph(n::Integer)
    g = Graph(n)
    for i = 1:n, j=1:n
        if i < j
            add_edge!(g, Edge(i,j))
        end
    end
    return g
end

"""Creates a complete bipartite graph with `n1+n2` vertices. It has edges
connecting each pair of vertices in the two sets.
"""
function CompleteBipartiteGraph(n1::Integer, n2::Integer)
    g = Graph(n1+n2)
    for i = 1:n1, j=n1+1:n1+n2
        add_edge!(g, Edge(i,j))
    end
    return g
end

"""Creates a complete digraph with `n` vertices. A complete digraph has edges
connecting each pair of vertices (both an ingoing and outgoing edge).
"""
function CompleteDiGraph(n::Integer)
    g = DiGraph(n)
    for i = 1:n, j=1:n
        if i != j
            add_edge!(g, Edge(i,j))
        end
    end
    return g
end

"""Creates a star graph with `n` vertices. A star graph has a central vertex
with edges to each other vertex.
"""
function StarGraph(n::Integer)
    g = Graph(n)
    for i = 2:n
        add_edge!(g, Edge(1,i))
    end
    return g
end

"""Creates a star digraph with `n` vertices. A star digraph has a central
vertex with directed edges to every other vertex.
"""
function StarDiGraph(n::Integer)
    g = DiGraph(n)
    for i = 2:n
        add_edge!(g, Edge(1,i))
    end
    return g
end

"""Creates a path graph with `n` vertices. A path graph connects each
successive vertex by a single edge."""
function PathGraph(n::Integer)
    g = Graph(n)
    for i = 2:n
        add_edge!(g, Edge(i-1, i))
    end
    return g
end

"""Creates a path digraph with `n` vertices. A path graph connects each
successive vertex by a single directed edge."""
function PathDiGraph(n::Integer)
    g = DiGraph(n)
    for i = 2:n
        add_edge!(g, Edge(i-1, i))
    end
    return g
end

"""Creates a cycle graph with `n` vertices. A cycle graph is a closed path graph.
"""
function CycleGraph(n::Integer)
    g = Graph(n)
    for i = 1:n-1
        add_edge!(g, Edge(i, i+1))
    end
    add_edge!(g, Edge(n, 1))
    return g
end

"""Creates a cycle digraph with `n` vertices. A cycle digraph is a closed path digraph.
"""
function CycleDiGraph(n::Integer)
    g = DiGraph(n)
    for i = 1:n-1
        add_edge!(g, Edge(i, i+1))
    end
    add_edge!(g, Edge(n, 1))
    return g
end

"""Creates a wheel graph with `n` vertices. A wheel graph is a star graph with
the outer vertices connected via a closed path graph.
"""
function WheelGraph(n::Integer)
    g = StarGraph(n)
    for i = 3:n
        add_edge!(g, Edge(i-1, i))
    end
    if n != 2
        add_edge!(g, Edge(n, 2))
    end
    return g
end

"""Creates a wheel digraph with `n` vertices. A wheel graph is a star digraph
with the outer vertices connected via a closed path graph.
"""
function WheelDiGraph(n::Integer)
    g = StarDiGraph(n)
    for i = 3:n
        add_edge!(g, Edge(i-1, i))
    end
    if n != 2
        add_edge!(g, Edge(n, 2))
    end
    return g
end

"""
    Grid{T<:Integer}(dims::AbstractVector{T}; periodic=false)

Creates a `d`-dimensional cubic lattice, with `d=length(dims)` and length  `dims[i]` in dimension `i`.
If `periodic=true` the resulting lattice will have periodic boundary condition in each dimension.
"""
function Grid{T<:Integer}(dims::AbstractVector{T}; periodic=false)
    func = periodic ? CycleGraph : PathGraph
    g = func(dims[1])
    for d in dims[2:end]
        g = cartesian_product(func(d), g)
    end
    return g
end

"""create a binary tree with k-levels vertices are numbered 1:2^levels-1"""
function BinaryTree(levels::Int)
    g = Graph(2^levels-1)
    for i in 0:levels-2
        for j in 2^i:2^(i+1)-1
            add_edge!(g, j, 2j)
            add_edge!(g, j, 2j+1)
        end
    end
    return g
end

"""create a double complete binary tree with k-levels
used as an example for spectral clustering by Guattery and Miller 1998."""
function DoubleBinaryTree(levels::Int)
    gl = BinaryTree(levels)
    gr = BinaryTree(levels)
    g = blkdiag(gl, gr)
    add_edge!(g,1, nv(gl)+1)
    return g
end

"""The Roach Graph from Guattery and Miller 1998"""
function RoachGraph(k::Int)
    dipole = CompleteGraph(2)
    nopole = Graph(2)
    antannae = crosspath(k, nopole)
    body = crosspath(k,dipole)
    roach = blkdiag(antannae, body)
    add_edge!(roach, nv(antannae)-1, nv(antannae)+1)
    add_edge!(roach, nv(antannae), nv(antannae)+2)
    return roach
end

"""This function generates `n` connected k-cliques """
function CliqueGraph(k::Integer, n::Integer)
    g = Graph(k*n)
    for c=1:n
        for i=(c-1)*k+1:c*k-1, j=i+1:c*k
            add_edge!(g, i, j)
        end
    end
    for i=1:n-1
        add_edge!(g, (i-1)*k+1, i*k+1)
    end
    add_edge!(g, 1, (n-1)*k+1)
    return g
end


text:
function selectionsort!(a::Vector, isless=isless)
    for i = 1:length(a)-1
        minindex = i
        for j = length(a):-1:i+1
            if isless(a[j], a[minindex])
                minindex = j
            end
        end
        swap!(a, i, minindex)
    end
    a
end

ext: jl


text:
type ArrayGenome{T,N} <: IndexableGenome{T}
    array::Array{T,N}
end

function ArrayGenome{T,N}(x0::Array{T,N})
    ArrayGenome{T,N}(x0)
end

function Base.show{T,N}(io::IO, genome::ArrayGenome{T,N})
    print(io, "ArrayGenome(", genome.array, ")")
end

Base.eltype(genome::ArrayGenome) = eltype(genome.array)
Base.size(genome::ArrayGenome) = size(genome.array)
Base.size(genome::ArrayGenome, i) = size(genome.array, i)
Base.length(genome::ArrayGenome) = length(genome.array)
Base.getindex(genome::ArrayGenome, i) = genome.array[i]
Base.setindex!(genome::ArrayGenome, v, i) = genome.array[i] = v

function crossover{T,N}(g1::ArrayGenome{T,N}, g2::ArrayGenome{T,N}, ::Random = Random())
    @assert length(g1) == length(g2)
    new_array = similar(g1.array)
    @inbounds @simd for i = 1:length(g1)
        new_array[i] = rand(Bool) ? g1[i] : g2[i]
    end
    ArrayGenome(new_array)::ArrayGenome{T,N}
end


text:
module HierarchialPerformanceTest

using HypothesisTests
using Statistics

const PerfMatrix = Array{Float64,2}

function hptRankSumTest(x::AbstractVector{S}, y::AbstractVector{T}, num_runs::Int) where {S<:Real,T<:Real}
    # Custom threshold as specified by the paper
    if num_runs < 12
        ExactMannWhitneyUTest(x, y)
    else
        ApproximateMannWhitneyUTest(x, y)
    end
end

function hptSignedRankTest(x::AbstractVector{T}, num_tests::Int) where T<:Real
    # Custom threshold as specified by the paper
    if num_tests < 25
        ExactSignedRankTest(x)
    else
        ApproximateSignedRankTest(x)
    end
end

function hpt(a::PerfMatrix, b::PerfMatrix; tail=:right)
    # Enforce identical and sane matrix sizes
    if size(a) != size(b)
        error("mismatching number of benchmarks/runs")
    elseif length(a) == 0
        error("no benchmarks given")
    end

    num_tests = size(a, 1) # n
    num_runs = size(a, 2) # m

    # Max significance for the *null* hypothesis, not the alternative
    maxNullSignificance = if (num_runs >= 5) 0.05 else 0.10 end

    # Calculate performance delta per benchmark
    test_deltas = Array{Float64}(undef, num_tests) # dT
    for test = 1:num_tests
        aRuns = a[test,:]
        bRuns = b[test,:]

        # Use SHT to determine which side is better
        shtP = pvalue(hptRankSumTest(aRuns, bRuns, num_runs), tail=tail)

        # Check is inverted because shtP = probability of null hypothesis (equal)
        test_deltas[test] = if (shtP < maxNullSignificance)
            median(aRuns) - median(bRuns)
        else
            0
        end
    end

    # Evaluate general performance using the selected SHT
    1 - pvalue(hptSignedRankTest(test_deltas, num_tests), tail=tail)
end

function hptSpeedup(a::PerfMatrix, b::PerfMatrix, minConfidence::Float64 = 0.95, digits::Int = 3, increment::Int = 1)
    # Calculate an extra digit internally to facilitate rounding
    baseSpeedup = 10 ^ digits
    # Use int for storing speedup to mitigate floating-point precision loss
    intSpeedup = baseSpeedup
    speedup = 1.0
    tail = :right

    # Invert if slower
    if hpt(a, b) < minConfidence
        increment *= -1
        tail = :left
    end

    while true
        confidence = hpt(a / speedup, b, tail=tail)
        if confidence >= minConfidence
            # Increment int and recalculate float to mitigate floating-point precision loss
            intSpeedup += increment
            speedup = intSpeedup / baseSpeedup
        else
            break
        end
    end

    speedup
end

function calcScore(results::PerfMatrix, reference::PerfMatrix, refScore::Int = 1000)
    # 95%-confidence speedup
    # We swap a and b because HPT considers higher values to be better, but
    # since our values are times rather than performance, lower is better
    round(Int, hptSpeedup(reference, results) * refScore)
end

export PerfMatrix, hpt, hptSpeedup, calcScore

end


text:
function compute(n::Int)::Int
    i = 1
    while true
        for j ∈ 10 ^ i:10 ^ (i + 1) ÷ n - 1
            j_digits = sort(collect(string(j)))
            is_permuted = true
            for k ∈ 2:n
                k_digits = sort(collect(string(j * k)))
                if j_digits != k_digits
                    is_permuted = false
                    break
                end
            end
            if is_permuted
                return j
            end
        end
    i += 1
    end
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
# Parameters
#	N::Int
#	- the problem size
# Return Values
#	YOR::Array{Array{Array{SparseMatrixCSC, 1}, 1}, 1} (Young's Orthogonal Representations)
#	- YOR[n][p][k] is Young's Orthogonal Representation for the Adjacent Transposition (K, K + 1) for the pth Partition of n
#	PT::Array{Array{Array{Int, 1}, 1}, 1} (Partition Tree)
#	- output1 from partition_tree()
function yor(N::Int)
	P, WI = partitions(N)
	PT = partition_tree(N, P, WI)
	YS = ys_symbols(N, P, PT)
	YOR = Array(Array{Array{SparseMatrixCSC, 1}, 1}, N)
	for n = 1:N
		YSn = YS[n]
		Pn = P[n]
		Pn_L = length(Pn)
		YORn = Array(Array{SparseMatrixCSC, 1}, Pn_L)
		for p = 1:Pn_L
			YSnp = YSn[p]
			Pnp = Pn[p]
			R = length(Pnp)
			YORn[p] = yor_p(n, R, YSnp)
		end
		YOR[n] = YORn
	end
	return YOR, PT
end

# Parameters:
#	N::Int
#	- the size of P
#	R::Int
#	- the length of P
#	P::Array{Int, 1}
#	- P is a Partition of N
# Return Values:
#	YORp::Array{SparseMatrixCSC, 1}
#	- YORp[k] is the Young's Orthogonal Representation of P for the Adjacent Transposistion (k, k + 1)
function yor_p(N::Int, R::Int, YSymbols::Array{Array{Int, 1}, 1})
	YORp = Array(SparseMatrixCSC, N - 1)
	L = length(YSymbols)
	DA, IA = ys_information(N, R, YSymbols)
	for k = 1:(N - 1)
		YORp[k] = yor_pk(DA, IA, L, k)
	end
	return YORp
end

#Parameters
#	DA::Array{Int, 2}
#	- output1 from ys_information()
#	IA::Array{Int, 2}
#	- outpu2 from ys_information()
#	L::Int
#	- the number of Yamanouchi Symbols for the current Partition
#	K::Int
#	- represents the Adjacent Transposition (K, K + 1)
#Return Values
#	YORpk::SparseMatrixCSC
#	- Young's Orthongal Representation of the Adjacent Transposition (K, K + 1) for the Partition p
function yor_pk(DA::Array{Int, 2}, IA::Array{Int, 2}, L::Int, k::Int)
	n = L
	for i = 1:L
		if IA[i,k] != 0
			n += 1
		end
	end
	colptr = Array(Int, L + 1)
	colptr[1] = 1
	rowval = Array(Int, n)
	nzval = Array(Float64, n)
	cp = 1
	for i = 1:L
		D = 1/DA[i, k]
		I = IA[i,k]
		if I == 0
			colptr[i + 1] = colptr[i] + 1
			rowval[cp] = i
			nzval[cp] = D
			cp += 1
		else
			colptr[i + 1] = colptr[i] + 2
			if I > i
				rowval[cp] = i
				nzval[cp] = D
				cp += 1
				rowval[cp] = I
				nzval[cp] = sqrt(1 - D * D)
				cp += 1
			else
				rowval[cp] = I
				nzval[cp] = sqrt(1 - D * D)
				cp += 1
				rowval[cp] = i
				nzval[cp] = D
				cp += 1
			end
		end
	end
	YORpk = SparseMatrixCSC(L, L, colptr, rowval, nzval)
	return YORpk
end


text:
import CUDAnative
using CuArrays

# const NON_DISPATCHED_OPS = [log, exp, sqrt, ^, ones]
# const CUDA_NATIVE_OPS = Dict{Function,Function}(op => op for op in NON_DISPATCHED_OPS)
const CUDA_NATIVE_OPS = Dict{Function,Function}()

CUDA_NATIVE_OPS[log] = CUDAnative.log
CUDA_NATIVE_OPS[exp] = CUDAnative.exp
CUDA_NATIVE_OPS[sqrt] = CUDAnative.sqrt
CUDA_NATIVE_OPS[^] = CUDAnative.pow
CUDA_NATIVE_OPS[ones] = CUDAnative.ones

to_device(device::GPU, x) = cu(x)
device_op(device::GPU, op) = get(CUDA_NATIVE_OPS, op, op)

# import Espresso: rewrite, rewrite_all

# const CUDA_NATIVE_RULES = [
#     :($log.(x)) => :(CUDAnative.log.(x)),
#     :($exp.(x)) => :(CUDAnative.exp.(x)),
#     :($sqrt.(x)) => :(CUDAnative.sqrt.(x)),
#     :($(^).(x, n)) => :(CUDAnative.pow.(x, Float32(n))),
#     :($ones(n)) => :(CuArray(ones(Float32, n))),
#     # :(transpose(x)) => :(permutedims(x, (2,1))),  -- seems to cauase segfault in complex cases
# ]

# # note: rewriting function objects to symbolic names
# to_cuda_ex(ex::Expr) = rewrite_all(ex, CUDA_NATIVE_RULES; phs=[:x, :n])

# with_cuda() = false
# has_cuda_inputs(tape::Tape) = false
# @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
#     with_cuda() = true

#     function has_cuda_inputs(tape::Tape)
#         res = false
#         for op in tape
#             if op isa Input && op.val isa CuArray
#                 res = true
#                 break
#             end
#         end
#         return res
#     end
# end


text:
@doc raw"""
    SymplecticStiefel{n, k, 𝔽} <: AbstractEmbeddedManifold{𝔽, DefaultIsometricEmbeddingType}

The symplectic Stiefel manifold consists of all
$2n × 2k, \; n \geq k$ matrices satisfying the requirement
````math
\operatorname{SpSt}(2n, 2k, ℝ)
    = \bigl\{ p ∈ ℝ^{2n × 2n} \, \big| \, p^TQ_{2n}p = Q_{2k} \bigr\},
````
where
````math
Q_{2n} =
\begin{bmatrix}
  0_n & I_n \\
 -I_n & 0_n
\end{bmatrix}.
````

The symplectic Stiefel tangent space at ``p`` can be parametrized as [^Bendokat2021]
````math
    \begin{align*}
    T_p\operatorname{SpSt}(2n, 2k)
    = \{&X \in \mathbb{R}^{2n \times 2k} \;|\; p^{T}Q_{2n}X + X^{T}Q_{2n}p = 0 \}, \\
    = \{&X = pΩ + p^sB \;|\;
        Ω ∈ ℝ^{2k × 2k}, Ω^+ = -Ω, \\
        &\; p^s ∈ \operatorname{SpSt}(2n, 2(n- k)), B ∈ ℝ^{2(n-k) × 2k}, \},
    \end{align*}
````
where ``Ω \in \mathfrak{sp}(2n,F)`` is Hamiltonian and ``p^s`` means
the symplectic complement of ``p`` s.t. ``p^{+}p^{s} = 0``.

# Constructor
    SymplecticStiefel(2n::Int, 2k::Int, field::AbstractNumbers=ℝ)
        -> SymplecticStiefel{div(2n, 2), div(2k, 2), field}()

Generate the (real-valued) symplectic Stiefel manifold of ``2n \times 2k``
matrices which span a ``2k`` dimensional symplectic subspace of ``ℝ^{2n \times 2n}``.
The constructor for the [`SymplecticStiefel`](@ref) manifold accepts the even column
dimension ``2n`` and an even number of columns ``2k`` for
the real symplectic Stiefel manifold with elements ``p \in ℝ^{2n × 2k}``.
"""
struct SymplecticStiefel{n,k,𝔽} <: AbstractDecoratorManifold{𝔽} end

function SymplecticStiefel(two_n::Int, two_k::Int, field::AbstractNumbers=ℝ)
    return SymplecticStiefel{div(two_n, 2),div(two_k, 2),field}()
end

function active_traits(f, ::SymplecticStiefel, args...)
    return merge_traits(IsEmbeddedManifold(), IsDefaultMetric(RealSymplecticMetric()))
end

function ManifoldsBase.default_inverse_retraction_method(::SymplecticStiefel)
    return CayleyInverseRetraction()
end

ManifoldsBase.default_retraction_method(::SymplecticStiefel) = CayleyRetraction()

@doc raw"""
    canonical_projection(::SymplecticStiefel, p_Sp)
    canonical_projection!(::SymplecticStiefel{n,k}, p, p_Sp)

Define the canonical projection from ``\operatorname{Sp}(2n, 2n)`` onto
``\operatorname{SpSt}(2n, 2k)``, by projecting onto the first ``k`` columns
and the ``n + 1``'th onto the ``n + k``'th columns [^Bendokat2021].

It is assumed that the point ``p`` is on ``\operatorname{Sp}(2n, 2n)``.
"""
function canonical_projection(M::SymplecticStiefel{n,k}, p_Sp) where {n,k}
    p_SpSt = similar(p_Sp, (2n, 2k))
    return canonical_projection!(M, p_SpSt, p_Sp)
end

function canonical_projection!(::SymplecticStiefel{n,k}, p, p_Sp) where {n,k}
    p[:, (1:k)] .= p_Sp[:, (1:k)]
    p[:, ((k + 1):(2k))] .= p_Sp[:, ((n + 1):(n + k))]
    return p
end

@doc raw"""
    check_point(M::SymplecticStiefel, p; kwargs...)

Check whether `p` is a valid point on the [`SymplecticStiefel`](@ref),
$\operatorname{SpSt}(2n, 2k)$ manifold.
That is, the point has the right [`AbstractNumbers`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#number-system) type and $p^{+}p$ is
(approximately) the identity,
where for $A \in \mathbb{R}^{2n \times 2k}$,
$A^{+} = Q_{2k}^TA^TQ_{2n}$ is the symplectic inverse, with
````math
Q_{2n} =
\begin{bmatrix}
0_n & I_n \\
 -I_n & 0_n
\end{bmatrix}.
````
The tolerance can be set with `kwargs...` (e.g. `atol = 1.0e-14`).
"""
function check_point(M::SymplecticStiefel{n,k}, p; kwargs...) where {n,k}
    # Perform check that the matrix lives on the real symplectic manifold:
    expected_zero = norm(inv(M, p) * p - I)
    if !isapprox(expected_zero, 0; kwargs...)
        return DomainError(
            expected_zero,
            (
                "The point p does not lie on $(M) because its symplectic" *
                " inverse composed with itself is not the identity."
            ),
        )
    end
    return nothing
end

@doc raw"""
    check_vector(M::Symplectic, p, X; kwargs...)

Checks whether `X` is a valid tangent vector at `p` on the [`SymplecticStiefel`](@ref),
``\operatorname{SpSt}(2n, 2k)`` manifold. First recall the definition of the symplectic
inverse for $A \in \mathbb{R}^{2n \times 2k}$,
$A^{+} = Q_{2k}^TA^TQ_{2n}$ is the symplectic inverse, with
````math
    Q_{2n} =
    \begin{bmatrix}
    0_n & I_n \\
     -I_n & 0_n
\end{bmatrix}.
````
The we check that ``H = p^{+}X \in 𝔤_{2k}``, where ``𝔤``
is the Lie Algebra of the symplectic group ``\operatorname{Sp}(2k)``,
characterized as [^Bendokat2021],
````math
    𝔤_{2k} = \{H \in ℝ^{2k \times 2k} \;|\; H^+ = -H \}.
````
The tolerance can be set with `kwargs...` (e.g. `atol = 1.0e-14`).
"""
check_vector(::SymplecticStiefel, ::Any...)

function check_vector(M::SymplecticStiefel{n,k,field}, p, X; kwargs...) where {n,k,field}
    # From Bendokat-Zimmermann: T_pSpSt(2n, 2k) = \{p*H | H^{+} = -H  \}
    H = inv(M, p) * X  # ∈ ℝ^{2k × 2k}, should be Hamiltonian.
    H_star = inv(Symplectic(2k, field), H)
    hamiltonian_identity_norm = norm(H + H_star)

    if !isapprox(hamiltonian_identity_norm, 0; kwargs...)
        return DomainError(
            hamiltonian_identity_norm,
            (
                "The matrix X is not in the tangent space at point p of the" *
                " $(M) manifold, as p^{+}X is not a Hamiltonian matrix."
            ),
        )
    end
    return nothing
end

@doc raw"""
    exp(::SymplecticStiefel, p, X)
    exp!(M::SymplecticStiefel, q, p, X)

Compute the exponential mapping
````math
    \operatorname{exp}\colon T\operatorname{SpSt}(2n, 2k)
    \rightarrow \operatorname{SpSt}(2n, 2k)
````
at a point ``p \in  \operatorname{SpSt}(2n, 2k)``
in the direction of ``X \in T_p\operatorname{SpSt}(2n, 2k)``.

The tangent vector ``X`` can be written in the form
``X = \bar{\Omega}p`` [^Bendokat2021], with
````math
    \bar{\Omega} = X (p^Tp)^{-1}p^T
        + Q_{2n}p(p^Tp)^{-1}X^T(I_{2n} - Q_{2n}^Tp(p^Tp)^{-1}p^TQ_{2n})Q_{2n}
        \in ℝ^{2n \times 2n},
````
where ``Q_{2n}`` is the [`SymplecticMatrix`](@ref). Using this expression for ``X``,
the exponential mapping can be computed as
````math
    \operatorname{exp}_p(X) = \operatorname{Exp}([\bar{\Omega} - \bar{\Omega}^T])
                              \operatorname{Exp}(\bar{\Omega}^T)p,
````
where ``\operatorname{Exp}(\cdot)`` denotes the matrix exponential.

Computing the above mapping directly however, requires taking matrix exponentials
of two ``2n \times 2n`` matrices, which is computationally expensive when ``n``
increases. Therefore we instead follow [^Bendokat2021] who express the above
exponential mapping in a way which only requires taking matrix exponentials
of an ``8k \times 8k`` matrix and a ``4k \times 4k`` matrix.

To this end, first define
````math
\bar{A} = Q_{2k}p^TX(p^Tp)^{-1}Q_{2k} +
            (p^Tp)^{-1}X^T(p - Q_{2n}^Tp(p^Tp)^{-1}Q_{2k}) \in ℝ^{2k \times 2k},
````
and
````math
\bar{H} = (I_{2n} - pp^+)Q_{2n}X(p^Tp)^{-1}Q_{2k} \in ℝ^{2n \times 2k}.
````
We then let ``\bar{\Delta} = p\bar{A} + \bar{H}``, and define the matrices
````math
    γ = \left[\left(I_{2n} - \frac{1}{2}pp^+\right)\bar{\Delta} \quad
              -p \right] \in ℝ^{2n \times 4k},
````
and
````math
    λ = \left[Q_{2n}^TpQ_{2k} \quad
        \left(\bar{\Delta}^+\left(I_{2n}
              - \frac{1}{2}pp^+\right)\right)^T\right] \in ℝ^{2n \times 4k}.
````
With the above defined matrices it holds that ``\bar{\Omega} = λγ^T``.
 As a last preliminary step, concatenate ``γ`` and ``λ`` to define the matrices
``Γ = [λ \quad -γ] \in ℝ^{2n \times 8k}`` and
``Λ = [γ \quad λ] \in ℝ^{2n \times 8k}``.

With these matrix constructions done, we can compute the
exponential mapping as
````math
    \operatorname{exp}_p(X) =
        Γ \operatorname{Exp}(ΛΓ^T)
        \begin{bmatrix}
            0_{4k} \\
            I_{4k}
        \end{bmatrix}
        \operatorname{Exp}(λγ^T)
        \begin{bmatrix}
            0_{2k} \\
            I_{2k}
        \end{bmatrix}.
````
which only requires computing the matrix exponentials of
``ΛΓ^T \in ℝ^{8k \times 8k}`` and ``λγ^T \in ℝ^{4k \times 4k}``.
"""
exp(::SymplecticStiefel, p, X)

function exp!(M::SymplecticStiefel{n,k}, q, p, X) where {n,k}
    Q = SymplecticMatrix(p, X)
    pT_p = lu(p' * p) # ∈ ℝ^{2k × 2k}

    C = pT_p \ X' # ∈ ℝ^{2k × 2n}

    # Construct A-bar:
    # First A-term: Q * (p^T * C^T) * Q
    A_bar = rmul!(lmul!(Q, (p' * C')), Q)
    A_bar .+= C * p

    # Second A-term, use C-memory:
    rmul!(C, Q') # C*Q^T -> C
    C_QT = C

    # Subtract C*Q^T*p*(pT_p)^{-1}*Q:
    # A_bar is "star-skew symmetric" (A^+ = -A).
    A_bar .-= rmul!((C_QT * p) / pT_p, Q)

    # Construct H_bar:
    # First H-term: Q * (C_QT * Q)' * Q -> Q * C' * Q = Q * (X / pT_p) * Q
    H_bar = rmul!(lmul!(Q, rmul!(C_QT, Q)'), Q)
    # Subtract second H-term:
    H_bar .-= p * symplectic_inverse_times(M, p, H_bar)

    # Construct Δ_bar in H_bar-memory:
    H_bar .+= p * A_bar

    # Rename H_bar -> Δ_bar.
    Δ_bar = H_bar

    γ_1 = Δ_bar - (1 / 2) .* p * symplectic_inverse_times(M, p, Δ_bar)
    γ = [γ_1 -p] # ∈ ℝ^{2n × 4k}

    Δ_bar_star = rmul!(Q' * Δ_bar', Q)
    λ_1 = lmul!(Q', p * Q)
    λ_2 = (Δ_bar_star .- (1 / 2) .* (Δ_bar_star * p) * λ_1')'
    λ = [λ_1 λ_2] # ∈ ℝ^{2n × 4k}

    Γ = [λ -γ] # ∈ ℝ^{2n × 8k}
    Λ = [γ λ] # ∈ ℝ^{2n × 8k}

    # At last compute the (8k × 8k) and (4k × 4k) matrix exponentials:
    q .= Γ * (exp(Λ' * Γ)[:, (4k + 1):end]) * (exp(λ' * γ)[:, (2k + 1):end])
    return q
end

get_embedding(::SymplecticStiefel{n,k,ℝ}) where {n,k} = Euclidean(2n, 2k; field=ℝ)

@doc raw"""
    gradient(::SymplecticStiefel, f, p, backend::RiemannianProjectionBackend)
    gradient!(::SymplecticStiefel, f, X, p, backend::RiemannianProjectionBackend)

Compute the gradient of
``f\colon \operatorname{SpSt}(2n, 2k) \rightarrow ℝ``
at ``p \in \operatorname{SpSt}(2n, 2k)``.

This function first computes the embedding gradient ``∇f(p) \in ℝ^{2n \times 2k}`` using
the [`AbstractRiemannianDiffBackend`](@ref) in the [`RiemannianProjectionBackend`](@ref).
Then it transforms the embedding gradient to the unique tangent vector space element
``\text{grad}f(p) \in T_p\operatorname{SpSt}(2n, 2k)``
which fulfills the variational equation
````math
    g_p(\text{grad}f(p), X)
    = \text{D}f(p)
    = \langle ∇f(p), X \rangle
    \quad\forall\; X \in T_p\operatorname{SpSt}(2n, 2k).
````
The manifold gradient ``\text{grad}f(p)`` is computed from ``∇f(p)`` as
````math
    \text{grad}f(p) = ∇f(p)p^Tp + Q_{2n}p∇f(p)^TQ_{2n}p,
````
where ``Q_{2n}`` is the [`SymplecticMatrix`](@ref).
"""
function gradient(::SymplecticStiefel, f, p, backend::RiemannianProjectionBackend)
    amb_grad = _gradient(f, p, backend.diff_backend)
    Q_p = SymplecticMatrix(p, amb_grad) * p
    return amb_grad * (p' * p) .+ Q_p * (amb_grad' * Q_p)
end

function gradient!(::SymplecticStiefel, f, X, p, backend::RiemannianProjectionBackend)
    _gradient!(f, X, p, backend.diff_backend)
    Q_p = SymplecticMatrix(p, X) * p
    X .= X * (p' * p) .+ Q_p * (X' * Q_p)
    return X
end

@doc raw"""
    inner(M::SymplecticStiefel{n, k}, p, X. Y)

Compute the Riemannian inner product ``g^{\operatorname{SpSt}}`` at
``p \in \operatorname{SpSt}`` between tangent vectors ``X, X \in T_p\operatorname{SpSt}``.
Given by Proposition 3.10 in [^Bendokat2021].
````math
g^{\operatorname{SpSt}}_p(X, Y)
    = \operatorname{tr}\left(X^T\left(I_{2n} -
        \frac{1}{2}Q_{2n}^Tp(p^Tp)^{-1}p^TQ_{2n}\right)Y(p^Tp)^{-1}\right).
````
"""
function inner(::SymplecticStiefel{n,k}, p, X, Y) where {n,k}
    Q = SymplecticMatrix(p, X, Y)
    # Procompute lu(p'p) since we solve a^{-1}* 3 times
    a = lu(p' * p) # note that p'p is symmetric, thus so is its inverse c=a^{-1}
    b = Q' * p
    # we split the original trace into two one with I->(X'Yc)
    # and the other with 1/2 X'b c b' Y c
    # 1) we permute X' and Y c to c^TY^TX = a\(Y'X) (avoids a large interims matrix)
    # 2) we permute Y c up front, the center term is symmetric, so we get cY'b c b' X
    # and (b'X) again avoids a large interims matrix, so does Y'b.
    return tr(a \ (Y' * X)) - (1 / 2) * tr(a \ ((Y' * b) * (a \ (b' * X))))
end

@doc raw"""
    inv(::SymplecticStiefel{n, k}, A)
    inv!(::SymplecticStiefel{n, k}, q, p)

Compute the symplectic inverse ``A^+`` of matrix ``A ∈ ℝ^{2n × 2k}``. Given a matrix
````math
A ∈ ℝ^{2n × 2k},\quad
A =
\begin{bmatrix}
A_{1, 1} & A_{1, 2} \\
A_{2, 1} & A_{2, 2}
\end{bmatrix},\; A_{i, j} \in ℝ^{2n × 2k}
````
the symplectic inverse is defined as:
````math
A^{+} := Q_{2k}^T A^T Q_{2n},
````
where
````math
Q_{2n} =
\begin{bmatrix}
0_n & I_n \\
 -I_n & 0_n
\end{bmatrix}.
````
For any ``p \in \operatorname{SpSt}(2n, 2k)`` we have that ``p^{+}p = I_{2k}``.

The symplectic inverse of a matrix A can be expressed explicitly as:
````math
A^{+} =
\begin{bmatrix}
  A_{2, 2}^T & -A_{1, 2}^T \\[1.2mm]
 -A_{2, 1}^T &  A_{1, 1}^T
\end{bmatrix}.
````
"""
function Base.inv(M::SymplecticStiefel{n,k}, p) where {n,k}
    q = similar(p')
    return inv!(M, q, p)
end

function inv!(::SymplecticStiefel{n,k}, q, p) where {n,k}
    checkbounds(q, 1:(2k), 1:(2n))
    checkbounds(p, 1:(2n), 1:(2k))
    @inbounds for i in 1:k, j in 1:n
        q[i, j] = p[j + n, i + k]
    end
    @inbounds for i in 1:k, j in 1:n
        q[i, j + n] = -p[j, i + k]
    end
    @inbounds for i in 1:k, j in 1:n
        q[i + k, j] = -p[j + n, i]
    end
    @inbounds for i in 1:k, j in 1:n
        q[i + k, j + n] = p[j, i]
    end
    return q
end

@doc raw"""
    inverse_retract(::SymplecticStiefel, p, q, ::CayleyInverseRetraction)
    inverse_retract!(::SymplecticStiefel, q, p, X, ::CayleyInverseRetraction)

Compute the Cayley Inverse Retraction ``X = \mathcal{L}_p^{\operatorname{SpSt}}(q)``
such that the Cayley Retraction from ``p`` along ``X`` lands at ``q``, i.e.
``\mathcal{R}_p(X) = q`` [^Bendokat2021].

First, recall the definition the standard symplectic matrix
````math
Q =
\begin{bmatrix}
 0    &  I \\
-I    &  0
\end{bmatrix}
````
as well as the symplectic inverse of a matrix ``A``, ``A^{+} = Q^T A^T Q``.

For ``p, q ∈ \operatorname{SpSt}(2n, 2k, ℝ)`` then, we can define the
inverse cayley retraction as long as the following matrices exist.
````math
    U = (I + p^+ q)^{-1} \in ℝ^{2k \times 2k},
    \quad
    V = (I + q^+ p)^{-1} \in ℝ^{2k \times 2k}.
````

If that is the case, the inverse cayley retration at ``p`` applied to ``q`` is
````math
\mathcal{L}_p^{\operatorname{Sp}}(q) = 2p\bigl(V - U\bigr) + 2\bigl((p + q)U - p\bigr)
                                        ∈ T_p\operatorname{Sp}(2n).
````

[^Bendokat2021]:
    > Bendokat, Thomas and Zimmermann, Ralf:
	> The real symplectic Stiefel and Grassmann manifolds: metrics, geodesics and applications
	> arXiv preprint arXiv:2108.12447, 2021 (https://arxiv.org/abs/2108.12447)
"""
inverse_retract(::SymplecticStiefel, p, q, ::CayleyInverseRetraction)

function inverse_retract_caley!(M::SymplecticStiefel, X, p, q)
    U_inv = lu!(add_scaled_I!(symplectic_inverse_times(M, p, q), 1))
    V_inv = lu!(add_scaled_I!(symplectic_inverse_times(M, q, p), 1))

    X .= 2 .* ((p / V_inv .- p / U_inv) .+ ((p + q) / U_inv) .- p)
    return X
end

@doc raw"""
    manifold_dimension(::SymplecticStiefel{n, k})

Returns the dimension of the symplectic Stiefel manifold embedded in ``ℝ^{2n \times 2k}``,
i.e. [^Bendokat2021]
````math
    \operatorname{dim}(\operatorname{SpSt}(2n, 2k)) = (4n - 2k + 1)k.
````
"""
manifold_dimension(::SymplecticStiefel{n,k}) where {n,k} = (4n - 2k + 1) * k

@doc raw"""
    project(::SymplecticStiefel, p, A)
    project!(::SymplecticStiefel, Y, p, A)

Given a point ``p \in \operatorname{SpSt}(2n, 2k)``,
project an element ``A \in \mathbb{R}^{2n \times 2k}`` onto
the tangent space ``T_p\operatorname{SpSt}(2n, 2k)`` relative to
the euclidean metric of the embedding ``\mathbb{R}^{2n \times 2k}``.

That is, we find the element ``X \in T_p\operatorname{SpSt}(2n, 2k)``
which solves the constrained optimization problem

````math
    \operatorname{min}_{X \in \mathbb{R}^{2n \times 2k}} \frac{1}{2}||X - A||^2, \quad
    \text{s.t.}\;
    h(X)\colon= X^T Q p + p^T Q X = 0,
````
where ``h : \mathbb{R}^{2n \times 2k} \rightarrow \operatorname{skew}(2k)`` defines
the restriction of ``X`` onto the tangent space ``T_p\operatorname{SpSt}(2n, 2k)``.
"""
project(::SymplecticStiefel, p, A)

function project!(::SymplecticStiefel, Y, p, A)
    Q = SymplecticMatrix(Y, p, A)
    Q_p = Q * p

    function h(X)
        XT_Q_p = X' * Q_p
        return XT_Q_p .- XT_Q_p'
    end

    # Solve for Λ (Lagrange mutliplier):
    pT_p = p' * p  # (2k × 2k)
    Λ = sylvester(pT_p, pT_p, h(A) ./ 2)

    Y[:, :] = A .- Q_p * (Λ .- Λ')
    return Y
end

@doc raw"""
    rand(M::SymplecticStiefel; vector_at=nothing,
        hamiltonian_norm=(vector_at === nothing ? 1/2 : 1.0))

Generate a random point ``p \in \operatorname{SpSt}(2n, 2k)`` or
a random tangent vector ``X \in T_p\operatorname{SpSt}(2n, 2k)``
if `vector_at` is set to a point ``p \in \operatorname{Sp}(2n)``.

A random point on ``\operatorname{SpSt}(2n, 2k)`` is found by first generating a
random point on the symplectic manifold ``\operatorname{Sp}(2n)``,
and then projecting onto the Symplectic Stiefel manifold using the
[`canonical_projection`](@ref) ``π_{\operatorname{SpSt}(2n, 2k)}``.
That is, ``p = π_{\operatorname{SpSt}(2n, 2k)}(p_{\operatorname{Sp}})``.

To generate a random tangent vector in ``T_p\operatorname{SpSt}(2n, 2k)``
this code exploits the second tangent vector space parametrization of
[`SymplecticStiefel`](@ref), showing that any ``X \in T_p\operatorname{SpSt}(2n, 2k)``
can be written as ``X = pΩ_X + p^sB_X``.
To generate random tangent vectors at ``p`` then, this function sets ``B_X = 0``
and generates a random Hamiltonian matrix ``Ω_X \in \mathfrak{sp}(2n,F)`` with
Frobenius norm of `hamiltonian_norm` before returning ``X = pΩ_X``.
"""
function Base.rand(
    M::SymplecticStiefel{n};
    vector_at=nothing,
    hamiltonian_norm=(vector_at === nothing ? 1 / 2 : 1.0),
) where {n}
    if vector_at === nothing
        return canonical_projection(
            M,
            rand(Symplectic(2n); hamiltonian_norm=hamiltonian_norm),
        )
    else
        return random_vector(M, vector_at; hamiltonian_norm=hamiltonian_norm)
    end
end

function random_vector(
    ::SymplecticStiefel{n,k},
    p::AbstractMatrix;
    hamiltonian_norm=1.0,
) where {n,k}
    Ω = rand_hamiltonian(Symplectic(2k); frobenius_norm=hamiltonian_norm)
    return p * Ω
end

@doc raw"""
    retract(::SymplecticStiefel, p, X, ::CayleyRetraction)
    retract!(::SymplecticStiefel, q, p, X, ::CayleyRetraction)

Compute the Cayley retraction on the Symplectic Stiefel manifold, computed inplace
of `q` from `p` along `X`.

Given a point ``p \in \operatorname{SpSt}(2n, 2k)``, every tangent vector
``X \in T_p\operatorname{SpSt}(2n, 2k)`` is of the form
``X = \tilde{\Omega}p``, with
````math
    \tilde{\Omega} = \left(I_{2n} - \frac{1}{2}pp^+\right)Xp^+ -
                     pX^+\left(I_{2n} - \frac{1}{2}pp^+\right) \in ℝ^{2n \times 2n},
````
as shown in Proposition 3.5 of [^Bendokat2021].
Using this representation of ``X``, the Cayley retraction
on ``\operatorname{SpSt}(2n, 2k)`` is defined pointwise as
````math
    \mathcal{R}_p(X) = \operatorname{cay}\left(\frac{1}{2}\tilde{\Omega}\right)p.
````
The operator ``\operatorname{cay}(A) = (I - A)^{-1}(I + A)`` is the Cayley transform.

However, the computation of an ``2n \times 2n`` matrix inverse in the expression
above can be reduced down to inverting a ``2k \times 2k`` matrix due to Proposition
5.2 of [^Bendokat2021].

Let ``A = p^+X`` and ``H = X - pA``. Then an equivalent expression for the Cayley
retraction defined pointwise above is
````math
    \mathcal{R}_p(X) = -p + (H + 2p)(H^+H/4 - A/2 + I_{2k})^{-1}.
````
It is this expression we compute inplace of `q`.
"""
retract(::SymplecticStiefel, p, X, ::CayleyRetraction)

function retract_caley!(M::SymplecticStiefel, q, p, X)
    # Define intermediate matrices for later use:
    A = symplectic_inverse_times(M, p, X)

    H = X .- p * A  # Allocates (2n × 2k).

    # A = I - A/2 + H^{+}H/4:
    A .= (symplectic_inverse_times(M, H, H) ./ 4) .- (A ./ 2)
    Manifolds.add_scaled_I!(A, 1.0)

    # Reuse 'H' memory:
    H .= H .+ 2 .* p
    r = lu!(A)
    q .= (-).(p) .+ rdiv!(H, r)
    return q
end

function Base.show(io::IO, ::SymplecticStiefel{n,k,𝔽}) where {n,k,𝔽}
    return print(io, "SymplecticStiefel{$(2n), $(2k), $(𝔽)}()")
end

@doc raw"""
    symplectic_inverse_times(::SymplecticStiefel, p, q)
    symplectic_inverse_times!(::SymplecticStiefel, A, p, q)

Directly compute the symplectic inverse of ``p \in \operatorname{SpSt}(2n, 2k)``,
multiplied with ``q \in \operatorname{SpSt}(2n, 2k)``.
That is, this function efficiently computes
``p^+q = (Q_{2k}p^TQ_{2n})q \in ℝ^{2k \times 2k}``,
where ``Q_{2n}, Q_{2k}`` are the [`SymplecticMatrix`](@ref)
of sizes ``2n \times 2n`` and ``2k \times 2k`` respectively.

This function performs this common operation without allocating more than
a ``2k \times 2k`` matrix to store the result in, or in the case of the in-place
function, without allocating memory at all.
"""
function symplectic_inverse_times(M::SymplecticStiefel{n,k}, p, q) where {n,k}
    A = similar(p, (2k, 2k))
    return symplectic_inverse_times!(M, A, p, q)
end

function symplectic_inverse_times!(::SymplecticStiefel{n,k}, A, p, q) where {n,k}
    # we write p = [p1 p2; p3 p4] (and q, too), then
    p1 = @view(p[1:n, 1:k])
    p2 = @view(p[1:n, (k + 1):(2k)])
    p3 = @view(p[(n + 1):(2n), 1:k])
    p4 = @view(p[(n + 1):(2n), (k + 1):(2k)])
    q1 = @view(q[1:n, 1:k])
    q2 = @view(q[1:n, (k + 1):(2k)])
    q3 = @view(q[(n + 1):(2n), 1:k])
    q4 = @view(q[(n + 1):(2n), (k + 1):(2k)])
    A1 = @view(A[1:k, 1:k])
    A2 = @view(A[1:k, (k + 1):(2k)])
    A3 = @view(A[(k + 1):(2k), 1:k])
    A4 = @view(A[(k + 1):(2k), (k + 1):(2k)])
    mul!(A1, p4', q1)           # A1  = p4'q1
    mul!(A1, p2', q3, -1, 1)    # A1 -= p2'p3
    mul!(A2, p4', q2)           # A2  = p4'q2
    mul!(A2, p2', q4, -1, 1)    # A2 -= p2'q4
    mul!(A3, p1', q3)           # A3  = p1'q3
    mul!(A3, p3', q1, -1, 1)    # A3 -= p3'q1
    mul!(A4, p1', q4)           # A4  = p1'q4
    mul!(A4, p3', q2, -1, 1)    #A4  -= p3'q2
    return A
end


text:
# Integrates the following set of equations outwards:
# dy1/dx =                y2
# dy2/dx = C1 * y1 + C2 * y2
function sch_rk4_step!(
    rmesh::Array{Float64,1}, # Grid
    y0::Array{Float64,1}, # Initial condition
    C1::Array{Float64,1},
    C2::Array{Float64,1},
    C1mid::Array{Float64,1},
    C2mid::Array{Float64,1},
    y1::Array{Float64,1},
    y2::Array{Float64,1}, # Solution y1 and y2
    maxval_stop::Float64,  # Maximum value (if y1 > max_val, the integration stops)
)

    Neqn = size(y0,1)  # should be hardcoded at 2 ?
    @assert Neqn == 2

    dym = zeros(Float64,2)
    dyt = zeros(Float64,2)
    yt = zeros(Float64,2)
    dydx = zeros(Float64,2)
    y = zeros(Float64, Neqn)

    Nr = size(rmesh,1)
    y = copy(y0)
    y1[1] = y[1]
    y2[1] = y[2]

    for i in 2:Nr
        # rk4 step size
        h = rmesh[i] - rmesh[i-1]

        # k1, evaluate F
        dydx[1] = y[2]
        dydx[2] = C1[i-1] * y[1] + C2[i-1] * y[2]
        #
        yt = y + h/2 * dydx  # yi + h/2*k1
        # k2, evaluate F at x_i + h/2 and yi + k1/2
        dyt[1] = yt[2]
        dyt[2] = C1mid[i-1] * yt[1] + C2mid[i-1] * yt[2]
        #
        yt = y + h/2 * dyt # yi + h/2*k2
        #
        dym[1] = yt[2]
        dym[2] = C1mid[i-1] * yt[1] + C2mid[i-1] * yt[2]
        #
        yt = y + h * dym
        dym = dyt + dym
        dyt[1] = yt[2]
        dyt[2] = C1[i] * yt[1] + C2[i] * yt[2]
        #
        y = y + h/6 * (dydx + dyt + 2*dym)
        #
        y1[i] = y[1]
        y2[i] = y[2]

        # The integration stops at R(imax)
        if abs(y[1]) >= maxval_stop
            imax = i
            return imax
        end
    end
    imax = Nr
    return imax
end


text:
# This type represents an AND-node (Q-node) of the belief tree, branching
# on observations. It maps the set of observations seen during simulations to
# the VNodes that the simulations reach at the next level

  # Fields:
  # obs_to_particles: a mapping from observations seen to the particles (after
  # the transition) that generated the observations. The particles for each
  # observation become the representative set for the corresponding v-node
  # at the next level.
  # depth: depth of the V-node *above* this node.
  # action: action that led to this q-node.
  # first_step_reward: The average first step reward of particles when they
  # took action 'action'.
  # history: history up to the V-node *above* this node.
  # debug: Flag controlling debugging output.

mutable struct _QNode{S,A,O,B,T} #Workaround to achieve a circular type definition
    obs_to_particles::Dict{O, Vector{DESPOTParticle{S}}}
    depth::Int64
    action::A
    first_step_reward::Float64
    history::History
    weight_sum::Float64
    obs_and_nodes::Vector{Tuple{O,T}}
    n_visits::Int64                # Needed for large problems
    bounds::B

      # default constructor
      function _QNode{S,A,O,B,T}(
                    pomdp::POMDP{S,A,O},
                    bounds::B,
                    obs_to_particles::Dict{O,Vector{DESPOTParticle{S}}},
                    depth::Int64,
                    action::A,
                    first_step_reward::Float64,
                    history::History{A,O},
                    config::DESPOTConfig) where {S,A,O,B,T}

            this = new()
            this.obs_to_particles = obs_to_particles
            this.depth = depth
            this.action = action
            this.first_step_reward = first_step_reward
            this.history = history
            this.weight_sum = 0
            this.obs_and_nodes = Tuple{O,T}[]
            this.n_visits = 0
            this.bounds = bounds

            for (obs, particles) in this.obs_to_particles
                obs_weight_sum = 0.0
                for p in particles
                    obs_weight_sum += p.weight
                end
                this.weight_sum += obs_weight_sum
                add(this.history, action, obs)
                remove_last(this.history)
                push!(this.obs_and_nodes, (obs,
                                           VNode{S,A,O,B}(pomdp,
                                                          particles,
                                                          bounds,
                                                          this.depth+1,  # TODO: check depth
                                                          obs_weight_sum,
                                                          false,
                                                          config)))
            end
            return this
        end
end

# This type ("Value Node") encapsulates a belief node (and recursively, a
# belief tree). It stores the set of particles associated with the node, an
# AND-node for each action, and some bookkeeping information.

mutable struct VNode{S,A,O,B}
  particles::Array{DESPOTParticle{S},1}
  lbound::Float64
  ubound::Float64
  depth::Int64
  default_value::Float64            # Value of the default policy (equals to lbound value
                                    # before any backups are performed)
  pruned_action::A                  # Best action at the node after pruning
  weight::Float64                   # Sum of particle weights at this belief
  best_ub_action::A                 # Action that gives the highest upper bound
  best_lb_action::A                 # Action that gives the highest lower bound
  in_tree::Bool                     # True if the node is visited by trial().
                                    # In order to determine if a node is a fringe node
                                    # of the belief tree, we need to expand it one level.
                                    # The nodes added during this expansion of a fringe
                                    # node are not considered to be within the tree unless
                                    # explicitly visited by trial(), so we use
                                    # this indicator variable.
  n_tree_nodes::Int64               # Number of nodes with in_tree == true in the subtree
                                    # rooted at this node
  q_nodes::Dict{A,_QNode{S,A,O,B,
                VNode{S,A,O,B}}}    # Dict of children q-nodes
  n_visits::Int64                   # Needed for large domains
  n_actions_allowed::Int64          # current number of action branches allowed in the node, needed for large domains
  q_star::Float64                   # best current Q-value, needed for large domains

  # default constructor
  function VNode{S,A,O,B}(
               pomdp::POMDP{S,A,O},
               particles::Vector{DESPOTParticle{S}},
               b::B,
               depth::Int64,
               weight::Float64,
               in_tree::Bool,
               config::DESPOTConfig) where {S,A,O,B}

        this = new()
        this.particles          = particles
        this.lbound,
        this.ubound             = bounds(b, pomdp, particles, config)
        this.depth              = depth
        this.default_value      = this.lbound
        this.weight             = weight
        this.in_tree            = in_tree
        this.n_tree_nodes       = in_tree ? 1:0
        this.q_nodes            = Dict{A,QNode{S,A,O,B}}()
        this.n_visits           = 0
        this.n_actions_allowed  = 0
        this.q_star             = -Inf

        # note these are left undefined for now. Are they needed?
        # this.pruned_action      = #undef
        # this.best_ub_action     = #undef
        # this.best_lb_action     = #undef

        validate_bounds(this.lbound, this.ubound, config)
        return this
  end
end

const QNode{S,A,O,B} = _QNode{S,A,O,B,VNode{S,A,O,B}}

function get_upper_bound{S,A,O,B}(qnode::QNode{S,A,O,B})
  ubound::Float64 = 0.0
  for (obs, node) in qnode.obs_and_nodes
      ubound += node.ubound * node.weight
  end
  return ubound/qnode.weight_sum
end

function get_lower_bound{S,A,O,B}(qnode::QNode{S,A,O,B})
  lbound::Float64 = 0.0
  for (obs, node) in qnode.obs_and_nodes
      lbound += node.lbound * node.weight
  end
  return lbound/qnode.weight_sum
end

#TODO: Fix this

# function prune{S,A,O,L,U}(qnode::QNode{S,A,O,L,U}, total_pruned::Int64, config::DESPOTConfig)
#   cost::Float64 = 0.0
#   total_pruned::Int64 = 0
#
#   for (obs,node) in qnode.obs_to_node
#     cost, total_pruned += prune(node, total_pruned, config)
#   end
#   return cost, total_pruned
# end

function get_lb_action{S,A,O,B}(node::VNode{S,A,O,B}, config::DESPOTConfig, discount::Float64)
  local a_star
  q_star::Float64 = -Inf
  remaining_reward::Float64 = 0.0

  for (a,q_node) in node.q_nodes
    remaining_reward = get_lower_bound(q_node)
    if q_node.first_step_reward + discount * remaining_reward > q_star + config.tiny
      q_star = q_node.first_step_reward + discount * remaining_reward
      a_star = a
    end
  end
  return a_star
end

#TODO: fix pruning

# function prune{S,A,O,L,U}(node::VNode{S,A,O,L,U}, total_pruned::Int64, config::DESPOTConfig)
#   # Cost if the node was pruned
#   cost = (config.discount^node.depth) * node.weight * node.default_value
#                 - config.pruning_constant
#   if !node.inTree # Leaf
#     @assert(nodes.n_tree_nodes == 0)
#     return cost, total_pruned
#   end
#
#   for q_node in node.q_nodes
#     first_step_reward = config.discount^depth * weight * q_node.first_step_reward
#     best_child_value, total_pruned = prune(q_node, total_pruned, config)
#
#     new_cost = first_step_reward + best_child_value - config.pruning_constant
#     if new_cost > cost
#       cost = new_cost
#       pruned_action = q_node.action
#     end
#   end
#   if pruned_action == -1
#     total_pruned +=1
#   end
#   return cost, total_pruned
# end


text:
consts=@equs(c=3e8, G=6.67e-11, massSun=2e30)
E=@equ m*v^2/2=G*M*m/R
R=E*:R/:m/:v^2*2
R_BH=R&@equs(v=c, M=massSun)&consts
@test floor(R_BH.rhs)==2964
R=U(6.67e-11,:m^3/:s^2/:kg)*U(2e30,:kg)/U(3e8,:m/:s)/U(3e8,:m/:s)*2
@test R.units==:m

include("../examples/LQG.jl")


text:
"""
al_methods.jl
"""

function solve!(solver::AugmentedLagrangianSolver{T,S}) where {T,S}
    # pull vars
    convals = solver.solver_uncon.obj.convals
    solver_uncon = solver.solver_uncon
    c_max = Inf
    J = 0.
    # intialize
    set_verbosity!(solver)
    clear_cache!(solver)
    reset!(solver)
    # run AL iterations
    for i = 1:solver.opts.al_max_iterations
        # solve the unconstrained problem
        J = solve!(solver.solver_uncon)
        # record the updated information
        record_iteration!(solver, J, max_violation_penalty(convals)...)
        # exit if converged
        exit = evaluate_convergence(solver, i)
        exit && break
        # outer loop update
        update_dual_penalty!(convals)
        # reset verbosity level after it's modified
        set_verbosity!(solver)
        reset!(solver_uncon)
    end
    terminate!(solver)
    return J
end

function evaluate_convergence(solver::AugmentedLagrangianSolver, al_iterations::Int)
    # exit if the maximum number of iterations has elapsed
    if al_iterations == solver.opts.al_max_iterations
        solver.stats.status = MAX_ITERATIONS_OUTER
        return true
    end

    # exit if the maximum total number of iterations has been reached
    if solver.stats.iterations >= solver.opts.iterations
        return true
    end

    # exit if violation tolerance satisfied or the maximum penalty has been reached
    if (solver.stats.c_max[solver.stats.iterations] < solver.opts.al_vtol ||
	    solver.stats.penalty_max[solver.stats.iterations] >= solver.opts.penalty_max)
        return true
    end

    # exit if iLQR converged
    if solver.stats.status > SOLVE_SUCCEEDED
        return true
    end

    return false
end

function record_iteration!(solver::AugmentedLagrangianSolver{T,S},
                           J::T, c_max::T, μ_max::T) where {T,S}
    # Just update constraint violation and max penalty
    record_iteration!(solver.stats, c_max=c_max, penalty_max=μ_max, is_outer=true)
    vmax, vmax_info = max_violation_info(solver.solver_uncon.obj.convals)
    j = solver.stats.iterations_outer::Int
    @logmsg OuterLoop :iter value=j
    @logmsg OuterLoop :total value=solver.stats.iterations
    @logmsg OuterLoop :cost value=J
    @logmsg OuterLoop :c_max value=c_max
    if is_verbose(solver)
	    print_level(OuterLoop, global_logger())
        println("vmax: $(vmax_info)")
    end
    return nothing
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 3786dc10-46bd-11eb-2bcb-1d2ef41250aa
md"## _Variable Scope_
In this lesson we will learn what is the scope of a variable and how scopes can be used to rule when a variable should be accessible in our program.

The scope of a variable is the region of a program where the variable is known and accessible. A variable may live in two kind of scopes: the global scope or a local scope.
"

# ╔═╡ 545598f2-46bd-11eb-22f2-d959fbadc549
md" ## Scope
A variable in the global scope is accessible everywhere in the program and can be modified by any part of the code. When we define a variable in the REPL or outside of a function, for example, we create a global variable.

A variable in a local scope is only accessible in that scope and in other scopes eventually defined inside it.

In Julia there are several constructs which introduces a scope:

[table]

As you can see, some constructs introduce a global scope (for example each module has its separate global scope) and others introduce a local scope (for example functions, for loops and let blocks).

Let’s now see in more details how to work with scopes and in which scope it is better to define a variable, depending on its usage.

"

# ╔═╡ 133c2e68-46be-11eb-2445-fdfd0bdaa763
md" ## Local Scope
Let’s start with a construct we are already familiar with: a function. A function declaration introduce a scope, which means that each variable declared inside a function lives inside the function: it can be accessed freely inside the function but cannot be accessed outside the function, which is good! Thanks to this property we can use the names most suitable for our variables (x, y, z, etc.) without the risk of clashing with the declaration of other functions.

If a value computed inside a function is needed outside the function it is a good idea to return that value instead of modifying a global constant external to the function. As a rule of thumb, a function should rely only on its input parameters and return only the variable which may be useful for further computations.

Let’s see an example of a variable which exist inside a function (local scope) but doesn’t exist in the global scope:
"

# ╔═╡ 24f7ee12-46be-11eb-3c97-05eb5e3a94c5
function example1()
    z = 42
    return
end

# ╔═╡ 2e56a278-46be-11eb-2592-6bc4a396d786
z

# ╔═╡ 3131345e-46be-11eb-1f2f-5d46726895dd
md" Although it is not advisable, it is possible to specify that a variable should be accessed in the global scope through:"

# ╔═╡ 3b64e420-46be-11eb-3a85-91a45197b5f7
function example2()
    global f = 42
    return
end

# ╔═╡ 4176df1e-46be-11eb-333e-9fbdf98b9ce9
example2()

# ╔═╡ 48822334-46be-11eb-0584-ed06686b4512
f

# ╔═╡ 56d6a844-46be-11eb-24a2-c912c5f6fae9
md"A better approach is instead to return z and let the user perform the allocation of z:"

# ╔═╡ 5e48ad82-46be-11eb-2af9-65d9f9dc3706
function example3()
    a = 42
    return a
end

# ╔═╡ 6628c890-46be-11eb-1b5f-a9dbe523ddb0
a = example3()

# ╔═╡ 6cea7e2c-46be-11eb-2a3e-b57d731a344c
a

# ╔═╡ 71dc1a96-46be-11eb-3b0e-e9ebcdb58206
md"In the case where it is necessary to distinguish between a variable which exists both in the local and global scope, it is possible to indicate the one that needs to be used through the global or local annotation before the desired variable."

# ╔═╡ 842e513c-46be-11eb-2754-fb974f8215e1
md" ## `Let` construct
The let construct can be used to introduce a new local scope. It is useful, for example, when you want to perform some computations but you don’t want the intermediate results/variables to pollute your current scope.

The let block will be able to access all the local (or global) variables available in its parent scope and will have its own set of local variables. It is also possible to specify some initial values to mimic the execution of a function once.
"

# ╔═╡ d303a000-46be-11eb-176d-65500635036f
b = let
    i=3
    i+=5
    i # the value returned from the computation
end

# ╔═╡ ddf0800a-46be-11eb-0caf-4fe37ba18934
b

# ╔═╡ e2dac2ce-46be-11eb-2e88-497a6c781222
d = let i=5
    i+=42
    i
end

# ╔═╡ ead5ea30-46be-11eb-2f46-658d41f10b20
d

# ╔═╡ f0dac34c-46be-11eb-14e8-e5ecb26ad036
e = let i=10
    i+=42
    i
end

# ╔═╡ f84bf7fe-46be-11eb-0416-79962d5dc4c6
e

# ╔═╡ fc3b5f26-46be-11eb-2655-6518d17e8ed4
md" As you can see, let blocks are pretty convenient when it comes to splitting computations over several lines, but there are other possible uses as shown here.

let blocks are somewhat similar to begin blocks, although begin blocks don’t introduce a local scope:"

# ╔═╡ 07e31486-46bf-11eb-11a0-ebb16fec9196
h = begin
    i=41
    i+=1
    i
end

# ╔═╡ 11e32cdc-46bf-11eb-105d-797889dcae00
h

# ╔═╡ 1ab5cdc4-46bf-11eb-0012-2b62db002c87
i

# ╔═╡ 1cc2fa2e-46bf-11eb-353d-bdd66ee27a35
md" ## Global Scope
Whenever we define a variable in the REPL or in general outside a construct which introduces a local scope, we place a variable in the global scope. The global scope is accessible everywhere in the program and a variable in the global scope can be modified by any part of the code. As such, it is generally advisable to avoid using global variables as much as possible, in fact since global variables can change their type and value at any time, they cannot be properly optimised by the compiler.

## Constants
One way to mitigate this performance issue is to define global variables as constants through the const annotation. When we think of a constant we generally imagine a variable which is immutable, for example the speed of light: the speed of light is a number and so it could be expressed by a Float64, there is no need to change its type, once it has been defined, to a String, in this case.

A constant in Julia is a variable which cannot change its type once it is defined (Julia will throw an error if there is an attempt to modify the type of a constant) and Julia will show a warning message if we try to modify the value of a constant. Since a constant is “type immutable” it can be properly optimised by the compiler and there are fewer performance issues.
"

# ╔═╡ 388c5d7c-46bf-11eb-1895-838457ae5907
const C = 299792458 # m / s, this is an Int

# ╔═╡ 442e44ba-46bf-11eb-21c3-171917b0c529
C1 = 300000000 # change the value of C

# ╔═╡ 5e76217e-46bf-11eb-1628-ef847aaa8267
C3 = 2.998 * 1e8 # change the type of C, not permitted

# ╔═╡ 7283fbde-46bf-11eb-2a2f-09036a699556
md" ## Modules
As we have anticipated, modules introduce separate global scope, which means that a variable which is known inside a module will not be accessible outside of the module unless it is exported or it is accessed through the ModuleName.varName notation.
"

# ╔═╡ 9970bdb8-46bf-11eb-3e48-0387b0d3f4f5
module ScopeTestModule
export a1
a1 = 25
b1 = 42
end

# ╔═╡ 9ca99a52-46bf-11eb-04b0-73214c9cc4cc
a1

# ╔═╡ a5560d4a-46bf-11eb-1c4b-6993cf3cbb7b
b1

# ╔═╡ c2cfd8ba-46bf-11eb-2072-e7b1f8a7f013
ScopeTestModule.b1

# ╔═╡ c6d83970-46bf-11eb-2ec1-cf03f0290a6b
ScopeTestModule.b1=26

# ╔═╡ d0caf5ee-46bf-11eb-1120-65d881285044
md"At line 9-10 we can see that the a1 variable, which is exported by the module, can be accessed directly without specifying the scope of the variable, while on line 12-13 and 15-16 we can see that b1 can only be accessed by specifying where the variable lives (i.e. inside ScopeTestModule). At line 18-19 we can see that it is not possible to directly modify a variable which is defined inside another module.
"

# ╔═╡ Cell order:
# ╟─3786dc10-46bd-11eb-2bcb-1d2ef41250aa
# ╟─545598f2-46bd-11eb-22f2-d959fbadc549
# ╟─133c2e68-46be-11eb-2445-fdfd0bdaa763
# ╠═24f7ee12-46be-11eb-3c97-05eb5e3a94c5
# ╠═2e56a278-46be-11eb-2592-6bc4a396d786
# ╟─3131345e-46be-11eb-1f2f-5d46726895dd
# ╠═3b64e420-46be-11eb-3a85-91a45197b5f7
# ╠═4176df1e-46be-11eb-333e-9fbdf98b9ce9
# ╠═48822334-46be-11eb-0584-ed06686b4512
# ╟─56d6a844-46be-11eb-24a2-c912c5f6fae9
# ╠═5e48ad82-46be-11eb-2af9-65d9f9dc3706
# ╠═6628c890-46be-11eb-1b5f-a9dbe523ddb0
# ╠═6cea7e2c-46be-11eb-2a3e-b57d731a344c
# ╟─71dc1a96-46be-11eb-3b0e-e9ebcdb58206
# ╟─842e513c-46be-11eb-2754-fb974f8215e1
# ╠═d303a000-46be-11eb-176d-65500635036f
# ╠═ddf0800a-46be-11eb-0caf-4fe37ba18934
# ╠═e2dac2ce-46be-11eb-2e88-497a6c781222
# ╠═ead5ea30-46be-11eb-2f46-658d41f10b20
# ╠═f0dac34c-46be-11eb-14e8-e5ecb26ad036
# ╠═f84bf7fe-46be-11eb-0416-79962d5dc4c6
# ╟─fc3b5f26-46be-11eb-2655-6518d17e8ed4
# ╠═07e31486-46bf-11eb-11a0-ebb16fec9196
# ╠═11e32cdc-46bf-11eb-105d-797889dcae00
# ╠═1ab5cdc4-46bf-11eb-0012-2b62db002c87
# ╟─1cc2fa2e-46bf-11eb-353d-bdd66ee27a35
# ╠═388c5d7c-46bf-11eb-1895-838457ae5907
# ╠═442e44ba-46bf-11eb-21c3-171917b0c529
# ╠═5e76217e-46bf-11eb-1628-ef847aaa8267
# ╠═7283fbde-46bf-11eb-2a2f-09036a699556
# ╠═9970bdb8-46bf-11eb-3e48-0387b0d3f4f5
# ╠═9ca99a52-46bf-11eb-04b0-73214c9cc4cc
# ╠═a5560d4a-46bf-11eb-1c4b-6993cf3cbb7b
# ╠═c2cfd8ba-46bf-11eb-2072-e7b1f8a7f013
# ╠═c6d83970-46bf-11eb-2ec1-cf03f0290a6b
# ╟─d0caf5ee-46bf-11eb-1120-65d881285044


text:
using DemultiplexFASTQ
using Test

@testset "DemultiplexFASTQ.jl" begin
    # Write your tests here.
end


text:
HR8799_desc = """
# HR8799 NIRC2/Vortex ADI Data

Website: https://doi.org/10.5281/zenodo.4142348

Author: Michael Bottom et al.

Date of Publication: October 27, 2020

License: https://creativecommons.org/licenses/by/4.0/legalcode

Processed data for HR8799 captured using Keck II's NIRC2 with the Vortex Coronagraph. The data was processed in preparation for a journal article and is further described on this page.

Please cite this dataset: Bottom, M., Mawet, D., and Ruane, G. (2020). *HR8799 NIRC2/Vortex Processed Data* [Data set]. Zenodo. https://doi.org/10.5281/ZENODO.4142348

## Keys
* `:cube` pre-processed ADI cube
* `:pa` parallactic angles
"""

HR8799_datadep = DataDep(
    "HR8799",
    HR8799_desc,
    [
        "https://zenodo.org/record/4142349/files/HR8799_cube.fits?download=1",
        "https://zenodo.org/record/4142349/files/HR8799_pa.fits?download=1"
    ],
    "950102a5133eaf08c9e7b15033f570df2eea968ab3db5e20e3ea6243708fb193"
)

struct HR8799Dataset <: HCIDataset end

@doc HR8799_desc
const HR8799 = HR8799Dataset()

Base.keys(::HR8799Dataset) = (:cube, :pa)

function Base.getindex(::HR8799Dataset, key::Symbol)
    if key === :cube
        getdata(datadep"HR8799/HR8799_cube.fits")
    elseif key === :pa
        getdata(datadep"HR8799/HR8799_pa.fits")
    else
        throw(KeyError(key))
    end
end


text:
using TcpInstruments
using Test

s = initialize(AgilentDSOX4034A)

timeit(func_time, title) = println("$title: $(func_time)")

for i in 1:2
println("Single Channe: Iteration $i")

@show @elapsed get_data(s, 1)

@show @elapsed get_data(s, [1])

@show @elapsed get_data(s, [1]; inbounds=true)

println("Two Channel Test: Iteration $i")

@show @elapsed get_data(s, [1, 2])

@show @elapsed get_data(s)

@show @elapsed get_data(s, [1, 2]; inbounds=true)

@show @elapsed get_data(s; inbounds=true)
end


text:
try
	using RCall
	RCall.eval("library(locfit);library(ggplot2)")
catch
	using Pkg
	ENV["R_HOME"]="*"

	Pkg.add("Conda")

	using Conda
	Conda.add("r-base",channel="conda-forge")
	Conda.add(["r-locfit","r-ggplot2","r-data.table","r-r.utils"],channel="conda-forge")

	Pkg.add("RCall")

	using RCall
end

"""
	Estimating and plotting MAP using locfit and ggplot2 in R. It assume your folder contains the posterior estimated through ABCreg
"""
function plot_map(;analysis_folder::String,weak::Bool=true,title::String="Posteriors")

	R"""library(locfit);library(ggplot2);library(data.table)""";

	out = filter(x -> occursin("post",x), readdir(analysis_folder,join=true))
	out          = filter(x -> !occursin(".1.",x),out)

	open(x)      = Array(CSV.read(x,DataFrame))

	# Control outlier inference. 2Nes non negative values
	flt(x)       = x[(x[:,4] .> 0) .& (x[:,1] .> 0) .& (x[:,2] .> 0) .& (x[:,3] .> 0) ,:]
	posteriors   = flt.(open.(out))

	R"""getmap <- function(df){
			temp = as.data.frame(df)
			d <-locfit(~temp[,1],temp);
			map<-temp[,1][which.max(predict(d,newdata=temp))]
		}"""

	getmap(x)    = rcopy(R"""suppressWarnings(matrix(apply($x,2,getmap),nrow=1))""")

	if !weak
		posteriors = [posteriors[i][:,3:end] for i in eachindex(posteriors)]
		tmp          = getmap.(posteriors)
		maxp         = DataFrame(vcat(tmp...),[:a,:gam_neg,:shape])
		al           = maxp[:,1:1]
		gam          = maxp[:,2:end]
	else
		tmp          = getmap.(posteriors)
		maxp         = DataFrame(vcat(tmp...),[:aw,:as,:a,:gam_neg,:shape])
		al           = maxp[:,1:3]
		gam          = maxp[:,4:end]
	end

	@rput al
	@rput title
	@rput posteriors
	@rput analysis_folder

	R"""al = as.data.table(al)
			lbls = if(ncol(al) > 1){c(expression(paste('Posterior ',alpha[w])), expression(paste('Posterior ',alpha[s])),expression(paste('Posterior ',alpha)))}else{c(expression(paste('Posterior ',alpha)))}
			clrs = if(ncol(al) > 1){c('#30504f', '#e2bd9a', '#ab2710')}else{c('#ab2710')}

			if(nrow(al) == 1){
				tmp = as.data.table(posteriors[1])
				dal = suppressWarnings(data.table::melt(tmp[,1:3]))
				pal = ggplot(dal) + geom_density(aes(x=value,fill=variable),alpha=0.75) + scale_fill_manual('Posterior distribution',values=clrs ,labels=lbls) + theme_bw() + ggtitle(title) + xlab(expression(alpha)) + ylab("")
			}else{
				dal = suppressWarnings(data.table::melt(al))
				pal = suppressWarnings(ggplot(dal) + geom_density(aes(x=value,fill=variable),alpha=0.75) + scale_fill_manual('Posterior distribution',values=clrs ,labels=lbls) + theme_bw() + ggtitle(title) + xlab(expression(alpha)) + ylab(""))
			}
			suppressWarnings(ggsave(pal,filename=paste0(analysis_folder,'/map.png'),dpi=600))
			"""
	CSV.write(analysis_folder * "/map.tsv",maxp,delim='\t',header=true)

	return(posteriors,maxp)
end


text:
using ForwardDiff, Roots, Plots, LinearAlgebra, Optim

f(x) = ((x[1]-1.0)^2) + (0.1*(x[2]-4)^2) + sin(5*pi*x[1])+sin(5*pi*x[2])
grad_f(x) = ForwardDiff.gradient(f,x)

# Ponto inicial
xk = [-1.0,1.0]

#direcao de busca
dk = -grad_f(xk)/norm(grad_f(xk))

beta = 0.1
alpha_f = 1.0

phi(alpha) = f(xk+alpha*dk)
d_phi(alpha) = grad_f(xk+alpha*dk)'dk

l_alpha(alpha) = phi(0)+alpha*beta*d_phi(0)
sigma = 0.2

wolf = []

for alpha = 0:0.01:alpha_f
    if (phi(alpha) <= l_alpha(alpha)) && (grad_f(xk+alpha*dk)'dk >= sigma*grad_f(xk)'dk)
        push!(wolf,alpha)
    end
end

plot(phi,0,alpha_f,lw=2,lab="phi")
plot!(l_alpha,0,alpha_f,lw=2,lab="L_phi",c=:black)
scatter!(wolf,phi.(wolf),lab="Wolfe",ms=3,c=:red)
plot!(legend=:outerbottomright)
plot!(title="beta = 0.1 e sigma = 0.2",xlabel = "alpha", ylabel = "phi(alpha), l_alpha(alpha")


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

using Markdown
using InteractiveUtils

# ╔═╡ 2ade0260-68da-11eb-3c47-51af562ca746
begin
	using PlutoUI
	using NativeSVG
end

# ╔═╡ 02c52b84-6991-11eb-1dbd-93af9a7ca7f1
let
	include("../src/chap04.jl")
	🐢 = Turtle()
	penup(🐢)
	forward(🐢, -200)
	pendown(🐢)
	flower(🐢, 7, 110.0, 60.0)
	penup(🐢)
	forward(🐢, 200)
	pendown(🐢)
	flower(🐢, 10, 80.0, 80.0)
	penup(🐢)
	forward(🐢, 200)
	pendown(🐢)
	flower(🐢, 20, 280.0, 20.0)
	Drawing(🐢, 720, 220)
end

# ╔═╡ 7933043e-6992-11eb-27b6-c10f5cb6dd29
let
	include("../src/chap04.jl")
	🐢 = Turtle()
	penup(🐢)
	forward(🐢, -280)
	pendown(🐢)
	size = 80
	drawpie(🐢, 5, size)
	drawpie(🐢, 6, size)
	drawpie(🐢, 7, size)
	drawpie(🐢, 8, size)
	Drawing(🐢, 720, 160)
end

# ╔═╡ 31bc6732-6993-11eb-295e-eb9f409a95fa
let
	include("../src/chap04.jl")
	🐢 = Turtle()
	penup(🐢)
	turn(🐢, 90)
	forward(🐢, 25)
	turn(🐢, -90)
	pendown(🐢)
  	spiral(🐢, 230, 6, 0.1, 0.0002)
	Drawing(🐢, 720, 220)
end

# ╔═╡ dc4299ea-68d9-11eb-273b-a597294ef196
md"""# Case Study: Interface Design

This chapter presents a case study that demonstrates a process for designing functions that work together.

It introduces turtle graphics, a way to create programmatic drawings. Turtle graphics are not included in the standard library, so to use them you’ll have to add the `NativeSVG` package to your Julia setup.

The examples in this chapter can be executed in a graphical web-based notebook using the `Pluto` package, which combines code, formatted text, math, and multimedia in a single document (see Appendix B).
"""

# ╔═╡ 3dcb939c-68da-11eb-2546-85b34837f29c
md"""## Turtles

A *module* contains a collection of related functions. Julia provides some modules in its standard library. Additional functionality can be added from a growing collection of *packages*.

Packages can be installed in the REPL:

```julia
julia> using Pkg

julia> pkg.add("https://github.com/BenLauwens/NativeSVG.jl.git")
```

This can take some time.

Before we can use the functions in a module, we have to import it with a using statement:

```julia
using NativeSVG
🐢 = Turtle()
```

The NativeSVG module provides a function called `Turtle` that creates a turtle object, which we assign to a variable named `🐢` (**`\:turtle: TAB`**).

Once you create a turtle, you can call a function to move it around. For example, to move the turtle forward:

```julia
forward(🐢, 100)
Drawing(🐢, 720, 10)
```

The `Drawing` function creates an SVG picture (Figure 4-1). Its first argument is the turtle object, the second the width of the drawing and the third the height.
"""

# ╔═╡ 319c52f2-68dd-11eb-130c-f3672af50055
let
	🐢 = Turtle()
	forward(🐢, 100)
	Drawing(🐢, 720, 10)
end

# ╔═╡ b1f3a4cc-68de-11eb-3de9-87b14a76490e
md"*Figure 4-1. Moving the turtle forward.*"

# ╔═╡ 09d4cbb2-68df-11eb-2eec-bbe39ef9880d
md"""The arguments of `forward` are the turtle and a distance in pixels, so the actual size of the line that’s drawn depends on your display.

!!! info
    Each turtle is holding a pen, which is either down or up; if the pen is down (the default), the turtle leaves a trail when it moves. Figure 4-1 shows the trail left behind by the turtle. To move the turtle without drawing a line, first call the function `penup`. To start drawing again, call `pendown`.

Another function you can call with a turtle as an argument is `turn` for turning. The second argument for turn is an angle in degrees.

To draw a right angle:

```julia
let
	🐢 = Turtle()
	forward(🐢, 100)
	turn(🐢, -90)
	forward(🐢, 100)
	Drawing(🐢, 720, 210)
end
```

The keyword `let` is used to group and isolate the statements. The result of the last statement in the let block is visualized in a Pluto notebook.
"""

# ╔═╡ c64b81dc-68df-11eb-1086-d57b61651791
md"""#### Exercise 4-1

Now modify the macro to draw a square. Don’t go on until you’ve got it working!
"""

# ╔═╡ 717ad810-6975-11eb-262d-7bc4027d58b6
md"""## Simple Repetition

Chances are you wrote something like this:

```julia
let
	🐢 = Turtle()
	forward(🐢, 100)
	turn(🐢, -90)
	forward(🐢, 100)
	turn(🐢, -90)
	forward(🐢, 100)
	turn(🐢, -90)
	forward(🐢, 100)
	Drawing(🐢, 720, 210)
end
```

We can do the same thing more concisely with a `for` statement:

```julia
julia> for i in 1:4
           println("Hello!")
       end
Hello!
Hello!
Hello!
Hello!
```

This is the simplest use of the `for` statement; we will see more later. But that should be enough to let you rewrite your square-drawing program. Don’t go on until you do.

Here is a `for` statement that draws a square:

```julia
let
	🐢 = Turtle()
	for i in 1:4
		forward(🐢, 100)
		turn(🐢, -90)
	end
	Drawing(🐢, 720, 210)
end
```

The syntax of a `for` statement is similar to a function definition. It has a header and a body that ends with the keyword `end`. The body can contain any number of state‐ ments.

A `for` statement is also called a *loop* because the flow of execution runs through the body and then loops back to the top. In this case, it runs the body four times.

This version is actually a little different from the previous square-drawing code because it makes another turn after drawing the last side of the square. The extra turn takes more time, but it simplifies the code if we do the same thing every time through the loop. This version also has the effect of leaving the turtle back in the starting position, facing in the starting direction.
"""

# ╔═╡ 81fd8b28-6976-11eb-2737-69f3b0b4624c
md"""
!!! languages
    In Python the `for` statement is very similar:

    ```python
    >>> for i in range(4):
    ...     print('Hello!')
    ...
    ```

    As in Julia the `for` statement ressembles the function definition with a header that ends with a colon and an indented body. An empty line closes the loop.

    The `for` statement is the C programming language is less trivial:

    ```c
    for (int i=0; i<4; i++) {
        printf("Hello!\n");
    }
    ```

    The header has three arguments: initial value, end condition and increment. The body is enclosed in brackets.
"""

# ╔═╡ dd37018e-6978-11eb-2885-29267aaac26a
md"""## Exercises

The following is a series of exercises using turtles. They are meant to be fun, but they have a point, too. While you are working on them, think about what the point is.

!!! tip
    The following sections contain solutions to the exercises, so don’t look until you have finished (or at least tried them).

#### Exercise 4-2

Write a function called `square` that takes a parameter named `t`, which is a turtle. It should use the turtle to draw a square.

Write a function call that passes `t` as an argument to `square`, and then call `Drawing`.

#### Exercise 4-3

Add another parameter, named `len`, to `square`. Modify the body so the length of the sides is `len`, and then modify the function call to provide a second argument. Call `Drawing again`. Test with a range of values for `len`.

#### Exercise 4-4

Make a copy of `square` and change the name to `polygon`. Add another parameter named `n` and modify the body so it draws an ``n``-sided regular polygon.

!!! tip
    The exterior angles of an ``n``-sided regular polygon are ``\frac{360}{n}`` degrees.

#### Exercise 4-5

Write a function called `circle` that takes a turtle, `t`, and radius, `r`, as parameters and that draws an approximate circle by calling polygon with an appropriate length and number of sides. Test your function with a range of values of `r`.

!!! tip
    Figure out the circumference of the circle and make sure that `len * n == circumference`.

#### Exercise 4-6

Make a more general version of circle called `arc` that takes an additional parameter `angle`, which determines what fraction of a circle to draw. `angle` is in units of degrees, so when `angle = 360`, `arc` should draw a complete circle.
"""

# ╔═╡ 25a280d6-697b-11eb-3279-95bb0e4d8229
md"""## Encapsulation

The first exercise asks you to put your square-drawing code into a function definition and then call the function, passing the turtle as a parameter. Here is a solution:

```julia
function square(t)
	for i in 1:4
		forward(t, 100)
		turn(t, -90)
	end
end
```

To draw the square, a call to `Drawing` is needed:

```julia
let
	🐢 = Turtle()
	square(🐢)
	Drawing(🐢, 720, 210)
end
```

The innermost statements, `forward` and `turn`, are indented twice to show that they are inside the for loop, which is inside the function definition.

Inside the function, `t` refers to the same turtle `🐢`, so `turn(t, -90)` has the same effect as `turn(🐢, -90)`. In that case, why not call the parameter `🐢`? The idea is that `t` can be any turtle, not just `🐢` so you could create a second turtle and pass it as an argument to `square`:

```julia
let
	🐫 = Turtle()
	square(🐫)
	Drawing(🐫, 720, 210)
end
```

Wrapping a piece of code up in a function is called *encapsulation*. One of the benefits of encapsulation is that it attaches a name to the code, which serves as a kind of documentation. Another advantage is that if you reuse the code, it is more concise to call a function twice than to copy and paste the body!
"""

# ╔═╡ c9fa13ec-697b-11eb-1a40-4b899f07dde7
md"""## Generalization

The next step is to add a `len` parameter to `square`. Here is a solution:

```julia
function square(t, len)
	for i in 1:4
    	forward(t, len)
	turn(t, -90)
end
```

Adding a parameter to a function is called *generalization* because it makes the function more general. In the previous version, the square is always the same size; in this version it can be any size.

The next step is also a generalization. Instead of drawing squares, polygon draws regular polygons with any number of sides. Here is a solution:

```julia
function polygon(t, n, len)
	angle = 360 / n
	for i in 1:n
		forward(t, len)
		turn(t, -angle)
	end
end
```

To draw a 7-sided polygon with side length 70:

```julia
let
	🐢 = Turtle()
	polygon(🐢, 7, 70)
	Drawing(🐢, 720, 320)
end
```
"""

# ╔═╡ 372310ae-697f-11eb-3217-f1afed271f4a
md"""## Interface Design

The next step is to write `circle`, which takes a radius, `r`, as a parameter. Here is a simple solution that uses `polygon` to draw a 50-sided polygon:

```julia
function circle(t, r)
	circumference = 2π * r
	n=50
	len = circumference / n
	polygon(t, n, len)
end
```

The first line computes the circumference of a circle with radius ``r`` using the formula ``2πr``. `n` is the number of line segments in our approximation of a circle, so len is the length of each segment. Thus, polygon draws a 50-sided polygon that approximates a circle with radius ``r``.

One limitation of this solution is that `n` is a constant, which means that for very big circles, the line segments are too long, and for small circles, we waste time drawing very small segments. One solution would be to generalize the function by taking `n` as a parameter. This would give the user (whoever calls `circle`) more control, but the interface would be less clean.

The *interface* of a function is a summary of how it is used: What are the parameters? What does the function do? And what is the return value? An interface is “clean” if it allows the caller to do what he wants without dealing with unnecessary details.

In this example, `r` belongs in the interface because it specifies the circle to be drawn. `n` is less appropriate because it pertains to the details of how the circle should be rendered.

Rather than cluttering up the interface, it is better to choose an appropriate value of `n` depending on circumference:

```julia
function circle(t, r)
	circumference = 2π * r
	n = trunc(circumference / 3) + 3
	len = circumference / n
	polygon(t, n, len)
end
```

Now the number of segments is an integer near `circumference/3`, so the length of each segment is approximately ``3``, which is small enough that the circles look good but big enough to be efficient, and acceptable for any size circle.

Adding `3` to `n` guarantees that the polygon has at least three sides.
"""

# ╔═╡ 1b946992-6980-11eb-14ee-a11b07a08c75
md"""## Refactoring

When I wrote `circle`, I was able to reuse `polygon` because a many-sided polygon is a good approximation of a circle. But `arc` is not as cooperative; we can’t use `polygon` or `circle` to draw an arc.

One alternative is to start with a copy of `polygon` and transform it into `arc`. The result might look like this:

```julia
function arc(t, r, angle)
	arc_len = 2π * r * angle / 360
	n = trunc(arc_len / 3) + 1
	step_len = arc_len / n
	step_angle = angle / n
	for i in 1:n
    	forward(t, step_len)
		turn(t, -step_angle)
	end
end
```

The second half of this function looks like `polygon`, but we can’t reuse `polygon` without changing the interface. We could generalize `polygon` to take an `angle` as a third argument, but then `polygon` would no longer be an appropriate name! Instead, let’s call the more general function `polyline`:

```julia
function polyline(t, n, len, angle)
	for i in 1:n
    	forward(t, len)
		turn(t, -angle)
	end
end
```

Now we can rewrite `polygon` and `arc` to use `polyline`:

```julia
function polygon(t, n, len)
	angle = 360 / n
	polyline(t, n, len, angle)
end

function arc(t, r, angle)
	arc_len = 2π * r * angle / 360
	n = trunc(arc_len / 3) + 1
	step_len = arc_len / n
	step_angle = angle / n
	polyline(t, n, step_len, step_angle)
end
```

Finally, we can rewrite `circle` to use `arc`:

```julia
function circle(t, r)
	arc(t, r, 360)
end
```

This process—rearranging a program to improve interfaces and facilitate code reuse —is called *refactoring*. In this case, we noticed that there was similar code in `arc` and `polygon`, so we “factored it out” into `polyline`.

If we had planned ahead, we might have written `polyline` first and avoided refactoring, but often you don’t know enough at the beginning of a project to design all the interfaces. Once you start coding, you understand the problem better. Sometimes refactoring is a sign that you have learned something.
"""

# ╔═╡ 090f60d4-6981-11eb-1c3d-df12052e6c4f
md"""## A Development Plan

A *development plan* is a process for writing programs. The process we used in this case study is “encapsulation and generalization.” The steps of this process are:

1. Start by writing a small program with no function definitions.
2. Once you get the program working, identify a coherent piece of it, encapsulate the piece in a function, and give it a name.
3. Generalize the function by adding appropriate parameters.
4. Repeat steps 1–3 until you have a set of working functions. Copy and paste working code to avoid retyping (and redebugging).
5. Look for opportunities to improve the program by refactoring. For example, if you have similar code in several places, consider factoring it into an appropriately general function.

This process has some drawbacks—we will see alternatives later—but it can be useful if you don’t know ahead of time how to divide the program into functions. This approach lets you design as you go along.
"""

# ╔═╡ 31ff7108-6981-11eb-1fb7-2ff7384538ab
md"""## Docstring

A *docstring* is a string before a function that explains the interface (“doc” is short for “documentation”). Here is an example:

```julia
\"\"\"
polyline(t, n, len, angle)

Draws n line segments with the given length and angle (in degrees) between them. t is a turtle.
\"\"\"
function polyline(t, n, len, angle)
	for i in 1:n
		forward(t, len)
		turn(t, -angle)
	end
end
```

Documentation can be accessed in the REPL or in a notebook by typing `?` followed by the name of a function, and pressing Enter:

```julia
help?> polyline
search: polyline

  polyline(t, n, len, angle)

  Draws n line segments with the given length and angle (in degrees) between them. t is a turtle.
```

Docstrings are often triple-quoted strings, also known as “multiline” strings because the triple quotes allow the string to span more than one line.

A docstring contains the essential information someone would need to use the function. It explains concisely what the function does (without getting into the details of how it does it). It explains what effect each parameter has on the behavior of the function and what type each parameter should be (if it is not obvious).

!!! tip
    Writing this kind of documentation is an important part of interface design. A well-designed interface should be simple to explain; if you have a hard time explaining one of your functions, maybe the interface could be improved.
"""

# ╔═╡ d24c63c8-6981-11eb-3be3-2f258a4436c5
md"""## Debugging

An interface is like a contract between a function and a caller. The caller agrees to provide certain parameters and the function agrees to do certain work.

For example, `polyline` requires four arguments: `t` has to be a turtle; `n` has to be an integer; `len` should be a positive number; and `angle` has to be a number, which is understood to be in degrees.

These requirements are called *preconditions* because they are supposed to be true before the function starts executing. Conversely, conditions at the end of the function are *postconditions*. Postconditions include the intended effect of the function (like drawing line segments) and any side effects (like moving the turtle or making other changes).

Preconditions are the responsibility of the caller. If the caller violates a (properly documented!) precondition and the function doesn’t work correctly, the bug is in the caller, not the function.

If the preconditions are satisfied and the postconditions are not, the bug is in the function. If your pre- and postconditions are clear, they can help with debugging.
"""

# ╔═╡ 4f879090-698e-11eb-2aa6-37a98c300038
md"""## Glossary

*module*:
A collection of related functions and other definitions.

*package*:
An external library with additional functionality.

*using statement*:
A statement that reads a module file and creates a module object.

*loop*:
A part of a program that can run repeatedly.

*encapsulation*:
The process of transforming a sequence of statements into a function definition.

*generalization*:
The process of replacing something unnecessarily specific (like a number) with something appropriately general (like a variable or parameter).

*interface*:
A description of how to use a function, including the name and descriptions of the arguments and return value.

*refactoring*:
The process of modifying a working program to improve function interfaces and other qualities of the code.

*development plan*:
A process for writing programs.

*docstring*:
A string that appears at the top of a function definition to document the func‐ tion’s

*precondition*:
A requirement that should be satisfied by the caller before a function starts.

*postcondition*:
A requirement that should be satisfied by the function before it ends.
"""

# ╔═╡ d2baf7ea-698e-11eb-1261-db00582907a3
md"## Exercises"

# ╔═╡ e3fd2546-698e-11eb-3031-c1cde742aa22
md"""#### Exercise 4-7

Enter the code in this chapter in a Pluto notebook.

1. Draw a stack diagram that shows the state of the program while executing `circle(🐢, radius)`. You can do the arithmetic by hand or add print statements to the code.

2. The version of `arc` in “Refactoring” is not very accurate because the linear approximation of the circle is always outside the true circle. As a result, the turtle ends up a few pixels away from the correct destination. The solution shown here illustrates a way to reduce the effect of this error. Read the code and see if it makes sense to you. If you draw a diagram, you might see how it works.

   ```julia
   \"\"\"
   arc(t, r, angle)

   Draws an arc with the given radius and angle:

       t: turtle
       r: radius
       angle: angle subtended by the arc, in degrees
   \"\"\"
   function arc(t, r, angle)
       arc_len = 2π * r * abs(angle) / 360
       n = trunc(arc_len / 4) + 3
       step_len = arc_len / n
       step_angle = angle / n

       # making a slight left turn before starting reduces
       # the error caused by the linear approximation of the arc
       turn(t, -step_angle/2)
       polyline(t, n, step_len, step_angle)
       turn(t, step_angle/2)
   end
   ```
"""

# ╔═╡ 9beec042-698f-11eb-3808-058f5e8e3bab
md"""#### Exercise 4-8
Write an appropriately general set of functions that can draw flowers as in Figure 4-2.
"""

# ╔═╡ 6168b66a-6991-11eb-1721-37fcb97a3b8b
md"*Figure 4-2. Turtle flowers.*"

# ╔═╡ 68fa7d90-6992-11eb-06f8-01089d2b49b9
md"""#### Exercise 4-9

Write an appropriately general set of functions that can draw shapes as in Figure 4-3.
"""

# ╔═╡ c791f734-6992-11eb-2187-8b79c9ae0e61
md"*Figure 4-3. Turtle pies.*"

# ╔═╡ 9e8cecc2-6992-11eb-39e9-fd8b983a88f7
md"""#### Exercise 4-10

The letters of the alphabet can be constructed from a moderate number of basic elements, like vertical and horizontal lines and a few curves. Design an alphabet that can be drawn with a minimal number of basic elements and then write functions that draw the letters.

You should write one function for each letter, with names `draw_a`, `draw_b`, etc., and put your functions in a file named *letters.jl*.
"""

# ╔═╡ 0e943bec-6993-11eb-39d4-b3350f7b238c
md"""#### Exercise 4-11

Read about spirals at [https://en.wikipedia.org/wiki/Spiral](https://en.wikipedia.org/wiki/Spiral); then write a program that draws an Archimedean spiral as in Figure 4-4.
"""

# ╔═╡ 03c84024-6994-11eb-1619-6b1d19725997
md"*Figure 4-4. Archimedean spiral.*"

# ╔═╡ Cell order:
# ╟─2ade0260-68da-11eb-3c47-51af562ca746
# ╟─dc4299ea-68d9-11eb-273b-a597294ef196
# ╟─3dcb939c-68da-11eb-2546-85b34837f29c
# ╟─319c52f2-68dd-11eb-130c-f3672af50055
# ╟─b1f3a4cc-68de-11eb-3de9-87b14a76490e
# ╟─09d4cbb2-68df-11eb-2eec-bbe39ef9880d
# ╟─c64b81dc-68df-11eb-1086-d57b61651791
# ╟─717ad810-6975-11eb-262d-7bc4027d58b6
# ╟─81fd8b28-6976-11eb-2737-69f3b0b4624c
# ╟─dd37018e-6978-11eb-2885-29267aaac26a
# ╟─25a280d6-697b-11eb-3279-95bb0e4d8229
# ╟─c9fa13ec-697b-11eb-1a40-4b899f07dde7
# ╟─372310ae-697f-11eb-3217-f1afed271f4a
# ╟─1b946992-6980-11eb-14ee-a11b07a08c75
# ╟─090f60d4-6981-11eb-1c3d-df12052e6c4f
# ╟─31ff7108-6981-11eb-1fb7-2ff7384538ab
# ╟─d24c63c8-6981-11eb-3be3-2f258a4436c5
# ╟─4f879090-698e-11eb-2aa6-37a98c300038
# ╟─d2baf7ea-698e-11eb-1261-db00582907a3
# ╟─e3fd2546-698e-11eb-3031-c1cde742aa22
# ╟─9beec042-698f-11eb-3808-058f5e8e3bab
# ╠═02c52b84-6991-11eb-1dbd-93af9a7ca7f1
# ╟─6168b66a-6991-11eb-1721-37fcb97a3b8b
# ╟─68fa7d90-6992-11eb-06f8-01089d2b49b9
# ╟─7933043e-6992-11eb-27b6-c10f5cb6dd29
# ╟─c791f734-6992-11eb-2187-8b79c9ae0e61
# ╟─9e8cecc2-6992-11eb-39e9-fd8b983a88f7
# ╟─0e943bec-6993-11eb-39d4-b3350f7b238c
# ╟─31bc6732-6993-11eb-295e-eb9f409a95fa
# ╟─03c84024-6994-11eb-1619-6b1d19725997


text:
# Implementation of the DeBoer Model

# include(joinpath(pwd(), "Tools/util.jl"))
# include(joinpath(pwd(), "Tools/tf2spec.jl"))

"""
    DeBoerModel()

A simple structure that stores a DeBoer-model.

The default initialization was taken from the original paper. When constructiong the model, single properties can be defined as keyword arguments.
The model has a state, if it was already used in `predict!`, otherwise a default is provided in `predict` or `predict!`.

# Examples
```julia
julia> model = DeBoerModel()
DeBoerModel
  a0: Int64 9
  ak: Array{Int64}((6,)) [0, 1, 2, 3, 2, 1]
  c1: Int64 -1360
  bk: Array{Int64}((6,)) [0, 2, 4, 6, 4, 2]
  TStar: Int64 3585
  γ: Float64 0.016
  c2: Float64 32.2
  c3: Float64 1.0957070684897647
  A: Int64 3
  fresp: Float64 0.3
  resp: #30 (function of type Main.CardioModels.var"#30#36"{Int64, Float64})
  noiseI: Distributions.Normal{Float64}
  noiseS: Distributions.Normal{Float64}
  operationPoint: Int64 120
  F: #31 (function of type Main.CardioModels.var"#31#37"{Int64})
  hasState: Bool false
  S: Nothing nothing
  Seff: Nothing nothing
  D: Nothing nothing
  I: Nothing nothing
  T: Nothing nothing
  ρ: Nothing nothing
  time: Int64 0
```
"""
@with_kw mutable struct DeBoerModel
    # Gains for RR from Seff
    a0::Real = 9 # vagus (fast) ms/mmHG
    ak::Vector{<:Real} = [0, 1, 2, 3, 2, 1] # sympathicus (slow) ms/mmHG
    c1::Real = 800 - 120a0 - sum(120 .* ak) # constant for mean RR interval
    # Gains for peripheral resistance
    bk::Vector{<:Real} = [0, 2, 4, 6, 4, 2]
    TStar::Real = 1425 + sum(120 .* bk) # Base aterial time constant
    # Baroreceptor Gain
    γ::Real = 0.016 # mmHg/ms
    c2::Real = 120 - 75 - 800γ
    #Windkessel
    c3::Real = 75 / 120exp(-800 / 1425)
    # Respiration
    A::Real = 3 # mmHg Amplitude of respirational influence
    fresp::Real = 0.3 # 18 breaths/min
    resp::Function = t -> A*sin(2pi * fresp * t / 1000) # Sinus Respiration
    # Noise (assumed as normally distributed with 0 mean)
    noiseI::Distribution = Normal(0, 25) # ms
    noiseS::Distribution = Normal(0, 2)  # mmHg
    # Function for effective BP Seff
    operationPoint::Real = 120
    F::Function = S -> operationPoint .+ 18atan.((S .- operationPoint) ./ 18)

    # model state
    hasState::Bool = false
    S::Union{Nothing, Vector{<:Real}} = nothing
    Seff::Union{Nothing, Vector{<:Real}} = nothing
    D::Union{Nothing, Vector{<:Real}} = nothing
    I::Union{Nothing, Vector{<:Real}} = nothing
    T::Union{Nothing, Vector{<:Real}} = nothing
    ρ::Union{Nothing, Vector{<:Real}} = nothing
    time::Real = 0
end

"""
    predict(model::DeBoerModel, n::Int; burnIn::Int = 0)

Predicts 'N' values of the cardiovascular variables for the respective model. With 'burnIn' a certain number af values can be dropped in the beginning.

# Examples
```julia
 julia> S, D, I, T, ρ = predict(model, 100)
 (S = [...], D = [...], I = [...], T = [...], ρ = [...])
 ```
 """
function predict(model::DeBoerModel, n::Int; burnIn::Int = 0)
    burnIn >= 0 || throw(DomainError(burnIn, "Burn in needs to >= 0"))
    n > 0 || throw(DomainError(n, "You need to predict at least one beat"))

    lena = length(model.ak) + 1
    lenb = length(model.bk)
    history = max(lena, lenb)
    if model.hasState && history <= length(model.S)
        # Initialization Phase - connect new arrays to state
        S = vcat(copy(model.S), Vector{Real}(undef, n + burnIn))
        Seff = vcat(copy(model.Seff), Vector{Real}(undef, n + burnIn))
        D = vcat(copy(model.D), Vector{Real}(undef, n + burnIn))
        I = vcat(copy(model.I), Vector{Real}(undef, n + burnIn))
        T = vcat(copy(model.T), Vector{Real}(undef, n + burnIn))
        ρ = vcat(copy(model.ρ), Vector{Real}(undef, n + burnIn))
        time = model.time
        burnIn += history
    else
        # Default Initialization
        if model.hasState
            @warn "State was too short (length of $(length(model.S)) but history needed is $history), switching to default initialization"
        end
        burnIn += history
        S = Vector{Real}(undef, n + burnIn)
        Seff = Vector{Real}(undef, n + burnIn)
        D = Vector{Real}(undef, n + burnIn)
        I = Vector{Real}(undef, n + burnIn)
        T = Vector{Real}(undef, n + burnIn)
        ρ = Vector{Real}(undef, n + burnIn)
        time = 0
        # fill in starting data
        for i in 1:history
            S[i] = model.operationPoint
            Seff[i] = model.F(model.operationPoint)
            D[i] = 75
            I[i] = 800
            T[i] = 1425
        end
    end

    # generate future data
    for i in history+1:n+burnIn
        time += I[i-1]
        D[i] = model.c3 * S[i-1] * exp(-I[i-1]/T[i-1])
        ρ[i] = model.resp(time)
        S[i] = D[i] + model.γ * I[i-1] + model.c2 + ρ[i] + rand(model.noiseS)
        Seff[i] = model.F(S[i])
        I[i] = model.c1 + model.a0 * Seff[i] + sum(model.ak .* reverse(Seff[i-(lena-1):i-1])) + rand(model.noiseI)
        T[i] = model.TStar - sum(model.bk .* reverse(Seff[i-lenb:i-1]))
    end
    return (S = S[burnIn+1:end], D = D[burnIn+1:end], I = I[burnIn+1:end], T = T[burnIn+1:end], ρ = ρ[burnIn+1:end])
end

# Statefull predict function
"""
    predict!(model::DeBoerModel, n::Int; burnIn::Int = 0)

Updates the state of the state of the model with the last predicted values.
Predicts 'N' values of the cardiovascular variables for the respective model and updates its state. With 'burnIn' a certain number af values can be dropped in the beginning.

# Examples
```julia
 julia> S, D, I, T, ρ = predict!(model, 100)
 (S = [...], D = [...], I = [...], T = [...], ρ = [...])
 julia> model.hasState
 true
 ```
"""
function predict!(model::DeBoerModel, n::Int; burnIn::Int = 0)
    S, D, I, T, ρ = predict(model, n + burnIn; burnIn = 0)

    lena = length(model.ak) + 1
    lenb = length(model.bk)
    history = max(lena, lenb)

    if n >= history
        # update / shorten  state
        model.hasState = true
        model.S = S[end-history+1:end]
        model.Seff = model.F(model.S)
        model.D = D[end-history+1:end]
        model.I = I[end-history+1:end]
        model.T = T[end-history+1:end]
        model.ρ = ρ[end-history+1:end]
        model.time += sum(I)
    elseif model.hasState
        model.hasState = true
        model.S = append!(model.S, S)[end-history+1:end]
        model.Seff = model.F(model.S)
        model.D = append!(model.D, D)[end-history+1:end]

        model.I = append!(model.I, I)[end-history+1:end]
        model.T = append!(model.T, T)[end-history+1:end]
        model.ρ = append!(model.ρ, ρ)[end-history+1:end]
        model.time += sum(I)
    else
        @warn "n was to short to update a valid state, it remains unchanged."
    end

    return (S = S[burnIn+1:end], D = D[burnIn+1:end], I = I[burnIn+1:end], T = T[burnIn+1:end], ρ = ρ[burnIn+1:end])
end

# simulates a phenylephrine injection by increaing the peripheral resistance over a defined time
function phenylephrine(model::DeBoerModel, TStar_increase::Real, n_increase::Int)
    n_increase > 0 || throw(DomainError(n_increase, "peripheral resistance needs to invrease over a minimum of one beat"))
    TStar_step = TStar_increase / n_increase

    S = Vector{Float64}(undef, n_increase)
    D = Vector{Float64}(undef, n_increase)
    I = Vector{Float64}(undef, n_increase)
    T = Vector{Float64}(undef, n_increase)
    ρ = Vector{Float64}(undef, n_increase)
    for i in 1:n_increase
        model.TStar += TStar_step
        res = predict!(model, 1)
        S[i] = res[:S][1]
        D[i] = res[:D][1]
        I[i] = res[:I][1]
        T[i] = res[:T][1]
        ρ[i] = res[:ρ][1]
    end
    return (S = S, D = D, I = I, T = T, ρ = ρ)
end

function BRSf(model::DeBoerModel)
    a = [model.a0; model.ak]
    A = impRes2tf(a)
    w = range(0, stop = π, length = 256)
    mag, phase, w = bode(A, w)
    f = w ./ 2π
    return f, vec(mag), vec(phase), A
end