text:
# For stochastic Runge-Kutta methods (SIRK and WIRK)
# x - the solution vector to be updated
# Vx - the matrix containing the drift vector evaluated at the internal stages v(Q_i) (SIRK) or v(Q0_i) (WIRK)
# Bx - the array containing the diffusion matrix evaluated at the internal stages B(Q_i) (SIRK) or B(Q1^(l)_i) (WIRK)
# bdrift - the Runge-Kutta coefficients for the drift part
# bdiff - the Runge-Kutta coefficients for the diffusion part
# Δt - the time step
# ΔW - the increments of the Brownian motion (SFIRK) or the increments represented by the random variables \hat I^(k) (WFIRK)
function update_solution!(x::Union{Vector{T}, Vector{TwicePrecision{T}}}, Vx::Vector{Vector{T}}, Bx::Vector{Matrix{T}}, bdrift::Vector{T}, bdiff::Vector{T}, Δt::T, ΔW::Vector{T}) where {T}

    @assert length(bdrift) == length(bdiff) == length(Vx) == length(Bx)

    for i in eachindex(Vx, Bx)
        @assert length(x) == length(Vx[i]) == size(Bx[i], 1)
        @assert length(ΔW)== size(Bx[i], 2)
    end

    local Δx::eltype(x)

    # Contribution from the drift part
    for k in eachindex(x)
        Δx = 0.
        for i in eachindex(Vx)
            Δx += bdrift[i] * Vx[i][k]
        end
        x[k] += Δt * Δx
    end

    local Δy = zero(ΔW)

    # Contribution from the diffusion part
    for k in eachindex(x)
        Δy .= 0
        for i in eachindex(Bx)
            Δy .+= bdiff[i] .* Bx[i][k,:]
        end
        x[k] += dot(Δy,ΔW)
    end
end

# For stochastic Runge-Kutta methods
# x - the solution vector to be updated
# Vx - the matrix containing the drift vector evaluated at the internal stages v(Q_i)
# Bx - the array containing the diffusion matrix evaluated at the internal stages B(Q_i)
# bdrift - the Runge-Kutta coefficients for the drift part
# bdiff - the Runge-Kutta coefficients for the ΔW terms of the diffusion part
# bdiff2- the Runge-Kutta coefficients for the ΔZ terms of the diffusion part
# Δt - the time step
# ΔW - the increments of the Brownian motion
# ΔZ - the integrals of the increments of the Brownian motion
function update_solution!(x::Union{Vector{T}, Vector{TwicePrecision{T}}}, Vx::Vector{Vector{T}}, Bx::Vector{Matrix{T}}, bdrift::Vector{T}, bdiff::Vector{T}, bdiff2::Vector{T}, Δt::T, ΔW::Vector{T}, ΔZ::Vector{T}) where {T}

    @assert length(bdrift) == length(bdiff) == length(bdiff2) == length(Vx) == length(Bx)

    for i in eachindex(Vx, Bx)
        @assert length(x) == length(Vx[i]) == size(Bx[i], 1)
        @assert length(ΔW) == length(ΔZ) == size(Bx[i], 2)
    end

    local Δx::eltype(x)

    # Contribution from the drift part
    for k in eachindex(x)
        Δx = 0.
        for i in eachindex(Vx)
            Δx += bdrift[i] * Vx[i][k]
        end
        x[k] += Δt * Δx
    end

    local Δy = zero(ΔW)

    # Contribution from the diffusion part (ΔW terms)
    for k in eachindex(x)
        Δy .= 0
        for i in eachindex(Bx)
            for l in eachindex(Δy)
                Δy[l] += bdiff[i] * Bx[i][k,l]
            end
        end
        x[k] += dot(Δy,ΔW)
    end

    # Contribution from the diffusion part (ΔZ terms)
    for k in eachindex(x)
        Δy .= 0
        for i in eachindex(Bx)
            for l in eachindex(Δy)
                Δy[l] += bdiff2[i] * Bx[i][k,l]
            end
        end
        x[k] += dot(Δy,ΔZ)/Δt
    end
end

# For stochastic partitioned Runge-Kutta methods
# q, p - the solution vector to be updated
# Vqp, Fqp - the matrix containing the drift vectors evaluated at the internal stages v(Q_i), f(Q_i)
# Bqp, Gqp - the array containing the diffusion matrices evaluated at the internal stages B(Q_i), G(Q_i)
# bqdrift, bpdrift - the Runge-Kutta coefficients for the drift parts of the q and p equations
# bqdiff, bpdiff - the Runge-Kutta coefficients for the diffusion parts of the q and p equations
# Δt - the time step
# ΔW - the increments of the Brownian motion
function update_solution!(q::Vector{VT}, p::Vector{VT},
                          Vqp::Vector{Vector{T}}, Fqp::Vector{Vector{T}},
                          Bqp::Vector{Matrix{T}}, Gqp::Vector{Matrix{T}},
                          bqdrift::Vector{T}, bqdiff::Vector{T},
                          bpdrift::Vector{T}, bpdiff::Vector{T},
                          Δt::T, ΔW::Vector{T}) where {T, VT <: Union{T, TwicePrecision{T}}}

    @assert length(bqdrift) == length(bqdiff) == length(bpdrift) == length(bpdiff) == length(Vqp) == length(Fqp) == length(Bqp) == length(Gqp)

    for i in eachindex(Vqp, Fqp, Bqp, Gqp)
        @assert length(q) == length(p) == length(Vqp[i]) == length(Fqp[i]) == size(Bqp[i], 1) == size(Gqp[i], 1)
        @assert length(ΔW) == size(Bqp[i], 2) == size(Gqp[i], 2)
    end

    local Δq::eltype(q)
    local Δp::eltype(p)

    # Contribution from the drift part
    for k in eachindex(q, p)
        Δq = 0.
        Δp = 0.
        for i in eachindex(Vqp, Fqp)
            Δq += bqdrift[i] * Vqp[i][k]
            Δp += bpdrift[i] * Fqp[i][k]
        end
        q[k] += Δt * Δq
        p[k] += Δt * Δp
    end

    local Δy = zero(ΔW)
    local Δz = zero(ΔW)

    # Contribution from the diffusion part
    for k in eachindex(q, p)
        Δy .= 0
        Δz .= 0

        for i in eachindex(Bqp, Gqp)
            Δy .+= bqdiff[i] .* Bqp[i][k,:]
            Δz .+= bpdiff[i] .* Gqp[i][k,:]
        end

        q[k] += dot(Δy,ΔW)
        p[k] += dot(Δz,ΔW)
    end
end

# For stochastic split partitioned Runge-Kutta methods
# q, p - the solution vector to be updated
# Vqp, Fqp1, Fqp2 - the matrix containing the drift vectors evaluated at the internal stages v(Q_i), fi(Q_i)
# Bqp, Gqp1, Gqp2 - the array containing the diffusion matrices evaluated at the internal stages B(Q_i), Gi(Q_i)
# bqdrift, bpdrift1, bpdrift2 - the Runge-Kutta coefficients for the drift parts of the q and p equations
# bqdiff, bpdiff1, bpdiff2 - the Runge-Kutta coefficients for the diffusion parts of the q and p equations
# Δt - the time step
# ΔW - the increments of the Brownian motion
function update_solution!(q::Union{Vector{T}, Vector{TwicePrecision{T}}}, p::Union{Vector{T}, Vector{TwicePrecision{T}}},
                            Vqp::Vector{Vector{T}}, Fqp1::Vector{Vector{T}}, Fqp2::Vector{Vector{T}},
                            Bqp::Vector{Matrix{T}}, Gqp1::Vector{Matrix{T}}, Gqp2::Vector{Matrix{T}},
                            bqdrift::Vector{T}, bqdiff::Vector{T},
                            bpdrift1::Vector{T}, bpdrift2::Vector{T},
                            bpdiff1::Vector{T}, bpdiff2::Vector{T},
                            Δt::T, ΔW::Vector{T}) where {T}

   @assert length(bqdrift) == length(bqdiff) == length(bpdrift1) == length(bpdrift2) == length(bpdiff1) == length(bpdiff2) == length(Vqp) == length(Fqp1) == length(Fqp2) == length(Bqp) == length(Gqp1) == length(Gqp2)

   for i in eachindex(Vqp, Fqp1, Fqp2, Bqp, Gqp1, Gqp2)
       @assert length(q) == length(p) == length(Vqp[i]) == length(Fqp1[i]) == length(Fqp2[i]) == size(Bqp[i], 1) == size(Gqp1[i], 1) == size(Gqp2[i], 1)
       @assert length(ΔW) == size(Bqp[i], 2) == size(Gqp1[i], 2) == size(Gqp2[i], 2)
   end

   local Δq::eltype(q)
   local Δp::eltype(p)

   # Contribution from the drift part
   for k in eachindex(q, p)
       Δq = 0.
       Δp = 0.
       for i in eachindex(Vqp, Fqp1, Fqp2)
           Δq += bqdrift[i] * Vqp[i][k]
           Δp += bpdrift1[i] * Fqp1[i][k] + bpdrift2[i] * Fqp2[i][k]
       end
       q[k] += Δt * Δq
       p[k] += Δt * Δp
   end

   local Δy = zero(ΔW)
   local Δz = zero(ΔW)

   # Contribution from the diffusion part
   for k in eachindex(q, p)
       Δy .= 0
       Δz .= 0

       for i in eachindex(Bqp, Gqp1, Gqp2)
           Δy .+= bqdiff[i]  .* Bqp[i][k,:]
           Δz .+= bpdiff1[i] .* Gqp1[i][k,:] .+ bpdiff2[i] .* Gqp2[i][k,:]
       end

       q[k] += dot(Δy,ΔW)
       p[k] += dot(Δz,ΔW)
   end
end

# For weak Runge-Kutta methods WERK
# x - the solution vector to be updated
# Vx - the matrix containing the drift vector evaluated at the internal stages v(Q_i)
# Bx1 - the array containing the diffusion matrix evaluated at the internal stages H^(l)_i, such that Bx1[:,l,i] is evaluated at H^(l)_i
# Bx2 - the array containing the diffusion matrix evaluated at the internal stages \hat H^(l)_i, such that Bx2[:,l,i] is evaluated at \hat H^(l)_i
# alpha - the Runge-Kutta coefficients for the drift part
# beta1 - the Runge-Kutta coefficients for the diffusion term with the random increments
# beta2 - the Runge-Kutta coefficients for the second diffusion term
# Δt - the time step
# ΔW - the increments of the Brownian motion represented by the random variables \hat I^(k)
function update_solution!(x::Union{Vector{T}, Vector{TwicePrecision{T}}}, Vx::Vector{Vector{T}}, Bx1::Vector{Matrix{T}}, Bx2::Vector{Matrix{T}}, alpha::Vector{T}, beta1::Vector{T}, beta2::Vector{T}, Δt::T, ΔW::Vector{T}) where {T}
    @assert length(alpha) == length(beta1) == length(beta2) == length(Vx) == length(Bx1) == length(Bx2)

    for i in eachindex(Vx, Bx1, Bx2)
        @assert length(x) == length(Vx[i]) == size(Bx1[i], 1) == size(Bx2[i], 1)
        @assert length(ΔW)== size(Bx1[i], 2) == size(Bx2[i], 2)
    end

    local Δx::eltype(x)

    # Contribution from the drift part
    for k in eachindex(x)
        Δx = 0.
        for i in eachindex(Vx)
            Δx += alpha[i] * Vx[i][k]
        end
        x[k] += Δt * Δx
    end

    local Δy = zero(ΔW)

    # Contribution from the diffusion term with the random variables I^(k)_i
    for k in eachindex(x)
        Δy .= 0
        for i in eachindex(Bx1)
            Δy .+= beta1[i] .* Bx1[i][k,:]
        end
        x[k] += dot(Δy,ΔW)
    end

    # Contribution from the second diffusion term
    for k in eachindex(x)
        Δx= 0.
        for i in eachindex(Bx2)
            for l in axes(Bx2[i], 2)
                Δx += beta2[i] * Bx2[i][k,l]
            end
        end
        x[k] += sqrt(Δt)*Δx
    end
end

# For stochastic Runge-Kutta methods
function update_solution!(x::Union{Vector{T}, Vector{TwicePrecision{T}}}, Vx::Vector{Vector{T}}, Bx::Vector{Matrix{T}},
                            bdrift::Vector{T}, b̂drift::Vector, bdiff::Vector{T}, b̂diff::Vector, Δt::T, ΔW::Vector{T}) where {T}
    update_solution!(x, Vx, Bx, bdrift, bdiff, Δt, ΔW)
    update_solution!(x, Vx, Bx, b̂drift, b̂diff, Δt, ΔW)
end

# For stochastic partitioned Runge-Kutta methods
function update_solution!(q::Union{Vector{T}, Vector{TwicePrecision{T}}}, p::Union{Vector{T}, Vector{TwicePrecision{T}}},
                            Vqp::Vector{Vector{T}}, Fqp::Vector{Vector{T}}, Bqp::Vector{Matrix{T}}, Gqp::Vector{Matrix{T}},
                            bqdrift::Vector{T}, b̂qdrift::Vector, bqdiff::Vector{T}, b̂qdiff::Vector,
                            bpdrift::Vector{T}, b̂pdrift::Vector, bpdiff::Vector{T}, b̂pdiff::Vector, Δt::T, ΔW::Vector{T}) where {T}
    update_solution!(q, p, Vqp, Fqp, Bqp, Gqp, bqdrift, bqdiff, bpdrift, bpdiff, Δt, ΔW)
    update_solution!(q, p, Vqp, Fqp, Bqp, Gqp, b̂qdrift, b̂qdiff, b̂pdrift, b̂pdiff, Δt, ΔW)
end

# For stochastic split partitioned Runge-Kutta methods
function update_solution!(q::Union{Vector{T}, Vector{TwicePrecision{T}}}, p::Union{Vector{T}, Vector{TwicePrecision{T}}},
                            Vqp::Vector{Vector{T}}, Fqp1::Vector{Vector{T}}, Fqp2::Vector{Vector{T}},
                            Bqp::Vector{Matrix{T}}, Gqp1::Vector{Matrix{T}}, Gqp2::Vector{Matrix{T}},
                            bqdrift::Vector{T}, b̂qdrift::Vector, bqdiff::Vector{T}, b̂qdiff::Vector,
                            bpdrift1::Vector{T}, b̂pdrift1::Vector, bpdrift2::Vector{T}, b̂pdrift2::Vector,
                            bpdiff1::Vector{T}, b̂pdiff1::Vector, bpdiff2::Vector{T}, b̂pdiff2::Vector,
                            Δt::T, ΔW::Vector{T}) where {T}
    update_solution!(q, p, Vqp, Fqp1, Fqp2, Bqp, Gqp1, Gqp2, bqdrift, bqdiff, bpdrift1, bpdrift2, bpdiff1, bpdiff2, Δt, ΔW)
    update_solution!(q, p, Vqp, Fqp1, Fqp2, Bqp, Gqp1, Gqp2, b̂qdrift, b̂qdiff, b̂pdrift1, b̂pdrift2, b̂pdiff1, b̂pdiff2, Δt, ΔW)
end

 # For weak Runge-Kutta methods WERK
 function update_solution!(x::Union{Vector{T}, Vector{TwicePrecision{T}}}, Vx::Vector{Vector{T}}, Bx1::Vector{Matrix{T}}, Bx2::Vector{Matrix{T}},
                             alpha::Vector{T}, âlpha::Vector, beta1::Vector{T}, b̂eta1::Vector, beta2::Vector{T}, b̂eta2::Vector, Δt::T, ΔW::Vector{T}) where {T}
     update_solution!(x, Vx, Bx1, Bx2, alpha, beta1, beta2, Δt, ΔW)
     update_solution!(x, Vx, Bx1, Bx2, âlpha, b̂eta1, b̂eta2, Δt, ΔW)
 end


text:
abstract type AbstractLatticeIterator end
# All indices are sites
abstract type DirectLatticeIterator <: AbstractLatticeIterator end
# first index is a meta index (e.g. direction), rest for sites
abstract type DeferredLatticeIterator <: AbstractLatticeIterator end
# Wraps a lattice iterator to do change what happens with the output
abstract type LatticeIterationWrapper{LI <: AbstractLatticeIterator} <: AbstractLatticeIterator end

################################################################################
### EachSiteAndFlavor
################################################################################

"""
    EachSiteAndFlavor(mc, model)

Creates an iterator which iterates through the diagonal of the Greensfunctio
"""
struct EachSiteAndFlavor <: DirectLatticeIterator
    N::Int64
end
function EachSiteAndFlavor(mc::MonteCarloFlavor, model::Model)
    EachSiteAndFlavor(length(lattice(model)) * nflavors(model))
end

@bm Base.iterate(iter::EachSiteAndFlavor, i=1) = i ≤ iter.N ? (i, i+1) : nothing
Base.length(iter::EachSiteAndFlavor) = iter.N
Base.eltype(::EachSiteAndFlavor) = Int64

################################################################################
### EachSite
################################################################################

"""
    EachSite(Nsites)
    EachSite(lattice)
    EachSite(mc, model)

Creates an iterator which iterates through every site of a given lattice.
"""
struct EachSite <: DirectLatticeIterator
    N::Int64
end
EachSite(l::AbstractLattice) = EachSite(length(l))
EachSite(mc::MonteCarloFlavor, model::Model) = EachSite(lattice(model))
eachsite(args...) = EachSite(args...)

@bm Base.iterate(iter::EachSite, i=1) = i ≤ iter.N ? (i, i+1) : nothing
Base.length(iter::EachSite) = iter.N
Base.eltype(::EachSite) = Int64

################################################################################
### OnSite
################################################################################

"""
    OnSite(Nsites)
    OnSite(lattice)
    OnSite(mc, model)

Creates an iterator which iterates through every site of a given lattice,
returning (site, site) at every step.
"""
struct OnSite <: DirectLatticeIterator
    N::Int64
end
OnSite(l::AbstractLattice) = OnSite(length(l))
OnSite(mc::MonteCarloFlavor, model::Model) = OnSite(lattice(model))
onsite(args...) = OnSite(args...)

@bm Base.iterate(iter::OnSite, i=1) = i ≤ iter.N ? ((i, i), i+1) : nothing
Base.length(iter::OnSite) = iter.N
Base.eltype(::OnSite) = Tuple{Int64, Int64}

################################################################################
### EachSitePair
################################################################################

"""
    EachSitePair(Nsites)
    EachSitePair(lattice)
    EachSitePair(mc, model)

Creates an iterator which returns every pair of sites `(s1, s2)` with
`s1, s2 ∈ 1:Nsites`.
"""
struct EachSitePair <: DirectLatticeIterator
    N::Int64
end
EachSitePair(l::AbstractLattice) = EachSitePair(length(l))
EachSitePair(mc::MonteCarloFlavor, model::Model) = EachSitePair(lattice(model))
eachsitepair(args...) = EachSitePair(args...)

@bm function Base.iterate(iter::EachSitePair, i=1)
    if i ≤ iter.N^2
        return ((div(i-1, iter.N)+1, mod1(i, iter.N)), i+1)
    else
        return nothing
    end
end
Base.length(iter::EachSitePair) = iter.N^2
Base.eltype(::EachSitePair) = NTuple{2, Int64}

################################################################################
### EachSitePairByDistance
################################################################################

"""
    EachSitePairByDistance(lattice)
    EachSitePairByDistance(mc, model)

Creates an iterator which returns triplets `(direction index, source, target)`
sorted by distance. The `direction index` identifies each unqiue direction
`position(target) - position(source)`.

Requires `lattice` to implement `positions` and `lattice_vectors`.
"""
struct EachSitePairByDistance <: DeferredLatticeIterator
    N::Int64
    pairs::Vector{Vector{Tuple{Int64, Int64}}}
end

# For shifting sites across periodic bounds
function generate_combinations(vs::Vector{<: Vector})
    out = [zeros(length(vs[1]))]
    for v in vs
        out = vcat([e.-v for e in out], out, [e.+v for e in out])
    end
    out
end

# norm + ϵ * angle(v, e_x)
function directed_norm(v, ϵ)
    l = norm(v)
    if length(v) == 2 && l > ϵ
        angle = acos(dot([1, 0], v) / l)
        v[2] < 0 && (angle = 2pi - angle)
        return l + ϵ * angle
    else
        return l
    end
end

function EachSitePairByDistance(lattice::AbstractLattice, ϵ = 1e-6)
    _positions = positions(lattice)
    wrap = generate_combinations(lattice_vectors(lattice))
    directions = Vector{Float64}[]
    # (src, trg), first index is dir, second index irrelevant
    bonds = [Tuple{Int64, Int64}[] for _ in 1:length(lattice)]

    for origin in 1:length(lattice)
        for (trg, p) in enumerate(_positions)
            d = _positions[origin] .- p .+ wrap[1]
            for v in wrap[2:end]
                new_d = _positions[origin] .- p .+ v
                if directed_norm(new_d, ϵ) + ϵ < directed_norm(d, ϵ)
                    d .= new_d
                end
            end
            # The rounding will allow us to use == here
            idx = findfirst(dir -> isapprox(dir, d, atol=ϵ), directions)
            if idx === nothing
                push!(directions, d)
                if length(bonds) < length(directions)
                    push!(bonds, Tuple{Int64, Int64}[])
                end
                push!(bonds[length(directions)], (origin, trg))
            else
                push!(bonds[idx], (origin, trg))
            end
        end
    end

    temp = sortperm(directions, by = v -> directed_norm(v, ϵ))

    EachSitePairByDistance(length(lattice)^2, bonds[temp])
end
function EachSitePairByDistance(mc::MonteCarloFlavor, model::Model)
    EachSitePairByDistance(lattice(model))
end

@bm function Base.iterate(iter::EachSitePairByDistance, state = (1, 1))
    dir, i = state
    # next_dir = dir + div(i, length(iter.pairs[dir]))
    # next_i = mod1(i+1, length(iter.pairs[dir]))
    # if next_dir > length(iter.pairs)
    #     return nothing
    # else
    #     src, trg = iter.pairs[dir][i]
    #     return (dir, src, trg), (next_dir, next_i)
    # end

    if dir ≤ length(iter.pairs)
        if i ≤ length(iter.pairs[dir])
            src, trg = iter.pairs[dir][i]
            return (dir, src, trg), (dir, i+1)
        else
            return iterate(iter, (dir+1, 1))
        end
    else
        return nothing
    end
end
ndirections(iter::EachSitePairByDistance) = size(iter.pairs)
Base.length(iter::EachSitePairByDistance) = iter.N
Base.eltype(::EachSitePairByDistance) = NTuple{3, Int64}

"""
    in_direction(iter::EachSitePairByDistance, dir::Int64)

Returns a vector of each `(src, trg)` pair in the given direction.
"""
in_direction(iter::EachSitePairByDistance, dir::Int64) = iter.pairs[dir]

################################################################################
### EachLocalQuadByDistance
################################################################################

"""
    EachLocalQuadByDistance{K}(lattice)
    EachLocalQuadByDistance{K}(mc, model)

Creates an iterator which returns a tuple
`(dir12, dir1, dir2, src1, trg1, src2, trg2)` with each iteration. The
directional index `dir12` identifies the direction `src2 - src1`.
The type parameter `K` is the maximum direction index `dir1` (`dir2`) a pair
`(src1, trg1)` (`(src2, trg2)`) may have. Picking for example `K = 5` for a
`SquareLattice` will result in the `src` site (`K = 1`) and the four nearest
neighbors of the `src` site being selected as targets.
Note that `K` is not always the number of nearest neighbors + 1. For example,
the Honeycomb lattice has 3 NNs, but 6 unique directions. To include all NN
here, one should pick `K = 7`.

The iterator is sorted by distance `src2 - src1`.

Requires `lattice` to implement `positions` and `lattice_vectors`.
"""
struct EachLocalQuadByDistance{K} <: DeferredLatticeIterator
    N::Int64
    implied::Array{Vector{NTuple{4, UInt16}}, 3}
end

function EachLocalQuadByDistance{K}(lattice::AbstractLattice) where {K}
    # (src1, src2) pairs from dir
    pairs_by_dir = EachSitePairByDistance(lattice)

    # (dir, trg1) from src1
    trg_from_src = [Tuple{Int64, Int64}[] for _ in 1:length(lattice)]
    for dir in 1:K
        for (src, trg) in in_direction(pairs_by_dir, dir)
            push!(trg_from_src[src], (dir, trg))
        end
    end

    # NOTE: Creating 7-Tuples is faster ~15% faster but I'm worried about memory
    # usuage. Caching 4-tuples is still ~36% faster than running the loop below
    # during measurements, and is not quite as demanding on memory.
    total_length = sum(length(x) for x in trg_from_src)^2
    # linear = Vector{NTuple{7, UInt16}}(undef, total_length)
    # k = 1
    implied = [NTuple{4, UInt16}[] for _ in 1:length(pairs_by_dir.pairs), _ in 1:K, _ in 1:K]
    dir12, idx, i, j = (1,1,1,1)
    while true
        if dir12 ≤ length(pairs_by_dir.pairs)
            if idx ≤ length(pairs_by_dir.pairs[dir12])
                src1, src2 = pairs_by_dir.pairs[dir12][idx]
                if i ≤ length(trg_from_src[src1])
                    dir1, trg1 = trg_from_src[src1][i]
                    if j ≤ length(trg_from_src[src2])
                        dir2, trg2 = trg_from_src[src2][j]
                        # linear[k] = UInt16.((dir12, dir1, dir2, src1, trg1, src2, trg2))
                        push!(
                            implied[dir12, dir1, dir2],
                            UInt16.((src1, trg1, src2, trg2))
                        )
                        # k += 1
                        j += 1
                    else
                        i += 1
                        j = 1
                    end
                else
                    idx += 1
                    i = j = 1
                end

            else
                dir12 += 1
                idx = i = j = 1
            end
        else
            break
        end
    end

    EachLocalQuadByDistance{K}(total_length, implied)
end
function EachLocalQuadByDistance{K}(mc::MonteCarloFlavor, model::Model) where {K}
    EachLocalQuadByDistance{K}(lattice(model))
end

# @bm function Base.iterate(iter::EachLocalQuadByDistance, state=1)
#     if state <= length(iter.direct)
#         return (iter.direct[state], state+1)
#     else
#         return nothing
#     end
# end
@bm function Base.iterate(iter::EachLocalQuadByDistance{K}, state=(1,1)) where {K}
    i, j = state
    if i <= length(iter.implied)
        N = length(iter.implied[i])

        # This is required for lattices with a basis
        j = ifelse(N == 0, 1, j)
        while N == 0
            i += 1
            i <= length(iter.implied) || return nothing
            N = length(iter.implied[i])
        end

        next_j = mod1(j+1, N)
        next_i = i + div(j, N)

        t = iter.implied[i][j]
        return ((i, t[1], t[2], t[3], t[4]), (next_i, next_j))
    else
        return nothing
    end
end

ndirections(iter::EachLocalQuadByDistance{K}) where {K} = size(iter.implied)
Base.length(iter::EachLocalQuadByDistance) = iter.N
Base.eltype(::EachLocalQuadByDistance) = Tuple{Int64, UInt16, UInt16, UInt16, UInt16}

################################################################################
### EachLocalQuadBySyncedDistance
################################################################################

"""
    EachLocalQuadBySyncedDistance{K}(lattice)
    EachLocalQuadBySyncedDistance{K}(mc, model)

Creates an iterator which returns a tuple
`(dir12, dir_ii, src1, trg1, src2, trg2)` with each iteration. The
directional index `dir12` identifies the direction `src2 - src1`.
The type parameter `K` is the maximum direction index `dir_ii` a pair
`(src1, trg1)` (`(src2, trg2)`) may have. Picking for example `K = 5` for a
`SquareLattice` will result in the `src` site (`K = 1`) and the four nearest
neighbors of the `src` site being selected as targets.
Note that `K` is not always the number of nearest neighbors + 1. For example,
the Honeycomb lattice has 3 NNs, but 6 unique directions. To include all NN
here, one should pick `K = 7`.

The iterator is sorted by distance `src2 - src1`.

Requires `lattice` to implement `positions` and `lattice_vectors`.
"""
struct EachLocalQuadBySyncedDistance{K} <: DeferredLatticeIterator
    N::Int64
    implied::Array{Vector{NTuple{4, UInt16}}, 2}
end

function EachLocalQuadBySyncedDistance{K}(lattice::AbstractLattice) where {K}
    # (src1, src2) pairs from dir
    pairs_by_dir = EachSitePairByDistance(lattice)

    # (dir, trg1) from src1
    trg_from_src = [Tuple{Int64, Int64}[] for _ in 1:length(lattice)]
    for dir in 1:K
        for (src, trg) in in_direction(pairs_by_dir, dir)
            push!(trg_from_src[src], (dir, trg))
        end
    end

    implied = [NTuple{4, UInt16}[] for _ in 1:length(pairs_by_dir.pairs), _ in 1:K]
    dir12, idx, i, j = (1,1,1,1)
    while true
        if dir12 ≤ length(pairs_by_dir.pairs)
            if idx ≤ length(pairs_by_dir.pairs[dir12])
                src1, src2 = pairs_by_dir.pairs[dir12][idx]
                if i ≤ length(trg_from_src[src1])
                    dir1, trg1 = trg_from_src[src1][i]
                    if j ≤ length(trg_from_src[src2])
                        dir2, trg2 = trg_from_src[src2][j]
                        j += 1
                        dir1 != dir2 && continue
                        push!(
                            implied[dir12, dir1],
                            UInt16.((src1, trg1, src2, trg2))
                        )
                    else
                        i += 1
                        j = 1
                    end
                else
                    idx += 1
                    i = j = 1
                end

            else
                dir12 += 1
                idx = i = j = 1
            end
        else
            break
        end
    end
    total_length = sum(length(x) for x in implied)

    EachLocalQuadBySyncedDistance{K}(total_length, implied)
end
function EachLocalQuadBySyncedDistance{K}(mc::MonteCarloFlavor, model::Model) where {K}
    EachLocalQuadBySyncedDistance{K}(lattice(model))
end

@bm function Base.iterate(iter::EachLocalQuadBySyncedDistance{K}, state = (1, 1)) where {K}
    i, j = state
    if i <= length(iter.implied)
        N = length(iter.implied[i])

        # This is required for lattices with a basis
        j = ifelse(N == 0, 1, j)
        while N == 0
            i += 1
            i <= length(iter.implied) || return nothing
            N = length(iter.implied[i])
        end

        next_j = mod1(j+1, N)
        next_i = i + div(j, N)
        t = iter.implied[i][j]
        return ((i, t[1], t[2], t[3], t[4]), (next_i, next_j))
    else
        return nothing
    end
end
ndirections(iter::EachLocalQuadBySyncedDistance{K}) where {K} = size(iter.implied)
Base.length(iter::EachLocalQuadBySyncedDistance) = iter.N
Base.eltype(::EachLocalQuadBySyncedDistance) = Tuple{Int64, UInt16, UInt16, UInt16, UInt16}

################################################################################
### Additonal stuff
################################################################################

# function directions(::EachSitePair, lattice::AbstractLattice)
#     pos = positions(lattice)
#     [p2 .- p1 for p2 in pos for p1 in pos]
# end

function directions(iter::EachSitePairByDistance, lattice::AbstractLattice, ϵ=1e-6)
    pos = MonteCarlo.positions(lattice)
    wrap = generate_combinations(lattice_vectors(lattice))

    dirs = map(iter.pairs) do pairs
        src, trg = pairs[1]
        _d = pos[src] - pos[trg]
        # Find lowest distance w/ periodic bounds
        d = _d .+ wrap[1]
        for v in wrap[2:end]
            new_d = _d .+ v
            if directed_norm(new_d, ϵ) + ϵ < directed_norm(d, ϵ)
                d .= new_d
            end
        end
        d
    end
end

directions(dqmc::MonteCarloFlavor, ϵ=1e-6) = directions(lattice(dqmc), ϵ)
directions(model::Model, ϵ=1e-6) = directions(lattice(model), ϵ)
function directions(lattice::AbstractLattice, ϵ = 1e-6)
    _positions = positions(lattice)
    wrap = generate_combinations(lattice_vectors(lattice))
    directions = Vector{Float64}[]
    for origin in 1:length(lattice)
        for (trg, p) in enumerate(_positions)
            d = _positions[origin] .- p .+ wrap[1]
            for v in wrap[2:end]
                new_d = _positions[origin] .- p .+ v
                if directed_norm(new_d, ϵ) + ϵ < directed_norm(d, ϵ)
                    d .= new_d
                end
            end
            idx = findfirst(dir -> isapprox(dir, d, atol=ϵ), directions)
            if idx === nothing
                push!(directions, d)
            end
        end
    end
    # temp = sortperm(directions, by=norm)
    # directions[temp]
    sort!(directions, by = v -> directed_norm(v, ϵ))
end

################################################################################
### Sum Wrapper
################################################################################

struct Sum{LI} <: LatticeIterationWrapper{LI}
    iter::LI
    function Sum{LI}(args...; kwargs...) where {LI <: AbstractLatticeIterator}
        new{LI}(LI(args...; kwargs...))
    end
    function Sum(iter::LI) where {LI <: AbstractLatticeIterator}
        new{LI}(iter)
    end
end

Base.iterate(s::Sum) = iterate(s.iter)
Base.iterate(s::Sum, state) = iterate(s.iter, state)
Base.length(s::Sum) = length(s.iter)
Base.eltype(s::Sum) = eltype(s.iter)

################################################################################
### Symmetry Wrapper
################################################################################

# This is a weird monster :(
# The idea is that we construct a thin wrapper
# li = ApplySymmetries{EachLocalQuadByDistance{5}}(sym1, sym2, ...)
# which contains weights for different neighbors in the given symmetries
# The backend then bundles all of these and constructs one thick wrapper by
# calling
# li(dqmc, model)
# This then actually contains the udnerlying lattice iterator

"""
    ApplySymmetries{lattice_iterator_type}(symmetries...)

`ApplySymmetries` is a wrapper for a `DeferredLatticeIterator`. It is meant to
specify how results from different directions are to be added up.

For example `T = EachLocalQuadByDistance{5}` specifies 4 site tuples where two
sites have set distances between them - one of the first 5 smallest ones. In a
square lattice these would be on-site (1) and the four nearest neighbors (2-5).
We may use `iter = ApplySymmetries{}([1], [0, 1, 1, 1, 1])` to specify how
results in these directions should be added up. The first rule would be s-wave,
the second extended s-wave.
These rules will be applied for DQMCMeasurements during the simulation. I.e.
first, the normal iteration and summation from `EachLocalQuadByDistance` is
performed. After that we have 5 values for each direction. These are then
weighted by each "symmetry" in `ApplySymmetries` to give the final result, saved
in the DQMCMeasurement.
"""
struct ApplySymmetries{LI <: DeferredLatticeIterator, N, T} <: LatticeIterationWrapper{LI}
    symmetries::NTuple{N, Vector{T}}
end
function ApplySymmetries{LI}(symmetries::Vector{T}...) where {LI <: DeferredLatticeIterator, T}
    ApplySymmetries{LI, length(symmetries), T}(symmetries)
end

struct _ApplySymmetries{LI <: DeferredLatticeIterator, N, T} <: LatticeIterationWrapper{LI}
    iter::LI
    symmetries::NTuple{N, T}
end
function (x::ApplySymmetries{LI})(mc, model) where {LI}
    iter = LI(mc, model)
    _ApplySymmetries(iter, x.symmetries)
end

Base.iterate(s::_ApplySymmetries) = iterate(s.iter)
Base.iterate(s::_ApplySymmetries, state) = iterate(s.iter, state)
Base.length(s::_ApplySymmetries) = length(s.iter)
Base.eltype(s::_ApplySymmetries) = eltype(s.iter)

################################################################################
### Symmetry Wrapper
################################################################################

#=
CCS[dr, (0, x, y, -x, -y)]
Λxx = CCS[:, 2]
Λxxq^L = dot(CCS[:, 2], exp(-dirs .* (1/L, 0))
Λxxq^T = dot(CCS[:, 2], exp(-dirs .* (0, 1/L))
ρs = Λxxq^L - Λxxq^T

# Save:
directions(iter)
# Supply/Pick:
dir_idx (xx)
L

So is longitudinal/transversal just
    1/L * normalize(NN_vector)
    1/L * normalize(⟂ NN_vector)
or is it a reciprocal lattice vector?
=#

"""
    SuperfluidDensity{lattice_iterator_type}(directions, longitudinal, transversal)

`SuperfluidDensity` works similarly to `ApplySymmetries` - it contains some
additional information that should be used to transform the result of of the
given lattice iterator. In this case the additional information is a set of
directional indices `directions` and the related `longitudinal` and `transversal`
vectors in reciprocal space.

The point of this wrapper is to calculate the Fourier transform
`O(q) = ∑_r0 ∑_Δr O(r_0, Δr) (exp(-i q_l Δr) - exp(-i q_t Δr))`
where q_l is `longitudinal` and q_t is transversal. The Δr ised in this formula
are restricted to the given `directions`.
"""
struct SuperfluidDensity{LI <: DeferredLatticeIterator, N} <: LatticeIterationWrapper{LI}
    dir_idxs::NTuple{N, Int}
    long_qs::NTuple{N, Vector{Float64}}
    trans_qs::NTuple{N, Vector{Float64}}
end
function SuperfluidDensity{LI}(directions, longitudinal, transversal) where {LI <: DeferredLatticeIterator}
    SuperfluidDensity{LI, length(directions)}(
        tuple(directions...), tuple(longitudinal...), tuple(transversal...)
    )
end

struct _SuperfluidDensity{LI <: DeferredLatticeIterator, N} <: LatticeIterationWrapper{LI}
    iter::LI

    dirs::Vector{Vector{Float64}}
    dir_idxs::NTuple{N, Int}
    long_qs::NTuple{N, Vector{Float64}}
    trans_qs::NTuple{N, Vector{Float64}}
end
function (x::SuperfluidDensity{LI})(mc, model) where {LI}
    iter = LI(mc, model)
    dirs = directions(lattice(model))
    _SuperfluidDensity(iter, dirs, x.dir_idxs, x.long_qs, x.trans_qs)
end

Base.iterate(s::_SuperfluidDensity) = iterate(s.iter)
Base.iterate(s::_SuperfluidDensity, state) = iterate(s.iter, state)
Base.length(s::_SuperfluidDensity) = length(s.iter)
Base.eltype(s::_SuperfluidDensity) = eltype(s.iter)

ext: jl


text:
defaultdataregistry() = Datasets.FASTAI_DATA_REGISTRY
defaulttaskregistry() = FASTAI_METHOD_REGISTRY

Datasets.listdatasources() = listdatasources(defaultdataregistry())
Datasets.finddatasets(; kwargs...) = finddatasets(defaultdataregistry(); kwargs...)
Datasets.loaddataset(args...; kwargs...) = loaddataset(defaultdataregistry(), args...; kwargs...)

findlearningtasks(args...; kwargs...) = findlearningtasks(defaulttaskregistry(), args...; kwargs...)


text:
using Roots

println(find_zero(x -> x^3 - 3x^2 + 2x, (-100, 100)))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie, Random
CairoMakie.activate!() # HIDE
let
    Random.seed!(124)
    n = 30
    x, y, z = randn(n), randn(n), randn(n)
    fig = Figure(resolution = (800, 450))
    ax = Axis(fig,xgridstyle=:dash,ygridstyle=:dash, xtickalign=1, ytickalign=1)
    pts1= scatter!(ax, x, y, color = z, colormap=:plasma, markersize = 15z,
                marker = :rect)
    pts2= scatter!(ax, 3x, 5y, color = z, colormap= (:thermal, 0.5),
                markersize = 45z, marker = :circle)

    leg1 = Colorbar(fig, pts1, label = "z1 value", ticklabelsize = 14,
                labelpadding = 5, width = 10)
    leg2= Colorbar(fig, pts2, label = "z2 value", ticklabelsize = 14,
                labelpadding = 5, width = 10)
    fig[1, 1] = ax
    fig[1, 2] = leg1
    fig[1, 3] = leg2
    fig
    #save("FilledLine.png"), fig, px_per_unit = 2.0) # HIDE
    save(joinpath(@__DIR__, "output", "BubblePlot.png"), fig, px_per_unit = 2.0) # HIDE
end

using Pkg # HIDE
Pkg.status(["CairoMakie","Random"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
# ---
# title: 687. Longest Univalue Path
# id: problem687
# author: Indigo
# date: 2021-06-29
# difficulty: Medium
# categories: Tree, Recursion
# link: <https://leetcode.com/problems/longest-univalue-path/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, return _the length of the longest path,
# where each node in the path has the same value_. This path may or may not pass
# through the root.
#
# **The length of the path** between two nodes is represented by the number of
# edges between them.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/13/ex1.jpg)
#
#
#
#     Input: root = [5,4,5,1,1,5]
#     Output: 2
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/13/ex2.jpg)
#
#
#
#     Input: root = [1,4,5,4,4,5]
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[0, 104]`.
#   * `-1000 <= Node.val <= 1000`
#   * The depth of the tree will not exceed `1000`.
#
#
## @lc code=start
using LeetCode

function longest_univalue_path(root::TreeNode{Int})
    res = 0
    function lup(root::TreeNode{Int})
        r2 = 1
        rl = lup(root.left)
        rr = lup(root.right)
        r = max(rl, rr) + 1
        if rl != 0 && root.val == root.left.val
            r2 += rl
        end
        if rr != 0 && root.val == root.right.val
            r2 += rr
        end
        res = max(r, res)
        return r
    end
    lup(::Nothing) = 0
    lup(root)
    return res - 1
end
## @lc code=end


text:
module Discord

using Base.Iterators: Pairs

using Dates: DateTime, ISODateTimeFormat, Millisecond, UTC, now, unix2datetime, year

using HTTP: HTTP, Form, Response, StatusError, escapeuri, header, request
using JSON3: JSON3, StructTypes
using Parameters: @with_kw
using SuperEnum: @se

const API_BASE = "https://discord.com/api"
const API_VERSION = 6
const USER_AGENT = let
    toml = read(joinpath(@__DIR__, "..", "Project.toml"), String)
    version = VersionNumber(match(r"version = \"(.*)\"", toml)[1])
    "Julia $VERSION / Discord.jl $version"
end

include("snowflake.jl")
include("objects.jl")
include("rate_limiter.jl")
include("clients.jl")
include("routes.jl")

end


text:
"""
    lenient_min_max_complex(A::AbstractArray)

Returns (min, max) of hypot.(A)
"""
function lenient_min_max_complex(A::AbstractArray)
    magnitude = hypot.(A)
    extrema(filter(x -> !isnan(x) && !isinf(x), magnitude))
end

"""
    lenient_min_max(A)

Neglecting NaN and Inf values, return
- Minimum and maximum value out of real element arrays.
- Minimum and maximum magnitude out of complex element arrays.
- Minimum and maximum magnitude out of tuple element arrays.
"""
lenient_min_max(A::AbstractArray{<:RealQuantity}) = extrema(filter(x-> !isnan(x) && !isinf(x), A))
lenient_min_max(A::AbstractArray{<:Real}) = extrema(filter(x-> !isnan(x) && !isinf(x), A))
lenient_min_max(A) = lenient_min_max_complex(A)

""""
    lenient_min_max(f_xy, xs, ys)

    f_xy is a function with domain defined by iterators (xs, ys)
    Neglecting NaN and Inf values, return
- Minimum and maximum value out of real valued function f_xy
- Minimum and maximum magnitude out of complex valued function f_xy
- Minimum and maximum magnitude out of tuple valued functions
"""
function lenient_min_max(f_xy, xs, ys)
    bigiterator = product(xs, ys)
    default = hypot(f_xy(first(bigiterator)...)) * 0.0
    function g(x, y)
        valu = hypot(f_xy(x, y))
        !isnan(valu) && !isinf(valu) ? valu : default
    end
    extrema(tu -> g(tu...), bigiterator)
end

"""
    lenient_min(A)
    lenient_min(f_xy, xs, ys)

Neglecting NaN and Inf values, return
- Minimum value out of real element arrays.
- Minimum magnitude out of complex element arrays.
- Minimum magnitude out of tuple element arrays.
- Minimum value out of real valued function f_xy on all of xs, ys
- Minimum magnitude out of complex valued function f_xy
- Minimum magnitude out of tuple valued functions
"""
lenient_min(A) = lenient_min_max(A)[1]
lenient_min(f_xy, xs, ys) = lenient_min_max(f_xy, xs, ys)[1]

"""
    lenient_max(A)
    lenient_max(f_xy, xs, ys)

Neglecting NaN and Inf values, return
- Maximum value out of real element arrays.
- Maximum magnitude out of complex element arrays.
- Maximum magnitude out of tuple element arrays.
- Maximum value out of real valued function f_xy
- Maximum magnitude out of complex valued function f_xy
- Maximum magnitude out of tuple valued functions
"""
lenient_max(A) = lenient_min_max(A)[2]
lenient_max(f_xy, xs, ys) = lenient_min_max(f_xy, xs, ys)[2]

magnitude(x::Complex) = hypot(x)
magnitude(x::ComplexQuantity) = hypot(x)
magnitude(x::NTuple{2, T}) where T = hypot(x[1], x[2])

qangle(x::Complex) = °(angle(x))
qangle(x::ComplexQuantity) = °(angle(x))
qangle(x::NTuple{2, T}) where T = °(atan(x[1], x[2]))


text:
function extract_all_lnns(x)
    lnns = Vector{LineNumberNode}(undef, 0)
    extract_all_lnns!(lnns, x)
    return lnns
end

function extract_all_lnns!(lnns::AbstractVector{<:LineNumberNode}, lnn::LineNumberNode)
    push!(lnns, lnn)
    return lnns
end
function extract_all_lnns!(lnns::AbstractVector{<:LineNumberNode}, ex::Expr)
    for arg in ex.args
        extract_all_lnns!(lnns::AbstractVector{<:LineNumberNode}, arg)
    end
    return lnns
end
function extract_all_lnns!(lnns::AbstractVector{<:LineNumberNode}, ::Any)
    return lnns
end

function prepend_lnns!(ex::Expr, lnns::AbstractVector{<:LineNumberNode})
    return prepend_lnns!(ex, lnns, Val(ex.head))
end
function prepend_lnns!(ex::Expr, lnns::AbstractVector{<:LineNumberNode}, ::Val{:block})
    for lnn in lnns
        pushfirst!(ex.args, Expr(:block, lnn, :(nothing)))
    end
    return ex
end


text:
abstract type NormalParams <: SolverParams end

abstract type PreallocatedDataNormal{T <: Real, S} <: PreallocatedData{T, S} end

abstract type PreallocatedDataNormalChol{T <: Real, S} <: PreallocatedDataNormal{T, S} end

include("K1CholDense.jl")

abstract type PreallocatedDataNormalKrylov{T <: Real, S} <: PreallocatedDataNormal{T, S} end

include("K1Krylov.jl")

abstract type PreallocatedDataNormalStructured{T <: Real, S} <: PreallocatedDataNormal{T, S} end

include("K1_1Structured.jl")
include("K1_2Structured.jl")


text:
module KPP

using OceanTurb

import OceanTurb: Constants

using Base: @propagate_inbounds

const nsol = 4

@solution U V T S

"""
    Parameters(; kwargs...)

Construct KPP parameters.
"""
@Base.kwdef struct Parameters{T<:AbstractFloat} <: AbstractParameters
    CSL   :: T  = 0.1   # Surface layer fraction
    Cτ    :: T  = 0.4   # Von Karman constant
    CNL   :: T  = 6.33  # Non-local flux proportionality constant

    Cstab :: T  = 2.0   # Stable buoyancy flux parameter for wind-driven turbulence
    Cunst :: T  = 6.4   # Unstable buoyancy flux parameter for wind-driven turbulence

       Cn :: T  = 1.0   # Exponent for effect of stable buoyancy forcing on wind mixing
    Cmτ_U :: T  = 0.25  # Exponent for effect of unstable buoyancy forcing on wind mixing of U
    Cmτ_T :: T  = 0.5   # Exponent for effect of unstable buoyancy forcing on wind mixing of T
    Cmb_U :: T  = 1/3   # Exponent for the effect of wind on convective mixing of U
    Cmb_T :: T  = 1/3   # Exponent for effect of wind on convective mixing of T

    Cd_U  :: T  = 0.5   # Wind mixing regime threshold for momentum
    Cd_T  :: T  = 2.5   # Wind mixing regime threshold for tracers

    Cb_U  :: T  = 0.599 # Buoyancy flux parameter for convective turbulence
    Cb_T  :: T  = 1.36  # Buoyancy flux parameter for convective turbulence
    Cτb_U :: T  = (Cτ / Cb_U)^(1/Cmb_U) * (1 + Cunst*Cd_U)^(Cmτ_U/Cmb_U) - Cd_U  # Wind stress parameter for convective turbulence
    Cτb_T :: T  = (Cτ / Cb_T)^(1/Cmb_T) * (1 + Cunst*Cd_T)^(Cmτ_T/Cmb_T) - Cd_T  # Wind stress parameter for convective turbulence

    CRi   :: T  = 0.3   # Critical bulk Richardson number
    CKE   :: T  = 4.32  # Unresolved turbulence parameter
    CKE2   :: T  = 0.0  # Unresolved turbulence parameter 2
    CKE3   :: T  = 0.0  # Unresolved turbulence parameter 3
    CKE4 :: T = 0.0 #Unresolved turbulence parameter 4
    CKE₀  :: T  = 1e-11 # Minimum unresolved turbulence kinetic energy

    KU₀   :: T  = 1e-6  # Interior viscosity for velocity
    KT₀   :: T  = 1e-7  # Interior diffusivity for temperature
    KS₀   :: T  = 1e-9  # Interior diffusivity for salinity
end

# Shape functions.
# 'd' is a non-dimensional depth coordinate.
default_NL_shape(d) = ifelse(0<d<1, d*(1-d)^2, -zero(d))
const default_K_shape = default_NL_shape

mutable struct State{T} <: FieldVector{6, T}
    Qu :: T
    Qv :: T
    Qθ :: T
    Qs :: T
    Qb :: T
    h  :: T
end

State(T=Float64) = State{T}(0, 0, 0, 0, 0, 0)

"""
    Forcing(; U=addzero, V=addzero, T=addzero, S=addzero)

Construct a `NamedTuple` of forcing functions for KPP `Model`s for each
field `U, V, T, S`. The functions must have the signature `forcing(model::Model, i)`,
where `i` is the vertical index at which the forcing is applied.
"""
Forcing(; U=addzero, V=addzero, T=addzero, S=addzero) = (U=U, V=V, T=T, S=S)

"""
    ModelBoundaryConditions([FT=Float64;] U = DefaultBoundaryConditions(FT),
                                          V = DefaultBoundaryConditions(FT),
                                          T = DefaultBoundaryConditions(FT),
                                          S = DefaultBoundaryConditions(FT))

Returns a `NamedTuple` of boundary conditions for a `KPP.Model` with solution
fields `U`, `V`, `T`, `S`.

Example
=======

julia> surface_temperature_flux(model) = cos(model.clock.time)

julia> T_bcs = BoundaryConditions(top = FluxBoundaryCondition(surface_flux))

julis> bcs = KPP.ModelBoundaryConditions(T=T_bcs)
"""
function ModelBoundaryConditions(FT=Float64; U = DefaultBoundaryConditions(FT),
                                             V = DefaultBoundaryConditions(FT),
                                             T = DefaultBoundaryConditions(FT),
                                             S = DefaultBoundaryConditions(FT))
    return (U=U, V=V, T=T, S=S)
end

addzero(args...) = 0

"""
    Model{S, G, T, U, B, F} <: AbstractModel{S, G, T}

Struct for KPP models.
"""
mutable struct Model{S, G, T, U, B, F} <: AbstractModel{S, G, T}
    clock       :: Clock{T}
    grid        :: G
    timestepper :: S
    solution    :: U
    bcs         :: B
    parameters  :: Parameters{T}
    constants   :: Constants{T}
    state       :: State{T}
    forcing     :: F
end

"""
    Model(; kwargs...)

Construct a KPP Model.
"""
function Model(; N=10, L=1.0,
            grid = UniformGrid(N, L),
       constants = Constants(),
      parameters = Parameters(),
         stepper = :ForwardEuler,
             bcs = ModelBoundaryConditions(eltype(grid)),
         forcing = Forcing()
    )

     K = (U=KU, V=KV, T=KT, S=KS)
     R = (U=RU, V=RV, T=RT, S=RS)
    eq = Equation(R=R, K=K, update=update_state!)

    solution = Solution(
        CellField(grid),
        CellField(grid),
        CellField(grid),
        CellField(grid)
    )

    lhs = OceanTurb.build_lhs(solution)
    timestepper = Timestepper(stepper, eq, solution, lhs)
    clock = Clock()
    state = State()

    return Model(clock, grid, timestepper, solution, bcs, parameters, constants, state,
                 forcing)
end

# Note: we use 'm' to refer to 'model' in function definitions below.

@inline Qb(g, α, β, Qθ, Qs) = g * (α*Qθ - β*Qs)

@propagate_inbounds d(m, i) = ifelse(m.state.h>0, -m.grid.zf[i]/m.state.h, -zero(m.state.h))

"Return the buoyancy gradient at face point i."
@propagate_inbounds ∂B∂z(T, S, g, α, β, i) = g * (α*∂z(T, i) - β*∂z(S, i))

"""
    update_state!(model)

Update the top flux conditions and mixing depth for `model`
and store in `model.state`.
"""
function update_state!(m)
    m.state.Qu = getbc(m, m.bcs.U.top)
    m.state.Qv = getbc(m, m.bcs.V.top)
    m.state.Qθ = getbc(m, m.bcs.T.top)
    m.state.Qs = getbc(m, m.bcs.S.top)
    m.state.Qb = Qb(m.constants.g, m.constants.α, m.constants.β, m.state.Qθ, m.state.Qs)
    m.state.h  = mixing_depth(m)
    return nothing
end

#
# Diagnosis of mixing depth "h"
#

"Returns the surface_layer_average for mixing depth h = -zf[i]."
@propagate_inbounds function surface_layer_average(c, CSL::T, i) where T
    if i > c.grid.N # Return surface value
        return onface(c, c.grid.N+1)
    else
        iε = length(c)+1 - CSL*(length(c)+1 - i) # (fractional) face "index" of the surface layer
        face = ceil(Int, iε)  # next cell face above the fractional depth
        frac = face - iε # fraction of lowermost cell in the surface layer.
        surface_layer_integral = zero(T)

        # Contribution of fractional cell to total integral
        if face > 1
            surface_layer_integral += frac * Δf(c, face-1) * c[face-1]
        else
            face = 1
        end

        # Add cells above face, if there are any.
        for j = face:length(c)
            surface_layer_integral += Δf(c, j) * c[j]
        end

        h = -c.grid.zf[i] # depth

        return surface_layer_integral / (CSL*h)
    end
end

"""
Return Δc(hᵢ), the difference between the surface-layer average of c and its value at depth hᵢ, where
i is a face index.
"""
@propagate_inbounds Δ(c, CSL, i) = surface_layer_average(c, CSL, i) - onface(c, i)

"Returns the parameterization for unresolved KE at face point i."
@inline function unresolved_kinetic_energy(h, Bz, Qb, CKE, CKE₀, g, α, β, CKE2, CKE3, CKE4)
    N² = max(0, Bz)
    w1 = max(0, h * Qb)^(1/3)
    w2 = sqrt(N² * h^2)
    tmp = CKE * w1 * w2 + CKE₀
    tmp += CKE2 * (CKE3 * h)^(1+CKE4) * w1^(1-CKE4) #CKE3 has to have units of N
    # tmp += CKE2 * w1^CKE3 * (w2)^(2-CKE3)
    return tmp
end

@inline function unresolved_kinetic_energy(h, Bz, Qb, CKE, CKE₀, g, α, β)
    tmp = CKE * h^(4/3) * sqrt(max(0, Bz)) * max(0, Qb)^(1/3) + CKE₀
    return tmp
end

"""
    bulk_richardson_number(model, i)

Returns the bulk Richardson number of `model` at face `i`.
"""
@propagate_inbounds function bulk_richardson_number(
            U, V, T, S, Qb::TT, CKE::TT, CKE₀::TT, CSL::TT,
            g::TT, α::TT, β::TT, CKE2::TT, CKE3::TT, CKE4::TT, i) where TT

    h = -U.grid.zf[i]
    # (h - hε) * ΔB
    # h⁺ΔB = h * (one(TT) - CSL/2) * g * (α*Δ(T, CSL, i) - β*Δ(S, CSL, i))
    # get rid of the CSL in the numerator
    h⁺ΔB = h * g * (α*Δ(T, CSL, i) - β*Δ(S, CSL, i))

    KE = (Δ(U, CSL, i)^2 + Δ(V, CSL, i)^2
              + unresolved_kinetic_energy(h, ∂B∂z(T, S, g, α, β, i), Qb, CKE, CKE₀, g, α, β, CKE2, CKE3, CKE4))

    if KE == 0 && h⁺ΔB == 0 # Alistar Adcroft's theorem
        return -zero(TT)
    else
        return h⁺ΔB / KE
    end
end

@propagate_inbounds bulk_richardson_number(m, i) = bulk_richardson_number(
    m.solution.U, m.solution.V, m.solution.T, m.solution.S,
    m.state.Qb, m.parameters.CKE, m.parameters.CKE₀, m.parameters.CSL, m.constants.g,
    m.constants.α, m.constants.β,m.parameters.CKE2,m.parameters.CKE3, m.parameters.CKE4, i)

"""
    mixing_depth(model)

Calculate the mixing depth 'h' for `model`.
"""
function mixing_depth(m)
    ih₁ = m.grid.N + 1 # start at top.
    @inbounds Ri₁ = bulk_richardson_number(m, ih₁) # should be 0.

    # Descend through grid until Ri rises above critical value
    while ih₁ > 1 && Ri₁ < m.parameters.CRi
        ih₁ -= 1 # descend
        @inbounds Ri₁ = bulk_richardson_number(m, ih₁)
    end

    # Edge cases:
    # 1. Mixing depth is at the top of the domain (z=0):
    if ih₁ == m.grid.N + 1
        @inbounds z★ = m.grid.zf[ih₁]

    # 2. Mixing depth is whole domain because Ri is always less than CRi:
    elseif ih₁ == 1 && Ri₁ < m.parameters.CRi
        @inbounds z★ = m.grid.zf[ih₁]

    # 3. Ri is infinite somewhere inside the domain.
    elseif !isfinite(Ri₁)
        @inbounds z★ = m.grid.zc[ih₁]

    # Main case: mixing depth is in the interior.
    else # Ri₁ > CRi
        ΔRi = bulk_richardson_number(m, ih₁+1) - Ri₁ # <0 linearly interpolate to find h.
        # x = x₀ + Δx * (y-y₀) / Δy
        @inbounds z★ = m.grid.zf[ih₁] + Δf(m.grid, ih₁) * (m.parameters.CRi - Ri₁) / ΔRi
    end

    -z★ < 0 && @warn "mixing depth $(-z★) is negative"

    return -z★ # "depth" is negative height.
end

#
# Vertical velocity scale
#

"Return true if the boundary layer is unstable and convecting."
@inline isunstable(model) = model.state.Qb > 0

"Return true if the boundary layer is forced."
@inline isforced(model) = model.state.Qu != 0 || model.state.Qv != 0 || model.state.Qb != 0

"Return the turbuent velocity scale associated with wind stress."
@inline ωτ(Qu, Qv) = (Qu^2 + Qv^2)^(1/4)
@inline ωτ(m::AbstractModel) = ωτ(m.state.Qu, m.state.Qv)

"Return the turbuent velocity scale associated with convection."
@inline ωb(Qb, h) = abs(h * Qb)^(1/3)
@inline ωb(m::AbstractModel) = ωb(m.state.Qb, m.state.h)

"Return the vertical velocity scale at depth d for a stable boundary layer."
@inline 𝒲_stable(Cτ, Cstab, Cn, ωτ, ωb, d) = Cτ * ωτ / (1 + Cstab * d * (ωb/ωτ)^3)^Cn

"Return the vertical velocity scale at scaled depth dϵ for an unstable boundary layer."
@inline function 𝒲_unstable(CSL, Cd, Cτ, Cunst, Cb, Cτb, Cmτ, Cmb, ωτ, ωb, d)
    dϵ = min(CSL, d)
    if dϵ * ωb^3 < Cd * ωτ^3
        return Cτ * ωτ * (1 + Cunst * dϵ * (ωb/ωτ)^3)^Cmτ
    else
        return Cb * ωb * (dϵ + Cτb * (ωτ/ωb)^3)^Cmb
    end
end

@propagate_inbounds function 𝒲_unstable_U(m, i)
    return 𝒲_unstable(m.parameters.CSL, m.parameters.Cd_U,
                            m.parameters.Cτ, m.parameters.Cunst,
                            m.parameters.Cb_U, m.parameters.Cτb_U,
                            m.parameters.Cmτ_U, m.parameters.Cmb_U,
                            ωτ(m), ωb(m), d(m, i)
                            )
end

@propagate_inbounds function 𝒲_unstable_T(m, i)
    return 𝒲_unstable(m.parameters.CSL, m.parameters.Cd_T,
                            m.parameters.Cτ, m.parameters.Cunst,
                            m.parameters.Cb_T, m.parameters.Cτb_T,
                            m.parameters.Cmτ_T, m.parameters.Cmb_T,
                            ωτ(m), ωb(m), d(m, i)
                            )
end

@propagate_inbounds function 𝒲_stable(m, i)
    return 𝒲_stable(m.parameters.Cτ, m.parameters.Cstab, m.parameters.Cn,
                          ωτ(m), ωb(m), d(m, i)
                          )
end

"Return the turbulent velocity scale for momentum at face point i."
@propagate_inbounds function 𝒲_U(m::AbstractModel{TS, G, T}, i) where {TS, G, T}
    if !isforced(m)
        return -zero(T)
    elseif isunstable(m)
        return 𝒲_unstable_U(m, i)
    else
        return 𝒲_stable(m, i)
    end
end

"Return the turbulent velocity scale for tracers at face point i."
@propagate_inbounds function 𝒲_T(m::AbstractModel{TS, G, T}, i) where {TS, G, T}
    if !isforced(m)
        return -zero(T)
    elseif isunstable(m)
        return 𝒲_unstable_T(m, i)
    else
        return 𝒲_stable(m, i)
    end
end

const 𝒲_V = 𝒲_U
const 𝒲_S = 𝒲_T

## ** The K-Profile-Parameterization **
K_KPP(h, 𝒲, d, shape=default_K_shape) = ifelse(0<d<1, max(zero(h), h*𝒲*shape(d)), -zero(h))

#
# Non-local flux
#

"""
    NL(CNL, flux, d, shape=default_shape)

Returns the nonlocal flux, N = CNL*flux*shape(d),
where `flux` is the flux of some quantity out of the surface,
`shape` is a shape function, and `d` is a non-dimensional depth coordinate
that increases from 0 at the surface to 1 at the bottom of the mixing layer.

Because flux is defined as pointing in the positive direction,
a positive surface flux implies negative surface flux divergence,
which implies a reduction to the quantity in question.
For example, positive heat flux out of the surface implies cooling.
"""
@inline NL(CNL, flux, d, shape=default_NL_shape) = CNL * flux * shape(d)

@inline function ∂NL∂z(CNL::T, Qϕ, dᵢ₊₁, dᵢ, Δf, m) where T
    if isunstable(m)
        return (NL(CNL, Qϕ, dᵢ₊₁) - NL(CNL, Qϕ, dᵢ)) / Δf
    else
        return -zero(T)
    end
end

@propagate_inbounds ∂NLT∂z(m, i) =
    ∂NL∂z(m.parameters.CNL, m.state.Qθ, d(m, i+1), d(m, i), Δf(m.grid, i), m)

@propagate_inbounds ∂NLS∂z(m, i) =
    ∂NL∂z(m.parameters.CNL, m.state.Qs, d(m, i+1), d(m, i), Δf(m.grid, i), m)

#
# Equation specification
#

# K_{U,V,T,S} is calculated at face points
@propagate_inbounds KU(m, i) = K_KPP(m.state.h, 𝒲_U(m, i), d(m, i)) + m.parameters.KU₀
@propagate_inbounds KT(m, i) = K_KPP(m.state.h, 𝒲_T(m, i), d(m, i)) + m.parameters.KT₀
@propagate_inbounds KS(m, i) = K_KPP(m.state.h, 𝒲_S(m, i), d(m, i)) + m.parameters.KS₀
const KV = KU

@propagate_inbounds RU(m, i) =   m.constants.f * m.solution.V[i] + m.forcing.U(m, i)
@propagate_inbounds RV(m, i) = - m.constants.f * m.solution.U[i] + m.forcing.V(m, i)
@propagate_inbounds RT(m, i) = - ∂NLT∂z(m, i) + m.forcing.T(m, i)
@propagate_inbounds RS(m, i) = - ∂NLS∂z(m, i) + m.forcing.S(m, i)

#####
##### Some utilities
#####

function nonlocal_salinity_flux!(flux, m)
    for i in interiorindices(flux)
        @inbounds flux[i] = NL(m.parameters.CNL, m.state.Qs, d(m, i))
    end
    return nothing
end

function nonlocal_temperature_flux!(flux, m)
    for i in interiorindices(flux)
        @inbounds flux[i] = NL(m.parameters.CNL, m.state.Qθ, d(m, i))
    end
    return nothing
end

function nonlocal_salinity_flux(model)
    flux = FaceField(model.grid)
    nonlocal_salinity_flux!(flux, model)
    return flux
end

function nonlocal_temperature_flux(model)
    flux = FaceField(model.grid)
    nonlocal_temperature_flux!(flux, model)
    return flux
end

end # module


text:
# copies

export append_copy!, append_fill!, append_prefetch!, append_advise!

append_copy!(list::ZeCommandList, dst::Union{Ptr,ZePtr}, src::Union{Ptr,ZePtr},
             size::Integer, signal_event::Union{ZeEvent,Nothing}=nothing,
             wait_events::ZeEvent...) =
    zeCommandListAppendMemoryCopy(list, dst, src, size, something(signal_event, C_NULL),
                                   length(wait_events), [wait_events...])

append_fill!(list::ZeCommandList, ptr::Union{Ptr,ZePtr}, pattern::Union{Ptr,ZePtr},
             pattern_size::Integer, size::Integer,
             signal_event::Union{ZeEvent,Nothing}=nothing, wait_events::ZeEvent...) =
    zeCommandListAppendMemoryFill(list, ptr, pattern, pattern_size, size,
                                  something(signal_event, C_NULL), length(wait_events),
                                  [wait_events...])

append_prefetch!(list::ZeCommandList, ptr::Union{Ptr,ZePtr}, size::Integer) =
    zeCommandListAppendMemoryPrefetch(list, ptr, size)

append_advise!(list::ZeCommandList, dev::ZeDevice, ptr::Union{Ptr,ZePtr}, size::Integer,
               advise::ze_memory_advice_t) =
    zeCommandListAppendMemAdvise(list, dev, ptr, size, advise)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "opusfile"
version = v"0.12.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/xiph/opusfile/releases/download/v$(version.major).$(version.minor)/opusfile-$(version.major).$(version.minor).tar.gz",
                  "118d8601c12dd6a44f52423e68ca9083cc9f2bfe72da7a8c1acb22a80ae3550b"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/opusfile-*

if [[ "${target}" == *mingw* ]]; then
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/no-openssl-wincert.patch"
fi

./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --with-pic
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(!Sys.isfreebsd, supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libopusurl", :libopusurl),
    LibraryProduct("libopusfile", :libopusfile)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Ogg_jll", uuid="e7412a2a-1a6e-54c0-be00-318e2571c051")),
    Dependency(PackageSpec(name="Opus_jll", uuid="91d4177d-7536-5919-b921-800302f37372")),
    Dependency(PackageSpec(name="OpenSSL_jll", uuid="458c3c95-2e84-50aa-8efc-19380b2a3a95")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
# # Radial Basis Function Surrogate Models

#src We generate a documentation page from this source file.
#src Let's tell the user:
#md # This file is automatically generated from source code.
#md # For usage examples refer to [Summary & Quick Examples](@ref rbf_summary).

# ## Intro and Prerequisites

# We want to offer radial basis function (RBF) surrogate models (implementing the `SurrogateModel` interface).
# To this end, we leverage the package [`RadialBasisFunctionModels.jl`](https://github.com/manuelbb-upb/RadialBasisFunctionModels.jl).
# A scalar RBF model consists of a ``n``-variate Polynomial and linear combination of shifted radial kernels.
# For more information, see [the documentation of `RadialBasisFunctionModels.jl`](https://manuelbb-upb.github.io/RadialBasisFunctionModels.jl/stable).

import RadialBasisFunctionModels as RBF
using LinearAlgebra: qr, Hermitian, cholesky, inv, I, givens, diag

# The polyonmials will have a degree of at most 1.
# To construct “good” linear polynomials, we need to make sure to have construction sites
# that span the decision space well. Such a set of construction sites is called Λ-poised or
# sufficiently affinely independent.
# The file `AffinelyIndependentPoints.jl` implements some helpers to find suitable points
# as described by Wild et. al.[^wild_diss]

include("AffinelyIndependentPoints.jl")

# ## Surrogate Interface Implementations

# The model used in our algorithm simply wraps an interpolation model from the `RBF` package.
@with_kw struct RbfModel{R} <: SurrogateModel
	model :: R

	## indicator: is the model fully linear?
	fully_linear :: Bool = false
end

fully_linear( rbf :: RbfModel ) :: Bool = rbf.fully_linear

# We offer a large range of configuration parameters in the `RBFConfig`, which implements
# a `SurrogateConfig`.
"""
    RbfConfig(; kwarg1 = val1, … )

Configuration type for local RBF surrogate models.

$(FIELDS)

"""
@with_kw mutable struct RbfConfig <: SurrogateConfig
	"(default `:cubic`) RBF kernel (Symbol), either `:cubic`, `:inv_multiquadric`, `:multiquadric`, `:exp` or `:thin_plate_spline`."
	kernel :: Symbol = :cubic

	"(default `NaN` for automatic) RBF shape paremeter, either a number or a string containing `Δ`."
	shape_parameter :: Union{String, Float64} = NaN

	"(default `1`) Degree of polynomial attached to RBF. `-1` means no polynomial."
	polynomial_degree :: Int64 = 1;

	"(default `2`) Local enlargment factor of trust region for sampling."
	θ_enlarge_1 :: Float64 = 2

	"(default `2`) Maximum enlargment factor of maximum trust region for sampling."
	θ_enlarge_2 :: Float64 = 2

	"(default `1/(2*θ_enlarge_1)` Sampling parameter to generate Λ-poised set. The higher, the more poised."
	θ_pivot :: Float64 = 1 / (2 * θ_enlarge_1)

	"(default `1e-7`) Parameter for 2nd sampling algorithm to ensure boundedness of Cholesky factors."
	θ_pivot_cholesky :: Float64 = 1e-7

	"(default `false`) Require models to be fully linear in each iteration."
	require_linear :: Bool = false

	"(default `-1`) Maximum number of training sites. `-1` is reset to `2n+1`."
	max_model_points :: Int64 = -1 # is probably reset in the algorithm
	"(default `false`) Sample new sites to always use the maximum number of points."
	use_max_points :: Bool = false

	"Whether or not to re-construct the training set in each iteration."
	optimized_sampling = true

	"(default `typemax(Int64)`) Maximum number of objective evaluations."
	max_evals :: Int64 = typemax(Int64)

	@assert θ_enlarge_1 * θ_pivot ≤ 1 "θ_pivot must be <= θ_enlarge_1^(-1)."

	##	@assert sampling_algorithm ∈ [:orthogonal, :monte_carlo] "Sampling algorithm must be either `:orthogonal` or `:monte_carlo`."
	@assert kernel ∈ Symbol.(["gaussian", "inv_multiquadric", "multiquadric", "cubic", "thin_plate_spline"]) "Kernel '$kernel' not supported yet."

	## Some sanity checks for the shape parameters:
	@assert kernel != :thin_plate_spline || ( isnan(shape_parameter) || shape_parameter % 1 == 0 && shape_parameter >= 1 ) "Invalid shape_parameter for :thin_plate_spline."
	@assert kernel != :cubic || ( isnan(shape_parameter) || shape_parameter % 1 == 0 && shape_parameter % 2 == 1 ) "Invalid shape_parameter for :cubic."
	@assert (isa( shape_parameter, String ) || isnan(shape_parameter)) || shape_parameter > 0 "Shape parameter must be strictly positive."
	@assert θ_enlarge_1 >=1 && θ_enlarge_2 >=1 "θ's must be >= 1."
end

# The required method implementations are straightforward.
# Note, thate we allow the models to be combined to vector functions if they
# share the same configuration to avoid redundant efforts whilst constructing models.

max_evals( cfg :: RbfConfig ) :: Int = cfg.max_evals
combinable( cfg :: RbfConfig ) :: Bool = true
combine(cfg1 :: RbfConfig, :: RbfConfig) :: RbfConfig = cfg1

# We also need to introduce our own implementation for `isequal` and `hash` for
# `RbfConfig`s to be combinable, see [the docs too](https://docs.julialang.org/en/v1/base/base/).
function Base.hash( cfg :: RbfConfig, h :: UInt )
	return hash( getfield.( cfg, Tuple( fn for fn ∈ fieldnames(RbfConfig) ) ), h )
end
function Base.isequal( cfg1 :: RbfConfig, cfg2 :: RbfConfig )
	all( isequal( getfield(cfg1, fn), getfield(cfg2, fn) ) for fn in fieldnames( RbfConfig) )
end

# To allow the user to set the shape parameter relative to the current trust region radius
# using a verbose string, we need this little helper function, which evaluates the string.

function parse_shape_param_string( Δ :: F, expr_str) :: F where F
    ex = Meta.parse(expr_str)
    sp = @eval begin
        let Δ=$Δ
            $ex
        end
    end
	return sp
end

# The `RbfMeta` is used to store construction and update data for the models.
# To be specific, we have several inidices lists that store database indices
# of (potentially unevaluated) results that are later used for fitting the model.
@with_kw mutable struct RbfMeta{F<:AbstractFloat} <: SurrogateMeta
    center_index :: Int = -1
    round1_indices :: Vector{Int} = []
    round2_indices :: Vector{Int} = []
    round3_indices :: Vector{Int} = []
    round4_indices :: Vector{Int} = []
    fully_linear :: Bool = false
	improving_directions :: Vector{Vector{F}} = []
end

get_saveable_type( meta :: T ) where {T<:RbfMeta} = T
get_saveable( meta :: RbfMeta ) = deepcopy(meta)

# A little helper to retrieve all those indices:
function _collect_indices( meta :: RbfMeta; include_x = true ) :: Vector{Int}
	return [
		include_x ? meta.center_index : [];
		meta.round1_indices;
		meta.round2_indices;
		meta.round3_indices;
		meta.round4_indices
	]
end

# And a helper, to partially copy some data from `src` to `dest`.
# This is due to the fact, that the first 3 rounds of construction
# data gathering are the same for all possible RBF models and we can safe
# some effort.
function copy_meta!(dest, src)
	dest.center_index = src.center_index
	for fn in [ Symbol("round$(i)_indices") for i = 1: 3 ]
		dest_arr = getfield(dest, fn)
		empty!( dest_arr )
		append!( dest_arr, getfield( src, fn) )
	end
	empty!(dest.improving_directions)
	append!(dest.improving_directions, src.improving_directions)
end

export RbfConfig, RbfMeta, RbfModel

# ## Construction

# The initial `prepare_init_model` function should return a meta object that can be used
# to build an initial surrogate model.
# We delegate the work to `prepare_update_model`.
function prepare_init_model( cfg :: RbfConfig, objf :: AbstractObjective, mop :: AbstractMOP,
	id :: AbstractIterData, db :: AbstractDB, ac :: AbstractConfig;
	ensure_fully_linear = true, kwargs...)
	F = eltype( get_x(id) )
	meta = RbfMeta{F}()
	return prepare_update_model(nothing, objf, meta, mop, id, db, ac; ensure_fully_linear = true, kwargs... )
end

# Usually, `prepare_update_model` would only accept a model as its first argument.
# Because of the trick from above, we actually allow `nothing`, too.
function prepare_update_model( mod :: Union{Nothing, RbfModel}, objf :: AbstractObjective, meta :: RbfMeta,
	mop :: AbstractMOP, iter_data :: AbstractIterData, db :: AbstractDB, algo_config :: AbstractConfig;
	ensure_fully_linear = false, force_rebuild = false, meta_array = nothing )

	!force_rebuild && @logmsg loglevel2 "Trying to find results for fitting an RBF model."

	## Retrieve current iteration information and some meta data.
	Δ = get_delta(iter_data)
	Δ_max = get_delta_max(algo_config)
	x = get_x(iter_data)
	x_index = get_x_index(iter_data)
	cfg = model_cfg( objf )

	F = eltype(x)
	n_vars = length(x)

	## Can we skip the first rounds? (Because we already found interpolation sets for other RBFModels?)
	all_objfs = list_of_objectives(mop)
	skip_first_rounds = false
	for (i,other_meta) in enumerate(meta_array)
		other_objf = all_objfs[i]
		if other_meta isa RbfMeta
			other_cfg = model_cfg(other_objf)
			if other_cfg.θ_pivot == cfg.θ_pivot && other_cfg.θ_enlarge_1 == cfg.θ_enlarge_1 &&
				other_cfg.θ_enlarge_2 == cfg.θ_enlarge_2 && other_cfg.optimized_sampling == cfg.optimized_sampling
				copy_meta!( meta, other_meta )
				skip_first_rounds = true
			end
		end
	end

	## By default, assume that our model is not fully linear
	meta.fully_linear = false

	## use center as first training site ⇒ at least `n_vars` required still
	meta.center_index = x_index

	## First round of sampling:
	### Try to find points in slightly enlarged trust region
	Δ_1 = F.(cfg.θ_enlarge_1 * Δ)
	lb_1, ub_1 = local_bounds( mop, x, Δ_1 )
	piv_val_1 = F.(cfg.θ_pivot * Δ_1) # threshold value for acceptance in filter

	### `Δ_2` is the maximum allowed trust region radius and used in rounds 2 & 4
	Δ_2 = F.(cfg.θ_enlarge_2 * Δ_max )
	lb_2, ub_2 = local_bounds( mop, x, Δ_2 )
	piv_val_2 = piv_val_1 # the pivot value stays the same

	skip_first_rounds && @goto round4

	if force_rebuild || !cfg.optimized_sampling
		### `force_rebuild` makes us skip the point searching procedures to …
		### … rebuild the model along the coordinate axes.
		filtered_indices_1 = Int[]
		improving_directions = [ [zeros(F,i-1); one(F); zeros(F,n_vars - i)] for i = 1:n_vars ]
	else
		@logmsg loglevel3 "Round1: Inspect box with radius $(Δ_1) and pivot value $(piv_val_1)."

		### only consider points from within current trust region …
		candidate_indices_1 = results_in_box_indices( db, lb_1, ub_1, [x_index],)

		### … and filter them to obtain affinely independent points.
		filter = AffinelyIndependentPointFilter(;
			x_0 = x,
			seeds = get_site.(db, candidate_indices_1),
			return_indices = true,
			pivot_val = piv_val_1
		)

		filtered_indices_1 = candidate_indices_1[ collect( filter ) ]
		### TODO should we rather use Z₂ to sample along unexplored directions? (for now, i simply reverse Z₁) #src
		improving_directions = reverse(collect(Vector{F}, eachcol(filter.Z)))

		@logmsg loglevel3 "Round1: Found $(length(filtered_indices_1)) sites in database."
	end
	### Store indices in meta data object:
	empty!(meta.round1_indices)
	append!(meta.round1_indices, filtered_indices_1)
	empty!(meta.improving_directions)
	append!(meta.improving_directions, improving_directions )

	## Second round of sampling:
	### If there are not enough sites to have a fully linear model …
	### … try to at least find more sites in maximum allowed radius
	n_missing = n_vars - length( meta.round1_indices )

	if n_missing == 0 || force_rebuild || !cfg.optimized_sampling || ensure_fully_linear || Δ ≈ Δ_max && cfg.θ_enlarge_1 == cfg.θ_enlarge_2
		@logmsg loglevel3 "Skipping round 2."

		meta.fully_linear = true
		filter_2 = filter
		empty!(meta.round2_indices)
	else
		### actually perform round 2

		@logmsg loglevel3 "Missing $(n_missing) sites still."
		@logmsg loglevel3 "Round2: Inspect box with radius $(Δ_2) and pivot value $(piv_val_1)."

		### as before, only consider points in box of radius `Δ_2`, but ignore `x` and the previous points
		candidate_indices_2 = results_in_box_indices( db, lb_2, ub_2, [x_index; candidate_indices_1])

		filter_2 = AffinelyIndependentPointFilter(;
			x_0 = x,
			seeds = get_site.(db, candidate_indices_2),
			Y = filter.Y,	# pass prior matrices, so that new points are projected onto span of Z
			Z = filter.Z,
			n = n_missing,
			return_indices = true,
			pivot_val = piv_val_2
		)

		filtered_indices_2 = candidate_indices_2[ collect(filter_2) ]

		### Store indices
		empty!(meta.round2_indices)
		append!(meta.round2_indices, filtered_indices_2)

		@logmsg loglevel3 "Round2: Found $(length(meta.round2_indices)) sites and model is $(meta.fully_linear ? "" : "not " )fully linear."
	end

	## Round 3:
	### If we still don't have enough sites, generate them
	### along model improving directions (from first round of sampling)

	n_missing -= length(meta.round2_indices)
	empty!(meta.round3_indices)
	if n_missing > 0

		@logmsg loglevel3 "Round3: Still missing $(n_missing). Sampling in box of radius $(Δ_1)."

		### If round 2 did not yield any new points, the model will hopefully be made fully linear now.
		if length(meta.round2_indices) == 0
			meta.fully_linear = true
		end

		### Take into consideration the maximum number of evaluations allowed:
		### TODO count unevaluated results too! This is a new requirement due to the 2-phase construction process. #src
		max_new = min( max_evals(algo_config), max_evals(cfg) ) - 1 - num_evals( objf )
		n_new = min(n_missing, max_new)

		new_points = Vector{F}[]
		while !isempty(meta.improving_directions) && length( new_points ) < n_new
			dir = popfirst!( meta.improving_directions )
			len = intersect_bounds( x, dir, lb_1, ub_1; return_vals = :absmax )
			offset = len .* dir
			if norm( offset, Inf ) <= piv_val_1
				### the new point does not pass the thresholding test
				if ensure_fully_linear && !force_rebuild
					### If we need a fully linear model, we dismiss the inidices gathered so far …
					### … and call for a rebuild along the coordinate axis:
					return prepare_update_model(mod, objf, meta, mop, iter_data, db, algo_config; ensure_fully_linear = true, force_rebuild = true)
				else
					### we include the point nonetheless, but the model will not qualify as fully linear...
					meta.fully_linear = false
				end
			end
			push!( new_points, x .+ offset )
		end

		### by adding the points to the database at this point in time we avoid
		### requesting unnecessary results from a round 3 interrupted by rebuilding
		new_indices = Int[]
		for p ∈ new_points
			new_id = new_result!( db, p, F[] )
			push!(new_indices, new_id)
		end

		append!(meta.round3_indices, new_indices)
	end

	@label round4

	## In round 4 we have found `n_vars + 1` training sites and try to find additional points within the
	## largest possible trust region.
	empty!(meta.round4_indices)

	if cfg.optimized_sampling

		max_points = cfg.max_model_points <= 0 ? 2 * n_vars + 1 : cfg.max_model_points
		indices_found_so_far = _collect_indices( meta )
		N = length(indices_found_so_far)

		candidate_indices_4 = results_in_box_indices( db, lb_2, ub_2, indices_found_so_far )

		max_tries = 10 * max_points
		num_tries = 0

		if N < max_points && ( !isempty(candidate_indices_4) || cfg.use_max_points )
			@logmsg loglevel3 "Round4: Can we find $(max_points - N) additional sites?"
			round4_indices = Int[]

			chol_pivot = cfg.θ_pivot_cholesky

			centers = get_site.(db, indices_found_so_far)
			φ = _get_radial_function( Δ, cfg )
			Φ, Π, kernels, polys = RBF.get_matrices( φ, centers; poly_deg = cfg.polynomial_degree )

			## prepare matrices as by Wild, R has to be augmented by rows of zeros
			Q, R = qr( transpose(Π) )
			R = [
				R;
				zeros( size(Q,1) - size(R,1), size(R,2) )
			]
			Z = Q[:, N + 1 : end ] ## columns of Z are orthogonal to Π

			## Note: usually, Z, ZΦZ and L should be empty (if N == n_vars + 1)
			ZΦZ = Hermitian(Z'Φ*Z)	## make sure, it is really symmetric
			L = cholesky( ZΦZ ).L   ## perform cholesky factorization
			L⁻¹ = inv(L)				 ## most likely empty at this point

			φ₀ = Φ[1,1]

			@logmsg loglevel3 "Round4: Considering $(length(candidate_indices_4)) candidates."

			while N < max_points && num_tries <= max_tries

				if !isempty( candidate_indices_4 )
					id = popfirst!( candidate_indices_4 )
					### get candidate site ξ ∈ ℝⁿ
					ξ = get_site( db, id )
				else
					if cfg.use_max_points
						### there are no more sites in the db, but we **want**
						### to use as many as possible
						id = -1
						ξ = _rand_box_point( lb_2, ub_2, F)
						num_tries += 1
					else
						break
					end
				end

				### apply all RBF kernels
				φξ = kernels( ξ )

				### apply polynomial basis system and augment polynomial matrix
				πξ = polys( ξ )
				Rξ = [ R; πξ' ]

				### perform Givens rotations to turn last row in Rξ to zeros
				row_index = size( Rξ, 1)
				G = Matrix(I, row_index, row_index) # whole orthogonal matrix
				for j = 1 : size(R,2)
					## in each column, take the diagonal as pivot to turn last elem to zero
					g = givens( Rξ[j,j], Rξ[row_index, j], j, row_index )[1]
					Rξ = g*Rξ;
					G = g*G;
				end

				### now, from G we can update the other matrices
				Gᵀ = transpose(G)
				g̃ = Gᵀ[1 : end-1, end]
				ĝ = Gᵀ[end, end]

				Qg = Q*g̃;
				v_ξ = Z'*( Φ*Qg + φξ .* ĝ )
				σ_ξ = Qg'*Φ*Qg + (2*ĝ) * φξ'*Qg + ĝ^2*φ₀

				τ_ξ² = σ_ξ - norm( L⁻¹ * v_ξ, 2 )^2
				## τ_ξ (and hence τ_ξ^2) must be bounded away from zero
				## for the model to remain fully linear
				if τ_ξ² > chol_pivot

					if id < 0
						id = new_result!( db, ξ, F[] )
					end
					push!(round4_indices, id)	# accept the result

					τ_ξ = sqrt(τ_ξ²)

					## zero-pad Q and multiply with Gᵗ
					Q = [
						Q 					zeros( size(Q,1), 1);
						zeros(1, size(Q,2)) 1
					] * Gᵀ

					Z = [
						Z  						Qg;
						zeros(1, size(Z,2)) 	ĝ
					]

					L = [
						L          zeros(size(L,1), 1) ;
						v_ξ'L⁻¹'   τ_ξ
					]

					L⁻¹ = [
						L⁻¹                zeros(size(L⁻¹,1),1);
						-(v_ξ'L⁻¹'L⁻¹)./τ_ξ   1/τ_ξ
					]

					R = Rξ

					## finally, augment basis matrices and add new kernel for next iteration
					Π = [ Π πξ ]

					Φ = [
						Φ   φξ;
						φξ' φ₀
					]
					push!( kernels, RBF.make_kernel(φ, ξ) )

					## assert all( diag( L * L⁻¹) .≈ 1 )
					N += 1
				end#if
			end#for
			append!(meta.round4_indices, round4_indices)
			@logmsg loglevel3 "Round4: found $(length(round4_indices)) additional sites."
		end#if
	end# if cfg.optimized_sampling

	return meta
end

# !!! note
#     At the moment, we do not store the matrices calculated in round 4 of the
#     update procedure. This could be done to save some work when actually
#     calculating the coefficients.

# In contrast to the old RBF mechanism, the models in `RadialBasisFunctionModels` sometimes
# accept 2 parameters for the kernel. We use this little helper, to get defaults from the shape parameter.
# Note, that sanity check are performed in the RbfConfig constructor.
function _get_kernel_params( Δ , cfg )

	sp = if cfg.shape_parameter isa String
		parse_shape_param_string( Δ, cfg.shape_parameter )
	else
		cfg.shape_parameter
	end

	isnan(sp) && return nothing

	kernel_name = cfg.kernel

	if kernel_name == :gaussian
		return sp
	elseif kernel_name == :inv_multiquadric
		return (sp, 1//2)
	elseif kernel_name == :multiquadric
		return (sp, 1//2)
	elseif kernel_name == :cubic
		return Int(sp)
	elseif kernel_name == :thin_plate_spline
		return Int(sp)
	else
		return sp
	end
end

function _get_radial_function( Δ, cfg )
	kernel_params = _get_kernel_params( Δ, cfg )

	return RBF._get_rad_func( cfg.kernel, kernel_params )
end

# An improvement step consists of adding a new site to the database, along an improving direction:
function prepare_improve_model( mod :: Union{Nothing, RbfModel}, objf :: AbstractObjective,
	meta :: RbfMeta, mop :: AbstractMOP, iter_data :: AbstractIterData, db :: AbstractDB,
	algo_config :: AbstractConfig; kwargs... )
	if !meta.fully_linear
		if isempty(meta.improving_directions)
			@warn "RBF model is not fully linear, but there are no improving directions."
		else
			cfg = model_cfg(objf)
			x = get_x(iter_data)
			fx = get_fx(iter_data)
			F = typeof(fx)
			Δ = get_delta(iter_data)
			Δ_1 = Δ * cfg.θ_enlarge_1
			lb_1, ub_1 = local_bounds(mop, x, Δ_1)
			piv_val_1 = Δ_1 * cfg.θ_pivot

			success = false
			dir = popfirst!( meta.improving_directions )
			len = intersect_bounds( x, dir, lb_1, ub_1; return_vals = :absmax )
			offset = len .* dir
			if norm( offset, Inf ) > piv_val_1
				new_id = new_result!( db, x .+ offset, F() )
				push!(meta.round1_indices, new_id)
				success = true
			end

			success && @logmsg loglevel3 "Performed an improvement step."
			if isempty( meta.improving_directions ) && success
				meta.fully_linear = true
				@logmsg loglevel3 "The RBF Model is now fully linear."
			end
		end
	end
	return meta
end

# Now, in the 2-phase construction process, first all `prepare_` functions are called for all surrogate models.
# Then, the unevaluated results are evaluated and we can proceed with the model building.
# As before, `_init_model` simply delegates work to `update_model`.

function _init_model( cfg :: RbfConfig, objf :: AbstractObjective, mop :: AbstractMOP,
	iter_data :: AbstractIterData, db :: AbstractDB, ac :: AbstractConfig, meta :: RbfMeta; kwargs... )
	return update_model( nothing, objf, meta, mop, iter_data, db, ac; kwargs... )
end

function update_model( mod::Union{Nothing,RbfModel}, objf:: AbstractObjective, meta :: RbfMeta,
	mop :: AbstractMOP, iter_data :: AbstractIterData, db :: AbstractDB, ac :: AbstractConfig;
	kwargs... )

	Δ = get_delta(iter_data)
	cfg = model_cfg( objf )

	kernel_params = _get_kernel_params( Δ, cfg )

	## get the training data from `meta` and the database `db`
	training_indices = _collect_indices( meta )
	training_results = get_result.( db, training_indices )
	training_sites = get_site.( training_results )
	oi = output_indices( objf, mop)	# only consider the objective output indices
	training_values = [ v[oi] for v in get_value.( training_results ) ]

	inner_model = RBF.RBFInterpolationModel( training_sites, training_values, cfg.kernel, kernel_params; save_matrices = false )

	@logmsg loglevel3 "The model is $(meta.fully_linear ? "" : "not ")fully linear."
	return RbfModel( inner_model, meta.fully_linear ), meta

end

# The improvement function also simply cals the update function:
function improve_model( mod::Union{Nothing,RbfModel}, objf:: AbstractObjective, meta :: RbfMeta,
	mop :: AbstractMOP, id :: AbstractIterData, db :: AbstractDB, ac :: AbstractConfig; kwargs... )
	return update_model( mod, objf, meta, mop, id, db, ac; kwargs... )
end

# ## Evaluation

# All the work is done by the inner model :)

"Evaluate `mod::RbfModel` at scaled site `x̂`."
function eval_models( mod :: RbfModel, x̂ :: Vec)
	return mod.model( x̂ )
end

"Evaluate output `ℓ` of `mod::RbfModel` at scaled site `x̂`."
function eval_models( mod :: RbfModel, x̂ :: Vec, ℓ :: Int)
	return mod.model( x̂, ℓ)
end

@doc "Gradient vector of output `ℓ` of `mod` at scaled site `x̂`."
function get_gradient( mod :: RbfModel, x̂ :: Vec, ℓ :: Int64)
    return RBF.grad( mod.model, x̂, ℓ )
end

@doc "Jacobian Matrix of ExactModel `em` at scaled site `x̂`."
function get_jacobian( mod :: RbfModel, x̂ :: Vec )
    return RBF.jac( mod.model, x̂ )
end

# ## [Summary & Quick Examples](@id rbf_summary)

# 1. To use the default configuration for a scalar objective `f` do
#    ```julia
#    add_objective!(mop, f, RbfConfig())
#    ```
# 2. For a vector valued objective do
#    ```julia
#    add_vector_objective!(mop, f, RbfConfig(); n_out = 2)
#    ```
# 3. If you don't want to use a polynomial:
#    ```julia
#    add_objective!(mop, f, RbfConfig(;kernel = :cubic, polynomial_degree = -1 ))
#    ```
#    This only works for certain kernels. `polynomial_degree = 0` will add a constant term.
# 4. To require sampling of the maximum number of allowed model points:
#    ```julia
#    RbfConfig(;use_max_points = true)
#    ```
# 5. To only sample along the coordinate axis:
#    ```julia
#    RbfConfig(;optimized_sampling = false)
#    ```
#    If `polynomial_degree == 1` the model will now be a linear interpolation model.

# ### Complete usage example
# ```julia
# using Morbit
# Morbit.print_all_logs()
# mop = MixedMOP(3)
#
# F = x -> [ sum( ( x .- 1 ).^2 ); sum( ( x .+ 1 ).^2 ) ]
#
# add_vector_objective!( mop, F, RbfConfig() )
#
# x_fin, f_fin, _ = optimize( mop, [-π, ℯ, 0])
# ```

# [^wild_diss]: “Derivative-Free Optimization Algorithms For Computationally Expensive Functions”, Stefan M. Wild, 2009


text:
using MIDI, Distributions

mutable struct Engine
    key::Key
    note_generator::Distribution
    beat_generator::Distribution
end

function generate_note(key::Key, initial::Note, generator::Distribution)::Note
    delta = key.ladder .- initial
    idx = findmin(abs.([d.steps for d in delta] .- trunc(Int, rand(generator))))[2]; increment = delta[idx]
    out = initial + delta[idx]
    out.position += initial.duration
    return out
end

function generate_notes(key::Key, initial::Note, generator::Distribution, n::Int)::Notes{Note}
    out = Notes()
    push!(out, initial)
    for i in 2:n
        push!(out, generate_note(key, out[i-1], generator))
    end
    return out
end

function generate_rhythm(notes::Notes{Note}, generator::Distribution, precision::Rational=SIXTEENTH)::Notes{Note}
    beats = round.(1 .// [ceil(Int, max(rand(generator), 1)) for i in 1:length(notes)], digits=1, base=Int(1/precision))
    notes[1].duration = beats[1]*TPQ
    for i in 2:length(notes)
        notes[i].duration = beats[i]*TPQ
        notes[i].position = notes[i-1].duration + notes[i-1].position
    end
    return notes
end

function generate(engine::Engine, n::Int, initial::Note=engine.key.ladder[rand(1:length(engine.key.ladder))], precision::Rational=SIXTEENTH)::Notes{Note}
    notes = generate_notes(engine.key, initial, engine.note_generator, n)
    notes = generate_rhythm(notes, engine.beat_generator, precision)
    return notes
end


text:
test_msh = (
    msh_version = "4.1",
    raw_PhysicalNames = [
        "9",
        "0 1 Point1",
        "0 2 Point2",
        "0 3 Point3",
        "0 4 Point4",
        "1 5 Curve1",
        "1 6 Curve2",
        "1 7 Curve3",
        "1 8 Curve4",
        "2 9 Surface1"
    ],
    raw_Entities = SubString{String}[
        "4 4 1 0",
        "1 0 0 0 1 1 ",
        "2 1 0 0 1 2 ",
        "3 1 1 0 1 3 ",
        "4 0 1 0 1 4 ",
        "1 0 0 0 1 0 0 1 5 2 1 -2 ",
        "2 1 0 0 1 1 0 1 6 2 2 -3 ",
        "3 0 1 0 1 1 0 1 7 2 3 -4 ",
        "4 0 0 0 0 1 0 1 8 2 4 -1 ",
        "1 0 0 0 1 1 0 1 9 4 1 2 3 4 "
    ],
    raw_Nodes = SubString{String}[
        "9 30 1 30",
        "0 1 0 1",
        "1",
        "0 0 0",
        "0 2 0 1",
        "2",
        "1 0 0",
        "0 3 0 1",
        "3",
        "1 1 0",
        "0 4 0 1",
        "4",
        "0 1 0",
        "1 1 0 3",
        "5",
        "6",
        "7",
        "0.2499999999994121 0 0",
        "0.499999999998694 0 0",
        "0.7499999999993416 0 0",
        "1 2 0 3",
        "8",
        "9",
        "10",
        "1 0.2499999999994121 0",
        "1 0.499999999998694 0",
        "1 0.7499999999993416 0",
        "1 3 0 3",
        "11",
        "12",
        "13",
        "0.7500000000003465 1 0",
        "0.5000000000020591 1 0",
        "0.2500000000010404 1 0",
        "1 4 0 3",
        "14",
        "15",
        "16",
        "0 0.7500000000003465 0",
        "0 0.5000000000020591 0",
        "0 0.2500000000010404 0",
        "2 1 0 14",
        "17",
        "18",
        "19",
        "20",
        "21",
        "22",
        "23",
        "24",
        "25",
        "26",
        "27",
        "28",
        "29",
        "30",
        "0.3640932128839348 0.7867687832230399 0",
        "0.2110423218192967 0.3753337416957934 0",
        "0.6146003733357942 0.2544869979106268 0",
        "0.8080509614599705 0.6243689539557383 0",
        "0.1979210047255827 0.6215671767908416 0",
        "0.4308090314147045 0.5056502726999197 0",
        "0.3712961061742487 0.2167990069479849 0",
        "0.6139924015723897 0.7478815024664035 0",
        "0.8112741566735278 0.3855439771726542 0",
        "0.6557453848912773 0.5035863408410685 0",
        "0.1830127018930926 0.816987298107279 0",
        "0.816987298107459 0.8169872981071898 0",
        "0.1830127018926974 0.1830127018931337 0",
        "0.8351749060017326 0.1780061950165386 0"
    ],
    raw_Elements = SubString{String}[
        "9 62 1 62",
        "0 1 15 1",
        "1 1 ",
        "0 2 15 1",
        "2 2 ",
        "0 3 15 1",
        "3 3 ",
        "0 4 15 1",
        "4 4 ",
        "1 1 1 4",
        "5 1 5 ",
        "6 5 6 ",
        "7 6 7 ",
        "8 7 2 ",
        "1 2 1 4",
        "9 2 8 ",
        "10 8 9 ",
        "11 9 10 ",
        "12 10 3 ",
        "1 3 1 4",
        "13 3 11 ",
        "14 11 12 ",
        "15 12 13 ",
        "16 13 4 ",
        "1 4 1 4",
        "17 4 14 ",
        "18 14 15 ",
        "19 15 16 ",
        "20 16 1 ",
        "2 1 2 42",
        "21 19 22 23 ",
        "22 17 22 24 ",
        "23 22 19 26 ",
        "24 24 22 26 ",
        "25 22 18 23 ",
        "26 17 21 22 ",
        "27 20 24 26 ",
        "28 12 13 17 ",
        "29 15 16 18 ",
        "30 9 10 20 ",
        "31 6 7 19 ",
        "32 19 25 26 ",
        "33 14 15 21 ",
        "34 11 12 24 ",
        "35 5 6 23 ",
        "36 8 9 25 ",
        "37 21 18 22 ",
        "38 15 18 21 ",
        "39 12 17 24 ",
        "40 6 19 23 ",
        "41 9 20 25 ",
        "42 4 14 27 ",
        "43 3 11 28 ",
        "44 13 4 27 ",
        "45 16 1 29 ",
        "46 1 5 29 ",
        "47 10 3 28 ",
        "48 7 2 30 ",
        "49 2 8 30 ",
        "50 14 21 27 ",
        "51 11 24 28 ",
        "52 5 23 29 ",
        "53 8 25 30 ",
        "54 17 13 27 ",
        "55 20 10 28 ",
        "56 18 16 29 ",
        "57 19 7 30 ",
        "58 21 17 27 ",
        "59 23 18 29 ",
        "60 24 20 28 ",
        "61 25 19 30 ",
        "62 25 20 26 "
    ]
)

test_physicals = (
    physicalTag0DtoName = Dict(
        4 => "Point4",
        2 => "Point2",
        3 => "Point3",
        1 => "Point1"
    ),
    physicalTag1DtoName = Dict(
        7 => "Curve3",
        5 => "Curve1",
        8 => "Curve4",
        6 => "Curve2"
    ),
    physicalTag2DtoName = Dict(
        9 => "Surface1"
    ),
    physicalTag3DtoName = Dict{Int64,String}()
)

test_entities = (
    pointTagtoPhysicalTag = Dict(
        4 => 4,
        2 => 2,
        3 => 3,
        1 => 1
    ),
    curveTagtoPhysicalTag = Dict(
        4 => 8,
        2 => 6,
        3 => 7,
        1 => 5
    ),
    surfaceTagtoPhysicalTag = Dict(
        1 => 9
    ),
    volumeTagtoPhysicalTag = Dict{Int64,Int64}()
)

test_nodes = (
    nodeTag = [
        1,
        2,
        3,
        4,
        5,
        6,
        7,
        8,
        9,
        10,
        11,
        12,
        13,
        14,
        15,
        16,
        17,
        18,
        19,
        20,
        21,
        22,
        23,
        24,
        25,
        26,
        27,
        28,
        29,
        30
    ],
    vx = [
        0.0,
        1.0,
        1.0,
        0.0,
        0.2499999999994121,
        0.499999999998694,
        0.7499999999993416,
        1.0,
        1.0,
        1.0,
        0.7500000000003465,
        0.5000000000020591,
        0.2500000000010404,
        0.0,
        0.0,
        0.0,
        0.3640932128839348,
        0.2110423218192967,
        0.6146003733357942,
        0.8080509614599705,
        0.1979210047255827,
        0.4308090314147045,
        0.3712961061742487,
        0.6139924015723897,
        0.8112741566735278,
        0.6557453848912773,
        0.1830127018930926,
        0.816987298107459,
        0.1830127018926974,
        0.8351749060017326
    ],
    vy = [
        0.0,
        0.0,
        1.0,
        1.0,
        0.0,
        0.0,
        0.0,
        0.2499999999994121,
        0.499999999998694,
        0.7499999999993416,
        1.0,
        1.0,
        1.0,
        0.7500000000003465,
        0.5000000000020591,
        0.2500000000010404,
        0.7867687832230399,
        0.3753337416957934,
        0.2544869979106268,
        0.6243689539557383,
        0.6215671767908416,
        0.5056502726999197,
        0.2167990069479849,
        0.7478815024664035,
        0.3855439771726542,
        0.5035863408410685,
        0.816987298107279,
        0.8169872981071898,
        0.1830127018931337,
        0.1780061950165386
    ],
    vz = [
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0
    ]
)

test_elements = (
    entityTag1D = [
        1,
        1,
        1,
        1,
        2,
        2,
        2,
        2,
        3,
        3,
        3,
        3,
        4,
        4,
        4,
        4
    ],
    entityTag2D = [
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1
    ],
    entityTag3D = Int64[],
    elementType1D = [
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1,
        1
    ],
    elementType2D = [
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2,
        2
    ],
    elementType3D = Int64[],
    elementTag1D = [
        5,
        6,
        7,
        8,
        9,
        10,
        11,
        12,
        13,
        14,
        15,
        16,
        17,
        18,
        19,
        20
    ],
    elementTag2D = [
        21,
        22,
        23,
        24,
        25,
        26,
        27,
        28,
        29,
        30,
        31,
        32,
        33,
        34,
        35,
        36,
        37,
        38,
        39,
        40,
        41,
        42,
        43,
        44,
        45,
        46,
        47,
        48,
        49,
        50,
        51,
        52,
        53,
        54,
        55,
        56,
        57,
        58,
        59,
        60,
        61,
        62
    ],
    elementTag3D = Int64[],
    nodeTags1D = [
        [1, 5],
        [5, 6],
        [6, 7],
        [7, 2],
        [2, 8],
        [8, 9],
        [9, 10],
        [10, 3],
        [3, 11],
        [11, 12],
        [12, 13],
        [13, 4],
        [4, 14],
        [14, 15],
        [15, 16],
        [16, 1]
    ],
    nodeTags2D = [
        [19, 22, 23],
        [17, 22, 24],
        [22, 19, 26],
        [24, 22, 26],
        [22, 18, 23],
        [17, 21, 22],
        [20, 24, 26],
        [12, 13, 17],
        [15, 16, 18],
        [9, 10, 20],
        [6, 7, 19],
        [19, 25, 26],
        [14, 15, 21],
        [11, 12, 24],
        [5, 6, 23],
        [8, 9, 25],
        [21, 18, 22],
        [15, 18, 21],
        [12, 17, 24],
        [6, 19, 23],
        [9, 20, 25],
        [4, 14, 27],
        [3, 11, 28],
        [13, 4, 27],
        [16, 1, 29],
        [1, 5, 29],
        [10, 3, 28],
        [7, 2, 30],
        [2, 8, 30],
        [14, 21, 27],
        [11, 24, 28],
        [5, 23, 29],
        [8, 25, 30],
        [17, 13, 27],
        [20, 10, 28],
        [18, 16, 29],
        [19, 7, 30],
        [21, 17, 27],
        [23, 18, 29],
        [24, 20, 28],
        [25, 19, 30],
        [25, 20, 26]
    ],
    nodeTags3D = Array{Int64,1}[]
)

using MeshCore.Exports
using StaticArrays
test_xyz = reshape([test_nodes.vx; test_nodes.vy], length(test_nodes.vx), 2)
test_N, test_T = size(test_xyz, 2), eltype(test_xyz)
test_locs = VecAttrib([SVector{test_N,test_T}(test_xyz[i, :]) for i in 1:size(test_xyz, 1)])
test_vertices = ShapeColl(P1, length(test_locs), "vertices")
test_vertices.attributes["geom"] = test_locs
test_shapes = ShapeColl(T3, size(test_elements.nodeTags2D, 1), "elements")
test_ir20 = IncRel(test_shapes, test_vertices, test_elements.nodeTags2D)