text:
#####
##### Differential operators for regular grids
#####

@inline ∂x_caa(i, j, k, grid, u::AbstractArray) = @inbounds (u[i+1, j, k] - u[i, j, k]) / grid.Δx
@inline ∂x_faa(i, j, k, grid, c::AbstractArray) = @inbounds (c[i, j, k] - c[i-1, j, k]) / grid.Δx

@inline ∂y_aca(i, j, k, grid, v::AbstractArray) = @inbounds (v[i, j+1, k] - v[i, j, k]) / grid.Δy
@inline ∂y_afa(i, j, k, grid, c::AbstractArray) = @inbounds (c[i, j, k] - c[i, j-1, k]) / grid.Δy

@inline ∂z_aac(i, j, k, grid, w::AbstractArray) = @inbounds (w[i, j, k] - w[i, j, k+1]) / grid.Δz
@inline ∂z_aaf(i, j, k, grid, c::AbstractArray) = @inbounds (c[i, j, k-1] - c[i, j, k]) / grid.Δz

#
# Differentiation and interpolation operators for functions
#
# The geometry of the Oceananigans grid (in one dimension) is
#
# face   cell   face   cell   face
#
#         i-1            i
#          ↓             ↓
#   |      ×      |      ×      |
#   ↑             ↑             ↑
#  i-1            i            i+1
#
# As a result the interpolation of a quantity cᵢ from a cell i to face i
# (this operation is denoted "▶x_faa" in the code below) is
#
# (▶c)ᵢ = (cᵢ + cᵢ₋₁) / 2 .
#
# Conversely, the interpolation of a quantity uᵢ from face i to cell i is given by
#
# (▶u)ᵢ = (uᵢ₊₁ + uᵢ) / 2.
#
# Derivative operators are defined similarly. Using the symbol "∂" to denote
# differentiation, we have
#
# (∂c)ᵢ = (cᵢ - cᵢ₋₁) / Δ
#
# for the differentation of a CellField ending on face i, and
#
# (∂u)ᵢ = (uᵢ₊₁ - uᵢ) / Δ
#
# for the differentation of a FaceField ending on cell i, where Δ is the spacing
# between cell points or face points, appropriately.
#

#####
##### Differential operators
#####

"""
    ∂x_faa(i, j, k, grid, F, args...)

Differentiate the function or callable object

    `F(i, j, k, grid, args...)`

located at `caa` in `x`, across `faa`.
"""
@inline ∂x_faa(i, j, k, grid, F::TF, args...) where TF<:Function =
    (F(i, j, k, grid, args...) - F(i-1, j, k, grid, args...)) / grid.Δx

"""
    ∂x_caa(i, j, k, grid, F, args...)

Differentiate the function or callable object

    `F(i, j, k, grid, args...)`

located at `faa` in `x`, across `caa`.
"""
@inline ∂x_caa(i, j, k, grid, F::TF, args...) where TF<:Function =
    (F(i+1, j, k, grid, args...) - F(i, j, k, grid, args...)) / grid.Δx

"""
    ∂y_afa(i, j, k, grid, F, args...)

Differentiate the function or callable object

    `F(i, j, k, grid, args...)`

located at `aca` in `y`, across `afa`.
"""
@inline ∂y_afa(i, j, k, grid, F::TF, args...) where TF<:Function =
    (F(i, j, k, grid, args...) - F(i, j-1, k, grid, args...)) / grid.Δy

"""
    ∂y_aca(i, j, k, grid, F, args...)

Differentiate the function or callable object

    `F(i, j, k, grid, args...)`

located at `afa` in `y`, across `aca`.
"""
@inline ∂y_aca(i, j, k, grid, F::TF, args...) where TF<:Function =
    (F(i, j+1, k, grid, args...) - F(i, j, k, grid, args...)) / grid.Δy

"""
    ∂z_aaf(i, j, k, grid, F, args...)

Differentiate the function or callable object

    `F(i, j, k, grid, args...)`

located at `aac` in `z`, across `aaf`.
"""
@inline ∂z_aaf(i, j, k, grid::AbstractGrid, F::TF, args...) where TF<:Function =
    (F(i, j, k-1, grid, args...) - F(i, j, k, grid, args...)) / grid.Δz

"""
    ∂z_aac(i, j, k, grid, F, args...)

Differentiate the function or callable object

    `F(i, j, k, grid, args...)`

located at `aaf` in `z`, across `aac`.
"""
@inline ∂z_aac(i, j, k, grid, F::TF, args...) where TF<:Function =
    (F(i, j, k, grid, args...) - F(i, j, k+1, grid, args...)) / grid.Δz

#####
##### Double differentiation
#####

@inline ∂x²_caa(i, j, k, grid, c) = ∂x_caa(i, j, k, grid, ∂x_faa, c)
@inline ∂x²_faa(i, j, k, grid, u) = ∂x_faa(i, j, k, grid, ∂x_caa, u)

@inline ∂y²_aca(i, j, k, grid, c) = ∂y_aca(i, j, k, grid, ∂y_afa, c)
@inline ∂y²_afa(i, j, k, grid, v) = ∂y_afa(i, j, k, grid, ∂y_aca, v)

@inline ∂z²_aac(i, j, k, grid, c) = ∂z_aac(i, j, k, grid, ∂z_aaf, c)
@inline ∂z²_aaf(i, j, k, grid, w) = ∂z_aaf(i, j, k, grid, ∂z_aac, w)

#####
##### Interpolation operations for functions
#####

"""
    ▶x_faa(i, j, k, grid, F, args...)

Interpolate the function or callable object

    `F(i, j, k, grid, args...)`

from `caa` to `faa`."
"""
@inline ▶x_faa(i, j, k, grid::RegularCartesianGrid{T}, F::TF, args...) where {T, TF<:Function} =
    T(0.5) * (F(i, j, k, grid, args...) + F(i-1, j, k, grid, args...))

"""
    ▶x_caa(i, j, k, grid, F, args...)

Interpolate the function or callable object

    `F(i, j, k, grid, args...)`

from `faa` to `caa`."
"""
@inline ▶x_caa(i, j, k, grid::RegularCartesianGrid{T}, F::TF, args...) where {T, TF<:Function} =
    return T(0.5) * (F(i+1, j, k, grid, args...) + F(i, j, k, grid, args...))

"""
    ▶y_afa(i, j, k, grid, F, args...)

Interpolate the function or callable object

    `F(i, j, k, grid, args...)`

from `aca` to `afa`.
"""
@inline ▶y_afa(i, j, k, grid::RegularCartesianGrid{T}, F::TF, args...) where {T, TF<:Function} =
    return T(0.5) * (F(i, j, k, grid, args...) + F(i, j-1, k, grid, args...))

"""
    ▶y_aca(i, j, k, grid, F, args...)

Interpolate the function or callable object

    `F(i, j, k, grid, args...)`

from `afa` to `aca`."
"""
@inline ▶y_aca(i, j, k, grid::RegularCartesianGrid{T}, F::TF, args...) where {T, TF<:Function} =
    T(0.5) * (F(i, j+1, k, grid, args...) + F(i, j, k, grid, args...))

"""
    ▶z_aaf(i, j, k, grid, F, args...)

Interpolate the function or callable object

    `F(i, j, k, grid, args...)`

from `aac` to `aaf`.
"""
@inline ▶z_aaf(i, j, k, grid::RegularCartesianGrid{T}, F::TF, args...) where {T, TF<:Function} =
    T(0.5) * (F(i, j, k, grid, args...) + F(i, j, k-1, grid, args...))

"""
    ▶z_aac(i, j, k, grid::RegularCartesianGrid{T}, F, args...) where T

Interpolate the function or callable object

    `F(i, j, k, grid, args...)`

from `aaf` to `aac`.
"""
@inline ▶z_aac(i, j, k, grid::RegularCartesianGrid{T}, F::TF, args...) where {T, TF<:Function} =
    T(0.5) * (F(i, j, k+1, grid, args...) + F(i, j, k, grid, args...))

# Convenience operators for "interpolating constants"
@inline ▶x_faa(i, j, k, grid, F::Number, args...) = F
@inline ▶x_caa(i, j, k, grid, F::Number, args...) = F
@inline ▶y_afa(i, j, k, grid, F::Number, args...) = F
@inline ▶y_aca(i, j, k, grid, F::Number, args...) = F
@inline ▶z_aaf(i, j, k, grid, F::Number, args...) = F
@inline ▶z_aac(i, j, k, grid, F::Number, args...) = F

@inline ▶x_faa(i, j, k, grid::RegularCartesianGrid{T}, F::AbstractArray, args...) where T =
    @inbounds T(0.5) * (F[i, j, k] + F[i-1, j, k])

@inline ▶x_caa(i, j, k, grid::RegularCartesianGrid{T}, F::AbstractArray, args...) where T =
    @inbounds T(0.5) * (F[i, j, k] + F[i+1, j, k])

@inline ▶y_afa(i, j, k, grid::RegularCartesianGrid{T}, F::AbstractArray, args...) where T =
    @inbounds T(0.5) * (F[i, j, k] + F[i, j-1, k])

@inline ▶y_aca(i, j, k, grid::RegularCartesianGrid{T}, F::AbstractArray, args...) where T =
    @inbounds T(0.5) * (F[i, j, k] + F[i, j+1, k])

@inline ▶z_aaf(i, j, k, grid::RegularCartesianGrid{T}, F::AbstractArray, args...) where T =
    @inbounds T(0.5) * (F[i, j, k] + F[i, j, k-1])

@inline ▶z_aac(i, j, k, grid::RegularCartesianGrid{T}, w::AbstractArray, args...) where T =
    @inbounds T(0.5) * (w[i, j, k] + w[i, j, k+1])

#####
##### Double interpolation: 12 operators
#####

"""
    ▶xy_cca(i, j, k, grid, F, args...)

Interpolate the function

    `F(i, j, k, grid, args...)`

in `x` and `y`, from `ffa` to `cca`.
"""
@inline ▶xy_cca(i, j, k, grid, F, args...) = ▶y_aca(i, j, k, grid, ▶x_caa, F, args...)

"""
    ▶xy_fca(i, j, k, grid, F, args...)

Interpolate the function

    `F(i, j, k, grid, args...)`

in `x` and `y`, from `cfa` to `fca`.
"""
@inline ▶xy_fca(i, j, k, grid, F, args...) = ▶y_aca(i, j, k, grid, ▶x_faa, F, args...)

"""
    ▶xy_ffa(i, j, k, grid, F, args...)

Interpolate the function

    `F(i, j, k, grid, args...)`

in `x` and `y`, from `cca` to `ffa`.
"""
@inline ▶xy_ffa(i, j, k, grid, F, args...) = ▶y_afa(i, j, k, grid, ▶x_faa, F, args...)

"""
    ▶xy_cfa(i, j, k, grid, F, args...)

Interpolate the function

    `F(i, j, k, grid, args...)`

in `x` and `y`, from `fca` to `cfa`.
"""
@inline ▶xy_cfa(i, j, k, grid, F, args...) = ▶y_afa(i, j, k, grid, ▶x_caa, F, args...)

"""
    ▶xz_cac(i, j, k, grid, F, args...)

Interpolate the function

    `F(i, j, k, grid, args...)`

in `x` and `z`, from `faf` to `cac`.
"""
@inline ▶xz_cac(i, j, k, grid, F, args...) = ▶z_aac(i, j, k, grid, ▶x_caa, F, args...)

"""
    ▶xz_fac(i, j, k, grid, F, args...)

Interpolate the function

    `F(i, j, k, grid, args...)`

in `x` and `z`, from `caf` to `fac`.
"""
@inline ▶xz_fac(i, j, k, grid, F, args...) = ▶z_aac(i, j, k, grid, ▶x_faa, F, args...)

"""
    ▶xz_faf(i, j, k, grid, F, args...)

Interpolate the function

    `F(i, j, k, grid, args...)`

in `x` and `z`, from `cac` to `faf`.
"""
@inline ▶xz_faf(i, j, k, grid, F, args...) = ▶z_aaf(i, j, k, grid, ▶x_faa, F, args...)

"""
    ▶xz_caf(i, j, k, grid, F, args...)

Interpolate the function

    `F(i, j, k, grid, args...)`

in `x` and `z`, from `fac` to `caf`.
"""
@inline ▶xz_caf(i, j, k, grid, F, args...) = ▶z_aaf(i, j, k, grid, ▶x_caa, F, args...)

"""
    ▶yz_acc(i, j, k, grid, F, args...)

Interpolate the function

    `F(i, j, k, grid, args...)`

in `y` and `z`, from `aff` to `acc`.
"""
@inline ▶yz_acc(i, j, k, grid, F, args...) = ▶z_aac(i, j, k, grid, ▶y_aca, F, args...)

"""
    ▶yz_afc(i, j, k, grid, F, args...)

Interpolate the function

    `F(i, j, k, grid, args...)`

in `y` and `z`, from `acf` to `afc`.
"""
@inline ▶yz_afc(i, j, k, grid, F, args...) = ▶z_aac(i, j, k, grid, ▶y_afa, F, args...)

"""
    ▶yz_aff(i, j, k, grid, F, args...)

Interpolate the function

    `F(i, j, k, grid, args...)`

in `y` and `z`, from `ffa` to `cca`.
"""
@inline ▶yz_aff(i, j, k, grid, F, args...) = ▶z_aaf(i, j, k, grid, ▶y_afa, F, args...)

"""
    ▶yz_acf(i, j, k, grid, F, args...)

Interpolate the function

    `F(i, j, k, grid, args...)`

in `y` and `z`, from `afc` to `acf`.
"""
@inline ▶yz_acf(i, j, k, grid, F, args...) = ▶z_aaf(i, j, k, grid, ▶y_aca, F, args...)

"""
    ▶xyc_ffc(i, j, k, grid, F, args...)

Interpolate the function

    `F(i, j, k, grid, args...)`

in `x`, `y`, and `z` from `ccf` to `ffc`.
"""
@inline ▶xyz_ffc(i, j, k, grid, F, args...) = ▶x_faa(i, j, k, grid, ▶y_afa, ▶z_aac, F, args...)
@inline ▶xyz_ccf(i, j, k, grid, F, args...) = ▶x_caa(i, j, k, grid, ▶y_aca, ▶z_aaf, F, args...)

"""
    ν_Σᵢⱼ(i, j, k, grid, ν, Σᵢⱼ, closure, buoyancy, u, v, w, T, S)

Multiply the viscosity function

    `ν(i, j, k, grid, closure, buoyancy, u, v, w, T, S)`

with the strain tensor component function

    `Σᵢⱼ(i, j, k, grid, u, v, w)`

at index `i, j, k`.
"""
@inline ν_Σᵢⱼ(i, j, k, grid, ν::TN, Σᵢⱼ::TS, closure, buoyancy, u, v, w, T, S) where {TN, TS} =
    ν(i, j, k, grid, closure, buoyancy, u, v, w, T, S) * Σᵢⱼ(i, j, k, grid, u, v, w)

@inline ν_Σᵢⱼ_ccc(i, j, k, grid, ν::TN, Σᵢⱼ::TS, u, v, w) where {TN<:AbstractArray, TS} =
    @inbounds ν[i, j, k] * Σᵢⱼ(i, j, k, grid, u, v, w)

@inline ν_Σᵢⱼ_ffc(i, j, k, grid, ν::TN, Σᵢⱼ::TS, u, v, w) where {TN<:AbstractArray, TS} =
    @inbounds ▶xy_ffa(i, j, k, grid, ν) * Σᵢⱼ(i, j, k, grid, u, v, w)

@inline ν_Σᵢⱼ_fcf(i, j, k, grid, ν::TN, Σᵢⱼ::TS, u, v, w) where {TN<:AbstractArray, TS} =
    @inbounds ▶xz_faf(i, j, k, grid, ν) * Σᵢⱼ(i, j, k, grid, u, v, w)

@inline ν_Σᵢⱼ_cff(i, j, k, grid, ν::TN, Σᵢⱼ::TS, u, v, w) where {TN<:AbstractArray, TS} =
    @inbounds ▶yz_aff(i, j, k, grid, ν) * Σᵢⱼ(i, j, k, grid, u, v, w)

#####
##### Stress divergences
#####

# At fcc
@inline ∂x_2ν_Σ₁₁(i, j, k, grid, closure, u, v, w, diffusivities) =
    2 * ∂x_faa(i, j, k, grid, ν_Σᵢⱼ_ccc, diffusivities.νₑ, Σ₁₁, u, v, w)

@inline ∂y_2ν_Σ₁₂(i, j, k, grid, closure, u, v, w, diffusivities) =
    2 * ∂y_aca(i, j, k, grid, ν_Σᵢⱼ_ffc, diffusivities.νₑ, Σ₁₂, u, v, w)

@inline ∂z_2ν_Σ₁₃(i, j, k, grid, closure, u, v, w, diffusivities) =
    2 * ∂z_aac(i, j, k, grid, ν_Σᵢⱼ_fcf, diffusivities.νₑ, Σ₁₃, u, v, w)

# At cfc
@inline ∂x_2ν_Σ₂₁(i, j, k, grid, closure, u, v, w, diffusivities) =
    2 * ∂x_caa(i, j, k, grid, ν_Σᵢⱼ_ffc, diffusivities.νₑ, Σ₂₁, u, v, w)

@inline ∂y_2ν_Σ₂₂(i, j, k, grid, closure, u, v, w, diffusivities) =
    2 * ∂y_afa(i, j, k, grid, ν_Σᵢⱼ_ccc, diffusivities.νₑ, Σ₂₂, u, v, w)

@inline ∂z_2ν_Σ₂₃(i, j, k, grid, closure, u, v, w, diffusivities) =
    2 * ∂z_aac(i, j, k, grid, ν_Σᵢⱼ_cff, diffusivities.νₑ, Σ₂₃, u, v, w)

# At ccf
@inline ∂x_2ν_Σ₃₁(i, j, k, grid, closure, u, v, w, diffusivities) =
    2 * ∂x_caa(i, j, k, grid, ν_Σᵢⱼ_fcf, diffusivities.νₑ, Σ₃₁, u, v, w)

@inline ∂y_2ν_Σ₃₂(i, j, k, grid, closure, u, v, w, diffusivities) =
    2 * ∂y_aca(i, j, k, grid, ν_Σᵢⱼ_cff, diffusivities.νₑ, Σ₃₂, u, v, w)

@inline ∂z_2ν_Σ₃₃(i, j, k, grid, closure, u, v, w, diffusivities) =
    2 * ∂z_aaf(i, j, k, grid, ν_Σᵢⱼ_ccc, diffusivities.νₑ, Σ₃₃, u, v, w)

#
# Without precomputed diffusivities
#

@inline ∂x_2ν_Σ₁₁(i, j, k, grid, closure, buoyancy, u, v, w, T, S) =
    2 * ∂x_faa(i, j, k, grid, ν_Σᵢⱼ, ν_ccc, Σ₁₁, closure, buoyancy, u, v, w, T, S)

@inline ∂y_2ν_Σ₁₂(i, j, k, grid, closure, buoyancy, u, v, w, T, S) =
    2 * ∂y_aca(i, j, k, grid, ν_Σᵢⱼ, ν_ffc, Σ₁₂, closure, buoyancy, u, v, w, T, S)

@inline ∂z_2ν_Σ₁₃(i, j, k, grid, closure, buoyancy, u, v, w, T, S) =
    2 * ∂z_aac(i, j, k, grid, ν_Σᵢⱼ, ν_fcf, Σ₁₃, closure, buoyancy, u, v, w, T, S)

# At cfc
@inline ∂x_2ν_Σ₂₁(i, j, k, grid, closure, buoyancy, u, v, w, T, S) =
    2 * ∂x_caa(i, j, k, grid, ν_Σᵢⱼ, ν_ffc, Σ₂₁, closure, buoyancy, u, v, w, T, S)

@inline ∂y_2ν_Σ₂₂(i, j, k, grid, closure, buoyancy, u, v, w, T, S) =
    2 * ∂y_afa(i, j, k, grid, ν_Σᵢⱼ, ν_ccc, Σ₂₂, closure, buoyancy, u, v, w, T, S)

@inline ∂z_2ν_Σ₂₃(i, j, k, grid, closure, buoyancy, u, v, w, T, S) =
    2 * ∂z_aac(i, j, k, grid, ν_Σᵢⱼ, ν_cff, Σ₂₃, closure, buoyancy, u, v, w, T, S)

# At ccf
@inline ∂x_2ν_Σ₃₁(i, j, k, grid, closure, buoyancy, u, v, w, T, S) =
    2 * ∂x_caa(i, j, k, grid, ν_Σᵢⱼ, ν_fcf, Σ₃₁, closure, buoyancy, u, v, w, T, S)

@inline ∂y_2ν_Σ₃₂(i, j, k, grid, closure, buoyancy, u, v, w, T, S) =
    2 * ∂y_aca(i, j, k, grid, ν_Σᵢⱼ, ν_cff, Σ₃₂, closure, buoyancy, u, v, w, T, S)

@inline ∂z_2ν_Σ₃₃(i, j, k, grid, closure, buoyancy, u, v, w, T, S) =
    2 * ∂z_aaf(i, j, k, grid, ν_Σᵢⱼ, ν_ccc, Σ₃₃, closure, buoyancy, u, v, w, T, S)

"""
    κ_∂x_c(i, j, k, grid, c, κ, closure, buoyancy, u, v, w, T, S)

Return `κ ∂x c`, where `κ` is a function that computes
diffusivity at cell centers (location `ccc`), and `c` is an array of scalar
data located at cell centers.
"""
@inline function κ_∂x_c(i, j, k, grid, c, κ, closure, args...)
    κ = ▶x_faa(i, j, k, grid, κ, closure, args...)
    ∂x_c = ∂x_faa(i, j, k, grid, c)
    return κ * ∂x_c
end

"""
    κ_∂y_c(i, j, k, grid, c, κ, closure, buoyancy, u, v, w, T, S)

Return `κ ∂y c`, where `κ` is a function that computes
diffusivity at cell centers (location `ccc`), and `c` is an array of scalar
data located at cell centers.
"""
@inline function κ_∂y_c(i, j, k, grid, c, κ, closure, args...)
    κ = ▶y_afa(i, j, k, grid, κ, closure, args...)
    ∂y_c = ∂y_afa(i, j, k, grid, c)
    return κ * ∂y_c
end

"""
    κ_∂z_c(i, j, k, grid, c, κ, closure, buoyancy, u, v, w, T, S)

Return `κ ∂z c`, where `κ` is a function that computes
diffusivity at cell centers (location `ccc`), and `c` is an array of scalar
data located at cell centers.
"""
@inline function κ_∂z_c(i, j, k, grid, c, κ, closure, args...)
    κ = ▶z_aaf(i, j, k, grid, κ, closure, args...)
    ∂z_c = ∂z_aaf(i, j, k, grid, c)
    return κ * ∂z_c
end

"""
    ∇_κ_∇_c(i, j, k, grid, c, closure, buoyancy, u, v, w, T, S)

Return the diffusive flux divergence `∇ ⋅ (κ ∇ c)` for the turbulence
`closure`, where `c` is an array of scalar data located at cell centers.
"""
@inline ∇_κ_∇c(i, j, k, grid, c, closure::IsotropicDiffusivity, buoyancy, u, v, w, T, S) = (
      ∂x_caa(i, j, k, grid, κ_∂x_c, c, κ_ccc, closure, buoyancy, u, v, w, T, S)
    + ∂y_aca(i, j, k, grid, κ_∂y_c, c, κ_ccc, closure, buoyancy, u, v, w, T, S)
    + ∂z_aac(i, j, k, grid, κ_∂z_c, c, κ_ccc, closure, buoyancy, u, v, w, T, S)
    )

"""
    ∇_κ_∇_c(i, j, k, grid, c, closure, diffusivities)

Return the diffusive flux divergence `∇ ⋅ (κ ∇ c)` for the turbulence
`closure`, where `c` is an array of scalar data located at cell centers.
"""
@inline ∇_κ_∇c(i, j, k, grid, c, closure::IsotropicDiffusivity, diffusivities) = (
      ∂x_caa(i, j, k, grid, κ_∂x_c, c, diffusivities.κₑ, closure)
    + ∂y_aca(i, j, k, grid, κ_∂y_c, c, diffusivities.κₑ, closure)
    + ∂z_aac(i, j, k, grid, κ_∂z_c, c, diffusivities.κₑ, closure)
)

"""
    ∂ⱼ_2ν_Σ₁ⱼ(i, j, k, grid, closure, u, v, w, diffusivities)

Return the ``x``-component of the turbulent diffusive flux divergence:

`∂x(2 ν Σ₁₁) + ∂y(2 ν Σ₁₁) + ∂z(2 ν Σ₁₁)`

at the location `fcc`.
"""
@inline ∂ⱼ_2ν_Σ₁ⱼ(i, j, k, grid, closure::IsotropicDiffusivity, u, v, w, diffusivities) = (
      ∂x_2ν_Σ₁₁(i, j, k, grid, closure, u, v, w, diffusivities)
    + ∂y_2ν_Σ₁₂(i, j, k, grid, closure, u, v, w, diffusivities)
    + ∂z_2ν_Σ₁₃(i, j, k, grid, closure, u, v, w, diffusivities)
)

"""
    ∂ⱼ_2ν_Σ₂ⱼ(i, j, k, grid, closure, u, v, w, diffusivities)

Return the ``y``-component of the turbulent diffusive flux divergence:

`∂x(2 ν Σ₂₁) + ∂y(2 ν Σ₂₂) + ∂z(2 ν Σ₂₂)`

at the location `ccf`.
"""
@inline ∂ⱼ_2ν_Σ₂ⱼ(i, j, k, grid, closure::IsotropicDiffusivity, u, v, w, diffusivities) = (
      ∂x_2ν_Σ₂₁(i, j, k, grid, closure, u, v, w, diffusivities)
    + ∂y_2ν_Σ₂₂(i, j, k, grid, closure, u, v, w, diffusivities)
    + ∂z_2ν_Σ₂₃(i, j, k, grid, closure, u, v, w, diffusivities)
)

"""
    ∂ⱼ_2ν_Σ₃ⱼ(i, j, k, grid, closure, diffusivities)

Return the ``z``-component of the turbulent diffusive flux divergence:

`∂x(2 ν Σ₃₁) + ∂y(2 ν Σ₃₂) + ∂z(2 ν Σ₃₃)`

at the location `ccf`.
"""
@inline ∂ⱼ_2ν_Σ₃ⱼ(i, j, k, grid, closure::IsotropicDiffusivity, u, v, w, diffusivities) = (
      ∂x_2ν_Σ₃₁(i, j, k, grid, closure, u, v, w, diffusivities)
    + ∂y_2ν_Σ₃₂(i, j, k, grid, closure, u, v, w, diffusivities)
    + ∂z_2ν_Σ₃₃(i, j, k, grid, closure, u, v, w, diffusivities)
)


text:
#=
function fprop!(t::Convexity, g::DAT, b::VexityCache)
end
=#

ext: jl

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
res = class_homogeneous_residuals();
newt = class_homogeneous_newton();
vars = class_variables(3,2);
mu(newt, vars);

settings = class_settings();
qp, A, b, c = qp_test();

validate_dimensions(qp, vars)
initialize_newton!(newt, qp, vars, settings)

max_forks_repo_name: ohinder/Everest


text:
export solve_vf_explore!

# must precede with running the following:
# learningUpdate!(ubV, dubV, dubV_σ, uex, duex, duexσ, EV, dEV, s2idx, βΠψ, βdΠψ, ψspace, vspace, σ, β)
# learningUpdate!(ubV, uex, EV, s2idx, βΠψ, ψspace, σ, β, v, h)

function solve_vf_explore!(evs::dcdp_Emax, t::dcdp_tmpvars, p::dcdp_primitives, θt::AbstractVector, σ::Real, dograd::Bool, itype::Tuple; maxit0::Integer=40, maxit1::Integer=20, vftol::Real=1e-9)

    EV         = evs.EV
    dEV        = evs.dEV
    dEVσ       = evs.dEVσ
    wp         = p.wp
    Πz         = p.Πz
    β          = p.β

    nz,nψ,nS = size(EV)
    nSexp, dmaxp1, nd = _nSexp(wp), _dmax(wp)+1, _dmax(wp)+1

    (nz,nψ,nd) == size(t.ubVfull)         || throw(DimensionMismatch())
    (nz,nz) == size(Πz)                   || throw(DimensionMismatch())
    (nz,nψ) == size(t.lse) == size(t.tmp) || throw(DimensionMismatch())

    if dograd
        nθ = size(dEV,3)
        (nz,nψ,nθ,nS)     == size(dEV)        || throw(DimensionMismatch())
        (nz,nψ,nSexp)     == size(dEVσ)       || throw(DimensionMismatch())
        (nz,nψ,dmaxp1)    == size(t.q)        || throw(DimensionMismatch())
        (nz,nψ,nθ,nd)     == size(t.dubVfull) || throw(DimensionMismatch())
        (nz,nψ,dmaxp1)    == size(t.dubV_σ)   || throw(DimensionMismatch())
    end

    # --------- VFIt --------------

    # Views of ubV so we can efficiently access them
    @views ubV0    = t.ubVfull[:,:,  1]
    @views ubV1    = t.ubVfull[:,:,  2:dmaxp1]
    @views ubV     = t.ubVfull[:,:,  1:dmaxp1]

    @views dubV0   = t.dubVfull[:,:,:,1]
    @views dubV1   = t.dubVfull[:,:,:,2:dmaxp1]
    @views dubV    = t.dubVfull[:,:,:,1:dmaxp1]

    @views dubV_σ0 = t.dubV_σ[  :,:,  1]
    @views dubV_σ1 = t.dubV_σ[  :,:,  2:dmaxp1]
    @views dubV_σ  = t.dubV_σ[  :,:,  1:dmaxp1]

    exp2lrn = exploratory_learning(wp)
    @views βEV1   =  EV[ :,:,  exp2lrn]
    @views βdEV1  = dEV[ :,:,:,exp2lrn]
    @views βdEVσ1 = dEVσ[:,:,  exp2lrn]

    for i in ind_exp(wp)
        ip = sprime(wp,i,0)
        horzn = _horizon(wp,i)

        @views EV0 = EV[:,:,ip]
        @views dEV0 = dEV[:,:,:,ip]
        @views dEVσ1 = dEVσ[ :,:,ip]

        # compute u + βEV(d) ∀ d ∈ actionspace(wp,i)
        fillflows!(ubV, flow, p, θt, σ, i, itype...)
        ubV0 .+= β .* EV0
        ubV1 .+= βEV1 # β already baked in

        if dograd
            fillflows_grad!(dubV, flowdθ, p, θt, σ, i, itype...)
            fillflows!(   dubV_σ, flowdσ, p, θt, σ, i, itype...)
            dubV0   .+= β .* dEV0
            dubV1   .+= βdEV1  # β already baked in
            dubV_σ0 .+= β .* dEVσ1
            dubV_σ1 .+= βdEVσ1 # β already baked in
        end

        if horzn == :Finite
            if dograd
                @views vfit!(EV[:,:,i], dEV[:,:,:,i], dEVσ[:,:,i], t, p)
            else
                @views vfit!(EV[:,:,i], t, p)
            end

        elseif horzn == :Infinite
            converged, iter, bnds =  solve_inf_vfit_pfit!(EV0, t, p; vftol=vftol, maxit0=maxit0, maxit1=maxit1)
            converged || @warn "Did not converge at state $i after $iter pfit. McQueen-Porteus bnds: $bnds. θt = $θt, σ = $σ"

            if dograd
                # TODO: only allows 0-payoff if no action
                ubV[:,:,1] .= β .* EV0
                gradinf!(dEV0, dEVσ1, t, p, true)
            end
        else
            throw(error("i = $i, horzn = $horzn but must be :Finite or :Infinite"))
        end
    end
end


text:
function Alhistory(; sessionprint = false)
   !sessionprint && (myhistory = activehistory)
   sessionprint  && (myhistory = sessionhistory)

   currentcommand =
     [priorcommands(myhistory)..., "<< present >>",
     futurecommands(myhistory)...]
   currentfeasible =
     [nfeasible.(priorlogicsets(myhistory))..., "...",
     nfeasible.(futurelogicsets(myhistory))...]

   txtout = "| Command | # feasible |\n| --- | --- |\n"
   for i in 1:length(currentcommand)
       txtout *= ("| $(markdownescape(currentcommand[i])) | $(currentfeasible[i]) |\n")
   end
   printmarkdown(txtout)
end


text:
function σx(n::Int,spinor::Spinor)
    newspinor = deepcopy(spinor)
    for state in newspinor.spins
        state[n] = (state[n] + 1) % 2
    end
    return newspinor
end

function σy(n::Int,spinor::Spinor)
    newspinor = deepcopy(spinor)
    for (i,state) in enumerate(newspinor.spins)
        state[n] = (state[n] + 1) % 2
        newspinor.coefficients[i] *= spinor.spins[i][n] == 1 ? im : -im
    end
    return newspinor
end

function σz(n::Int,spinor::Spinor)
    newspinor = deepcopy(spinor)
    for i in 1:length(newspinor.coefficients)
        newspinor.coefficients[i] *= ((spinor.spins[i][n] == 1) ? 1 : -1)
    end
    return newspinor
end

σiσj(i,j,spinor::Spinor) = σx(i,σx(j, spinor)) + σy(i,σy(j, spinor)) + σz(i,σz(j, spinor))

dipoledipole(i,j,spinor::Spinor) = 1/abs(i-j)^3 * ((-1)*σx(i,σx(j, spinor)) + (-1)*σy(i,σy(j, spinor)) + 2*σz(i,σz(j, spinor)))


text:
export absolute_primitive_element, absolute_simple_field

################################################################################
#
#  Absolute primitive element
#
################################################################################

@doc Markdown.doc"""
    absolute_primitive_element(K::NumField) -> NumFieldElem
Given a number field $K$, this function returns an element $\gamma \in K$
such that $K = \mathbf{Q}(\gamma)$.
"""
absolute_primitive_element(::NumField)

function absolute_primitive_element(K::NfRelNS)
  k = base_field(K)
  a = primitive_element(K)
  if is_absolutely_primitive(a)
    return a
  end
  gk = absolute_primitive_element(k)
  a += gk
  while !is_absolutely_primitive(a)
    a += gk
  end
  return a
end

function absolute_primitive_element(K::NfRel)
  k = base_field(K)
  gk = absolute_primitive_element(k)
  a = gen(K)
  while !is_absolutely_primitive(a)
    a += gk
  end
  return a
end

function absolute_primitive_element(K::NfAbsNS)
  return primitive_element(K)
end

function absolute_primitive_element(K::AnticNumberField)
  return gen(K)
end

################################################################################
#
#  Is absolutely primitive test
#
################################################################################

#TODO: Put some more thought in it.
#In practice, we should check this modularly
function is_absolutely_primitive(a::NumFieldElem)
  c = conjugates_arb(a, 16)
  is_primitive = true
  for i = 2:length(c)
    for j = i+1:length(c)
      if overlaps(c[i], c[i])
        is_primitive = false
        break
      end
    end
    if !is_primitive
      break
    end
  end
  if is_primitive
    return true
  end
  return degree(absolute_minpoly(a)) == absolute_degree(parent(a))
end

function is_absolutely_primitive(a::T) where T <: Union{NfRelNSElem{nf_elem}, NfRelElem{nf_elem}}
  L = parent(a)
  rt, PR, tmp = _setup_block_system(L)
  if _is_primitive_via_block(a, rt, PR, tmp)
    return true
  end
  return degree(absolute_minpoly(a)) == absolute_degree(parent(a))
end

################################################################################
#
#  Absolute field
#
################################################################################

@doc Markdown.doc"""
    absolute_simple_field(K::NumField) -> NumField, Map

Given a number field $K$, this function returns an absolute simple number field
$M/\mathbf{Q}$ together with a $\mathbf{Q}$-linear isomorphism $M \to K$.
"""
absolute_simple_field(::NumField)

function absolute_simple_field(K::AnticNumberField)
  return K, id_hom(K)
end

function absolute_simple_field(K::NfAbsNS; cached::Bool = true, simplify::Bool = false)
  return simple_extension(K, cached = cached, simplified = simplify)
end

function absolute_simple_field(K::NumField; cached::Bool = false, simplify::Bool = false)
  if simplify
    return simplified_absolute_field(K, cached = cached)
  end
  el = absolute_primitive_element(K)
  f = absolute_minpoly(el)
  Kabs, gKabs = number_field(f, cached = cached, check = false)
  mp = hom(Kabs, K, el)
  embed(mp)
  return Kabs, mp
end

#Special function for NfRel{nf_elem}. In this case, we can easily construct the
#inverse of the isomorphism, so we do it separately
function absolute_simple_field(K::NfRel{nf_elem}; cached::Bool = false, simplify::Bool = false)
  if simplify
    return simplified_absolute_field(K, cached = cached)
  end
  Ka, a, b, c = _absolute_field(K, cached = cached)
  h1 = hom(Ka, K, c, inverse = (a, b))
  embed(h1)
  embed(MapFromFunc(x->preimage(h1, x), K, Ka))
  return Ka, h1
end

#Trager: p4, Algebraic Factoring and Rational Function Integration
function _absolute_field(K::NfRel; cached::Bool = false)
  f = K.pol
  kx = parent(f)
  k = base_ring(kx)
  Qx = parent(k.pol)

  l = 0
  #TODO: Until we don't find a primitive element, we should check
  #primitiveness instead of computing the full norm.
  g = f
  N = norm(g)

  while true
    @assert degree(N) == degree(g) * degree(k)
    if !is_constant(N) && is_squarefree(N)
      break
    end
    l += 1
    g = compose(f, gen(kx) - l*gen(k))
    N = norm(g)
  end

  Ka, gKa = NumberField(N, "x", cached = cached, check = false)
  KaT, T = PolynomialRing(Ka, "T", cached = false)

  # map Ka -> K: gen(Ka) -> gen(K)+ k gen(k)

  # gen(k) -> Root(gcd(g, poly(k)))  #gcd should be linear:
  # g in kx = (Q[a])[x]. Want to map x -> gen(Ka), a -> T

  gg = zero(KaT)
  for i=degree(g):-1:0
    auxp = change_base_ring(Ka, Qx(coeff(g, i)), parent = KaT)
    gg = gg*gKa
    add!(gg, gg, auxp)
    #gg = gg*gKa + auxp
  end

  q = gcd(gg, change_base_ring(Ka, k.pol, parent = KaT))
  @assert degree(q) == 1
  al = -constant_coefficient(q)//leading_coefficient(q)
  be = gKa - l*al
  ga = gen(K) + l*gen(k)

  #al -> gen(k) in Ka
  #be -> gen(K) in Ka
  #ga -> gen(Ka) in K
  return Ka, al, be, ga
end

################################################################################
#
#  Collapse_top_layer function
#
################################################################################

function collapse_top_layer(K::NfRel{T}; cached::Bool = false) where T <: SimpleNumFieldElem
  Ka, a, b, c = _absolute_field(K)
  h1 = hom(Ka, K, c, inverse = (a, b))
  h2 = hom(base_field(K), Ka, a, check = false)
  embed(h1)
  embed(MapFromFunc(x->preimage(h1, x), K, Ka))
  embed(h2)
  return Ka, h1, h2
end


text:
# Autogenerated wrapper script for OpenFAST_jll for x86_64-apple-darwin-libgfortran5
export aerodyn_driver, beamdyn_driver, dwm_driver_wind_farm, feam_driver, hydrodyn_driver, inflowwind_driver, libaeroacoustics, libaerodyn14lib, libaerodynlib, libafinfolib, libawaelib, libbeamdynlib, libelastodynlib, libextptfm_mckflib, libfeamlib, libfoamfastlib, libfvwlib, libhydrodynlib, libicedynlib, libicefloelib, libifwlib, libmapcpplib, libmaplib, libmoordynlib, libnwtclibs, libopenfast_postlib, libopenfast_prelib, libopenfastlib, libopenfoamtypeslib, liborcaflexlib, libscdataexlib, libscdataexttypeslib, libscfastlib, libsctypeslib, libservodynlib, libsubdynlib, libuaaerolib, libversioninfolib, libwdlib, moordyn_driver, openfast, orca_driver, servodyn_driver, subdyn_driver, turbsim, unsteadyaero_driver

using OpenBLAS32_jll
using LAPACK_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("OpenFAST")
JLLWrappers.@declare_library_product(libaeroacoustics, "@rpath/libaeroacoustics.dylib")
JLLWrappers.@declare_library_product(libaerodyn14lib, "@rpath/libaerodyn14lib.dylib")
JLLWrappers.@declare_library_product(libaerodynlib, "@rpath/libaerodynlib.dylib")
JLLWrappers.@declare_library_product(libafinfolib, "@rpath/libafinfolib.dylib")
JLLWrappers.@declare_library_product(libawaelib, "@rpath/libawaelib.dylib")
JLLWrappers.@declare_library_product(libbeamdynlib, "@rpath/libbeamdynlib.dylib")
JLLWrappers.@declare_library_product(libelastodynlib, "@rpath/libelastodynlib.dylib")
JLLWrappers.@declare_library_product(libextptfm_mckflib, "@rpath/libextptfm_mckflib.dylib")
JLLWrappers.@declare_library_product(libfeamlib, "@rpath/libfeamlib.dylib")
JLLWrappers.@declare_library_product(libfoamfastlib, "@rpath/libfoamfastlib.dylib")
JLLWrappers.@declare_library_product(libfvwlib, "@rpath/libfvwlib.dylib")
JLLWrappers.@declare_library_product(libhydrodynlib, "@rpath/libhydrodynlib.dylib")
JLLWrappers.@declare_library_product(libicedynlib, "@rpath/libicedynlib.dylib")
JLLWrappers.@declare_library_product(libicefloelib, "@rpath/libicefloelib.dylib")
JLLWrappers.@declare_library_product(libifwlib, "@rpath/libifwlib.dylib")
JLLWrappers.@declare_library_product(libmapcpplib, "@rpath/libmapcpplib.dylib")
JLLWrappers.@declare_library_product(libmaplib, "@rpath/libmaplib.dylib")
JLLWrappers.@declare_library_product(libmoordynlib, "@rpath/libmoordynlib.dylib")
JLLWrappers.@declare_library_product(libnwtclibs, "@rpath/libnwtclibs.dylib")
JLLWrappers.@declare_library_product(libopenfast_postlib, "@rpath/libopenfast_postlib.dylib")
JLLWrappers.@declare_library_product(libopenfast_prelib, "@rpath/libopenfast_prelib.dylib")
JLLWrappers.@declare_library_product(libopenfastlib, "@rpath/libopenfastlib.dylib")
JLLWrappers.@declare_library_product(libopenfoamtypeslib, "@rpath/libopenfoamtypeslib.dylib")
JLLWrappers.@declare_library_product(liborcaflexlib, "@rpath/liborcaflexlib.dylib")
JLLWrappers.@declare_library_product(libscdataexlib, "@rpath/libscdataexlib.dylib")
JLLWrappers.@declare_library_product(libscdataexttypeslib, "@rpath/libscdataextypeslib.dylib")
JLLWrappers.@declare_library_product(libscfastlib, "@rpath/libscfastlib.dylib")
JLLWrappers.@declare_library_product(libsctypeslib, "@rpath/libsctypeslib.dylib")
JLLWrappers.@declare_library_product(libservodynlib, "@rpath/libservodynlib.dylib")
JLLWrappers.@declare_library_product(libsubdynlib, "@rpath/libsubdynlib.dylib")
JLLWrappers.@declare_library_product(libuaaerolib, "@rpath/libuaaerolib.dylib")
JLLWrappers.@declare_library_product(libversioninfolib, "@rpath/libversioninfolib.dylib")
JLLWrappers.@declare_library_product(libwdlib, "@rpath/libwdlib.dylib")
JLLWrappers.@declare_executable_product(aerodyn_driver)
JLLWrappers.@declare_executable_product(beamdyn_driver)
JLLWrappers.@declare_executable_product(dwm_driver_wind_farm)
JLLWrappers.@declare_executable_product(feam_driver)
JLLWrappers.@declare_executable_product(hydrodyn_driver)
JLLWrappers.@declare_executable_product(inflowwind_driver)
JLLWrappers.@declare_executable_product(moordyn_driver)
JLLWrappers.@declare_executable_product(openfast)
JLLWrappers.@declare_executable_product(orca_driver)
JLLWrappers.@declare_executable_product(servodyn_driver)
JLLWrappers.@declare_executable_product(subdyn_driver)
JLLWrappers.@declare_executable_product(turbsim)
JLLWrappers.@declare_executable_product(unsteadyaero_driver)
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS32_jll, LAPACK_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libaeroacoustics,
        "lib/libaeroacoustics.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libaerodyn14lib,
        "lib/libaerodyn14lib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libaerodynlib,
        "lib/libaerodynlib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libafinfolib,
        "lib/libafinfolib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libawaelib,
        "lib/libawaelib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libbeamdynlib,
        "lib/libbeamdynlib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libelastodynlib,
        "lib/libelastodynlib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libextptfm_mckflib,
        "lib/libextptfm_mckflib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libfeamlib,
        "lib/libfeamlib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libfoamfastlib,
        "lib/libfoamfastlib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libfvwlib,
        "lib/libfvwlib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libhydrodynlib,
        "lib/libhydrodynlib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libicedynlib,
        "lib/libicedynlib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libicefloelib,
        "lib/libicefloelib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libifwlib,
        "lib/libifwlib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmapcpplib,
        "lib/libmapcpplib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmaplib,
        "lib/libmaplib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmoordynlib,
        "lib/libmoordynlib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnwtclibs,
        "lib/libnwtclibs.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopenfast_postlib,
        "lib/libopenfast_postlib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopenfast_prelib,
        "lib/libopenfast_prelib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopenfastlib,
        "lib/libopenfastlib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopenfoamtypeslib,
        "lib/libopenfoamtypeslib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        liborcaflexlib,
        "lib/liborcaflexlib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libscdataexlib,
        "lib/libscdataexlib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libscdataexttypeslib,
        "lib/libscdataextypeslib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libscfastlib,
        "lib/libscfastlib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsctypeslib,
        "lib/libsctypeslib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libservodynlib,
        "lib/libservodynlib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsubdynlib,
        "lib/libsubdynlib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libuaaerolib,
        "lib/libuaaerolib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libversioninfolib,
        "lib/libversioninfolib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libwdlib,
        "lib/libwdlib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        aerodyn_driver,
        "bin/aerodyn_driver",
    )

    JLLWrappers.@init_executable_product(
        beamdyn_driver,
        "bin/beamdyn_driver",
    )

    JLLWrappers.@init_executable_product(
        dwm_driver_wind_farm,
        "bin/dwm_driver_wind_farm",
    )

    JLLWrappers.@init_executable_product(
        feam_driver,
        "bin/feam_driver",
    )

    JLLWrappers.@init_executable_product(
        hydrodyn_driver,
        "bin/hydrodyn_driver",
    )

    JLLWrappers.@init_executable_product(
        inflowwind_driver,
        "bin/inflowwind_driver",
    )

    JLLWrappers.@init_executable_product(
        moordyn_driver,
        "bin/moordyn_driver",
    )

    JLLWrappers.@init_executable_product(
        openfast,
        "bin/openfast",
    )

    JLLWrappers.@init_executable_product(
        orca_driver,
        "bin/orca_driver",
    )

    JLLWrappers.@init_executable_product(
        servodyn_driver,
        "bin/servodyn_driver",
    )

    JLLWrappers.@init_executable_product(
        subdyn_driver,
        "bin/subdyn_driver",
    )

    JLLWrappers.@init_executable_product(
        turbsim,
        "bin/turbsim",
    )

    JLLWrappers.@init_executable_product(
        unsteadyaero_driver,
        "bin/unsteadyaero_driver",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#######################################################################################################################################################################################################
#
# Changes to the function
# General
#     2022-Feb-03: add recursive variable query feature
#
#######################################################################################################################################################################################################
"""

    find_variable(ds::Dataset, var_name::String)

Return the path to dataset if it exists, given
- `ds` NCDatasets.Dataset type dataset
- `var_name` Variable to read
"""
function find_variable(ds::Dataset, var_name::String)
    # if var_name is in the current dataset, return it
    if var_name in keys(ds)
        return ds[var_name]
    end;

    # loop through the groups and find the data
    _dvar = nothing;
    for _group in keys(ds.group)
        _dvar = find_variable(ds.group[_group], var_name)
        if _dvar !== nothing
            break;
        end;
    end;

    # return the variable if it exists, otherwise return nothing
    return _dvar
end

#######################################################################################################################################################################################################
#
# Changes to the function
# General
#     2021-Dec-24: migrate the function from PkgUtility to NetcdfIO
#     2022-Jan-28: fix documentation
#
#######################################################################################################################################################################################################
"""
NCDatasets.jl and NetCDF.jl both provide function to read data out from NC dataset. However, while NetCDF.jl is more convenient to use (less lines of code to read data), NCDatasets.jl is better to
    read a subset from the dataset and is able to detect the scale factor and offset. Here, we used a wrapper function to read NC dataset using NCDatasets.jl:

$(METHODLIST)

"""
function read_nc end

#######################################################################################################################################################################################################
#
# Changes to the function
# General
#     2021-Dec-24: migrate the function from PkgUtility to NetcdfIO
#     2022-Jan-28: fix documentation
#     2022-Feb-03: add recursive variable query feature
#     2022-Feb-03: add option to read raw data to avoid NCDatasets transform errors
# Bug fixes
#     2021-Dec-24: fix the bug that reads integer as float (e.g., ind)
#
#######################################################################################################################################################################################################
"""
When only file name and variable label are provided, `read_nc` function reads out all the data:

    read_nc(file::String, var_name::String; transform::Bool = true)

Read data from NC file, given
- `file` Path of the netcdf dataset
- `var_name` Variable to read
- `transform` If true, transform the data using NCDatasets rules, otherwise read the raw data

Note that the missing data will be labeled as NaN.

---
# Examples
```julia
# read data labeled as test from test.nc
save_nc!("test.nc", "test", rand(36,18,12), Dict("description" => "Random randoms"));
data = read_nc("test.nc", "test");
```
"""
read_nc(file::String, var_name::String; transform::Bool = true) = (
    _dset = Dataset(file, "r");

    _fvar = find_variable(_dset, var_name);
    if _fvar === nothing
        close(_dset)
        @error "$(var_name) does not exist in $(file)!";
    end;

    if transform
        _dvar = _fvar[:,:];
    else
        _dvar = _fvar.var[:,:];
    end;
    close(_dset);

    if sum(ismissing.(_dvar)) == 0
        return _dvar
    end;

    return replace(_dvar, missing=>NaN)
);

#######################################################################################################################################################################################################
#
# Changes to the function
# General
#     2021-Dec-24: migrate the function from PkgUtility to NetcdfIO
#     2022-Jan-28: fix documentation
#     2022-Feb-03: add option to read raw data to avoid NCDatasets transform errors
#
#######################################################################################################################################################################################################
"""
If a float type is given, the data will be converted to T, namely the output will be an array of T type numbers:

    read_nc(T, file::String, var_name::String; transform::Bool = true)

Read data from nc file, given
- `T` Number type
- `file` Path of the netcdf dataset
- `var_name` Variable name
- `transform` If true, transform the data using NCDatasets rules, otherwise read the raw data

---
# Examples
```julia
# read data labeled as test from test.nc as Float32
save_nc!("test.nc", "test", rand(36,18,12), Dict("description" => "Random randoms"));
data = read_nc(Float32, "test.nc", "test");
```
"""
read_nc(T, file::String, var_name::String; transform::Bool = true) = T.(read_nc(file, var_name; transform = transform));

#######################################################################################################################################################################################################
#
# Changes to the function
# General
#     2021-Dec-24: migrate the function from PkgUtility to NetcdfIO
#     2022-Jan-28: fix documentation
#     2022-Feb-03: add recursive variable query feature
#     2022-Feb-03: add option to read raw data to avoid NCDatasets transform errors
#     2022-Feb-04: allow to read value from 1D array as well
# Bug fixes
#     2021-Dec-24: fix the bug that reads integer as float (e.g., ind)
#     2022-Jan-20: add dimension control to avoid errors
#
#######################################################################################################################################################################################################
"""
In many cases, the NC dataset can be very huge, and reading all the data points into one array could be time and memory consuming. In this case, reading a subset of data would be the best option:

    read_nc(file::String, var_name::String, indz::Int; transform::Bool = true)

Read a subset from nc file, given
- `file` Path of the netcdf dataset
- `var_name` Variable name
- `indz` The 3rd index of subset data to read
- `transform` If true, transform the data using NCDatasets rules, otherwise read the raw data

Note that the dataset must be a 3D array to use this method.

---
# Examples
```julia
# read 1st layer data labeled as test from test.nc
save_nc!("test.nc", "test", rand(36,18,12), Dict("description" => "Random randoms"));
data = read_nc("test.nc", "test", 1);
```
"""
read_nc(file::String, var_name::String, indz::Int; transform::Bool = true) = (
    _ndim = size_nc(file, var_name)[1];
    @assert _ndim in [1,3] "The dataset must be a 1D or 3D array to use this method!";

    _dset = Dataset(file, "r");
    _fvar = find_variable(_dset, var_name);
    if transform
        _dvar = (_ndim == 1 ? _fvar[indz] : _fvar[:,:,indz]);
    else
        _dvar = (_ndim == 1 ? _fvar.var[indz] : _fvar.var[:,:,indz]);
    end;
    close(_dset);

    if sum(ismissing.(_dvar)) == 0
        return _dvar
    end;

    return replace(_dvar, missing=>NaN)
);

#######################################################################################################################################################################################################
#
# Changes to the function
# General
#     2021-Dec-24: migrate the function from PkgUtility to NetcdfIO
#     2022-Jan-28: fix documentation
#     2022-Feb-03: add option to read raw data to avoid NCDatasets transform errors
#
#######################################################################################################################################################################################################
"""
Similarly, one may want to read the subset as a certain type using

    read_nc(T, file::String, var_name::String, indz::Int; transform::Bool = true)

Read a subset from nc file, given
- `T` Number type
- `file` Path of the netcdf dataset
- `var_name` Variable name
- `indz` The 3rd index of subset data to read
- `transform` If true, transform the data using NCDatasets rules, otherwise read the raw data

---
# Examples
```julia
# read 1st layer data labeled as test from test.nc as Float32
save_nc!("test.nc", "test", rand(36,18,12), Dict("description" => "Random randoms"));
data = read_nc(Float32, "test.nc", "test", 1);
```
"""
read_nc(T, file::String, var_name::String, indz::Int; transform::Bool = true) = T.(read_nc(file, var_name, indz; transform = transform));

#######################################################################################################################################################################################################
#
# Changes to the function
# General
#     2021-Dec-24: migrate the function from PkgUtility to NetcdfIO
#     2022-Jan-28: fix documentation
#     2022-Feb-03: add recursive variable query feature
#     2022-Feb-03: add option to read raw data to avoid NCDatasets transform errors
# Bug fixes
#     2021-Dec-24: fix the bug that reads integer as float (e.g., ind)
#     2022-Jan-20: add dimension control to avoid errors
#
#######################################################################################################################################################################################################
"""
Another convenient wrapper is to read all the data for given index in x and y, for example, if one wants to read the time series of data at a given site:

    read_nc(file::String, var_name::String, indx::Int, indy::Int; transform::Bool = true)

Read the time series of data for a site, given
- `file` Path of the netcdf dataset
- `var_name` Variable name
- `indx` The 1st index of subset data to read, typically longitude
- `indy` The 2nd index of subset data to read, typically latitude
- `transform` If true, transform the data using NCDatasets rules, otherwise read the raw data

---
# Examples
```julia
save_nc!("test1.nc", "test", rand(36,18), Dict("description" => "Random randoms"));
save_nc!("test2.nc", "test", rand(36,18,12), Dict("description" => "Random randoms"));
data1 = read_nc("test1.nc", "test", 1, 1);
data2 = read_nc("test2.nc", "test", 1, 1);
```
"""
read_nc(file::String, var_name::String, indx::Int, indy::Int; transform::Bool = true) = (
    _ndim = size_nc(file, var_name)[1];
    @assert 2 <= _ndim <= 3 "The dataset must be a 2D or 3D array to use this method!";

    _dset = Dataset(file, "r");
    _fvar = find_variable(_dset, var_name);
    if transform
        _dvar = (_ndim==2 ? _fvar[indx,indy] : _fvar[indx,indy,:]);
    else
        _dvar = (_ndim==2 ? _fvar.var[indx,indy] : _fvar.var[indx,indy,:]);
    end;
    close(_dset);

    if sum(ismissing.(_dvar)) == 0
        return _dvar
    end;

    return replace(_dvar, missing=>NaN)
);

#######################################################################################################################################################################################################
#
# Changes to the function
# General
#     2021-Dec-24: migrate the function from PkgUtility to NetcdfIO
#     2022-Jan-28: fix documentation
#     2022-Feb-03: add option to read raw data to avoid NCDatasets transform errors
#
#######################################################################################################################################################################################################
"""
Similarly, one may want to read the subset as a certain type using

    read_nc(T, file::String, var_name::String, indx::Int, indy::Int; transform::Bool = true)

Read the time series of data for a site, given
- `T` Number type
- `file` Path of the netcdf dataset
- `var_name` Variable name
- `indx` The 1st index of subset data to read, typically longitude
- `indy` The 2nd index of subset data to read, typically latitude
- `transform` If true, transform the data using NCDatasets rules, otherwise read the raw data

---
# Examples
```julia
save_nc!("test1.nc", "test", rand(36,18), Dict("description" => "Random randoms"));
save_nc!("test2.nc", "test", rand(36,18,12), Dict("description" => "Random randoms"));
data1 = read_nc(Float32, "test1.nc", "test", 1, 1);
data2 = read_nc(Float32, "test2.nc", "test", 1, 1);
```
"""
read_nc(T, file::String, var_name::String, indx::Int, indy::Int; transform::Bool = true) = T.(read_nc(file, var_name, indx, indy; transform = transform));

#######################################################################################################################################################################################################
#
# Changes to the function
# General
#     2021-Dec-24: migrate the function from PkgUtility to NetcdfIO
#     2022-Jan-28: fix documentation
#     2022-Feb-03: add recursive variable query feature
#     2022-Feb-03: add option to read raw data to avoid NCDatasets transform errors
# Bug fixes
#     2021-Dec-24: fix the bug that reads integer as float (e.g., ind)
#     2022-Jan-20: add dimension control to avoid errors
#
#######################################################################################################################################################################################################
"""
Another convenient wrapper is to read the data for given index in x, y, and z, for example, if one wants to read the time series of data at a given site:

    read_nc(file::String, var_name::String, indx::Int, indy::Int, indz::Int; transform::Bool = true)

Read the time series of data for a site, given
- `file` Path of the netcdf dataset
- `var_name` Variable name
- `indx` The 1st index of subset data to read, typically longitude
- `indy` The 2nd index of subset data to read, typically latitude
- `indz` The 3rd index of subset data to read, typically time
- `transform` If true, transform the data using NCDatasets rules, otherwise read the raw data

---
# Examples
```julia
save_nc!("test.nc", "test", rand(36,18,12), Dict("description" => "Random randoms"));
data = read_nc("test.nc", "test", 1, 1, 1);
```
"""
read_nc(file::String, var_name::String, indx::Int, indy::Int, indz::Int; transform::Bool = true) = (
    @assert size_nc(file, var_name)[1] == 3 "The dataset must be a 3D array to use this method!";

    _dset = Dataset(file, "r");
    _fvar = find_variable(_dset, var_name);
    if transform
        _dvar = _fvar[indx,indy,indz];
    else
        _dvar = _fvar.var[indx,indy,indz];
    end;
    close(_dset);

    return ismissing(_dvar) ? NaN : _dvar
);

#######################################################################################################################################################################################################
#
# Changes to the function
# General
#     2021-Dec-24: migrate the function from PkgUtility to NetcdfIO
#     2022-Jan-28: fix documentation
#     2022-Feb-03: add option to read raw data to avoid NCDatasets transform errors
#
#######################################################################################################################################################################################################
"""
Similarly, one may want to read the data as a certain type using

    read_nc(T, file::String, var_name::String, indx::Int, indy::Int, indz::Int; transform::Bool = true)

Read the time series of data for a site, given
- `T` Number type
- `file` Path of the netcdf dataset
- `var_name` Variable name
- `indx` The 1st index of subset data to read, typically longitude
- `indy` The 2nd index of subset data to read, typically latitude
- `indz` The 3rd index of subset data to read, typically time
- `transform` If true, transform the data using NCDatasets rules, otherwise read the raw data

---
# Examples
```julia
save_nc!("test.nc", "test", rand(36,18,12), Dict("description" => "Random randoms"));
data = read_nc(Float32, "test.nc", "test", 1, 1, 1);
```
"""
read_nc(T, file::String, var_name::String, indx::Int, indy::Int, indz::Int; transform::Bool = true) = T.(read_nc(file, var_name, indx, indy, indz; transform = transform));

#######################################################################################################################################################################################################
#
# Changes to the function
# General
#     2022-Jan-04: add this function to read netcdf as DataFrame
#     2022-Jan-28: fix documentation
#     2022-Feb-03: add option to read raw data to avoid NCDatasets transform errors
# Bug fixes
#     2021-Dec-24: fix the bug that reads integer as float (e.g., ind)
#
#######################################################################################################################################################################################################
"""
The method below reads all the 1D data (with the same length) into a DataFrame

    read_nc(file::String, selections::Vector{String} = varname_nc(file); transform::Bool = true)

Read the selected variables from a netcdf file as a DataFrame, given
- `file` Path of the netcdf dataset
- `selections` Variables to read from the file
- `transform` If true, transform the data using NCDatasets rules, otherwise read the raw data

---
# Examples
```julia
df_raw = DataFrame();
df_raw[!,"A"] = rand(5);
df_raw[!,"B"] = rand(5);
df_raw[!,"C"] = rand(5);
save_nc!("test.nc", df_raw);
df_new = read_nc("test.nc");
df_new = read_nc("test.nc", ["A", "B"]);
```
"""
read_nc(file::String, selections::Vector{String} = varname_nc(file); transform::Bool = true) = (
    _dims = [size_nc(file, _var)[1] for _var in selections];
    _lens = [size_nc(file, _var)[2][1] for _var in selections];
    @assert all(_dims .== 1) "All variables need to be 1D!";
    @assert all(_lens .== _lens[1]) "Dimensions of the variables need to be the same!";

    return DataFrame( [Pair(_var, read_nc(file, _var; transform = transform)) for _var in selections] )
);


text:
export SpfPosteriorMcmcOutput, SpfPosteriorModeConditionalMcmcOutput, SpfPosteriorDispersionConditionalMcmcOutput
export SisPosteriorMcmcOutput, SisPosteriorModeConditionalMcmcOutput, SisPosteriorDispersionConditionalMcmcOutput
export SimPosteriorModeConditionalMcmcOutput, SimPosteriorDispersionConditionalMcmcOutput

# ==========
#    SPF
# ==========

struct SpfPosteriorMcmcOutput{T<:Union{Int, String}}
    I_sample::Vector{Path{T}}
    γ_sample::Vector{Float64}
    log_post::Dict # Had to do dict since might want different output when (Iᵐ, γ) updated jointly
    dist::Metric
    I_prior::SPF{T}
    γ_prior::ContinuousUnivariateDistribution
    data::Vector{Path{T}}
    performance_measures::Dict
end

struct SpfPosteriorModeConditionalMcmcOutput{T<:Union{Int, String}}
    γ_fixed::Float64
    I_sample::Vector{Path{T}}
    dist::Metric
    I_prior::SPF{T}
    data::Vector{Path{T}}
    performance_measures::Dict
end

struct SpfPosteriorDispersionConditionalMcmcOutput{T<:Union{Int, String}}
    I_fixed::Path{T}
    γ_sample::Vector{Float64}
    γ_prior::ContinuousUnivariateDistribution
    data::Vector{Path{T}}
    performance_measures::Dict
end

function Base.show(io::IO, output::T) where {T<:SpfPosteriorMcmcOutput}
    title = "MCMC Sample for SPF Posterior"
    n = length(title)
    println(io, title)
    println(io, "-"^n)
    for (key, value) in output.performance_measures
        println(io, key, ": ", value)
    end
end

function Base.show(io::IO, output::T) where {T<:SpfPosteriorModeConditionalMcmcOutput}
    title = "MCMC Sample for SPF Posterior (Mode Conditional)"
    n = length(title)
    println(io, title)
    println(io, "-"^n)
    for (key, value) in output.performance_measures
        println(io, key, ": ", value)
    end
end

function Base.show(io::IO, output::T) where {T<:SpfPosteriorDispersionConditionalMcmcOutput}
    title = "MCMC Sample for SPF Posterior (Dispersion Conditional)"
    n = length(title)
    println(io, title)
    println(io, "-"^n)
    for (key, value) in output.performance_measures
        println(io, key, ": ", value)
    end
end

struct SisPosteriorMcmcOutput
    S_sample::Vector{InteractionSequence{Int}}
    γ_sample::Vector{Float64}
    posterior::SisPosterior
    suff_stat_trace::Vector{Float64}
    performace_measures::Dict
end

struct SisPosteriorModeConditionalMcmcOutput
    γ_fixed::Float64
    S_sample::Vector{Vector{Path{Int}}}
    posterior::SisPosterior
    suff_stat_trace::Vector{Float64}
    performance_measures::Dict
end

struct SisPosteriorDispersionConditionalMcmcOutput
    S_fixed::Vector{Path{Int}}
    γ_sample::Vector{Float64}
    posterior::SisPosterior
    performance_measures::Dict
end

# ==========
#    SIM
# ==========

struct SimPosteriorMcmcOutput
    S_sample::Vector{InteractionSequence{Int}}
    γ_sample::Vector{Float64}
    posterior::SimPosterior
    suff_stat_trace::Vector{Float64}
    performace_measures::Dict
end

struct SimPosteriorModeConditionalMcmcOutput
    γ_fixed::Float64
    S_sample::Vector{Vector{Path{Int}}}
    posterior::SimPosterior
    suff_stat_trace::Vector{Float64}
    performance_measures::Dict
end

struct SimPosteriorDispersionConditionalMcmcOutput
    S_fixed::Vector{Path{Int}}
    γ_sample::Vector{Float64}
    γ_prior::ContinuousUnivariateDistribution
    data::Vector{Vector{Path{Int}}}
    performance_measures::Dict
end


text:
abstract type AbstractSym{TAG} end

const DEFAULT_TAG = Any

mutable struct Sym{TAG} <: AbstractSym{TAG}
    head::Symbol
    args::Vector{Any}

    function Sym{TAG}(head::Symbol, @nospecialize args...) where {TAG}
        return new{TAG}(head, collect(args))
    end

    function Sym{TAG}(::Val{:noargs}, head::Symbol, @nospecialize args...) where {TAG}
        return new{TAG}(head, collect(args))
    end
end

# Constructors

# Syms
Sym(x::Sym) = x
Sym{TAG}(x::Sym{TAG}) where {TAG} = x
Sym{TAG}(x::Sym) where {TAG} = Sym{TAG}(Val(:noargs), gethead(x), getargs(x)...)

# Generic
Sym(x) = Sym{tagof(x)}(x)
Sym{TAG}(x) where {TAG} = Sym{TAG}(:object, x)
Sym{TAG}(x::Type) where {TAG} = Sym{TAG}(:type, x)
Sym{TAG}(x::Symbol) where {TAG} = Sym{TAG}(:symbol, x)
Sym{TAG}(x::Function) where {TAG} = Sym{TAG}(:fn, x)

# Expressions
Sym(x::Expr) = Sym(Val(x.head), x)
Sym(::Val{head}, x::Expr) where {head} = Sym{tagof(x)}(Val(head), x)
Sym{TAG}(x::Expr) where {TAG} = convert(Sym{TAG}, Sym(x))

function Sym{TAG}(::Val{head}, x::Expr) where {TAG,head}
    @assert(x.head === head)
    args = map(Sym, filter(!symignore, x.args))
    return Sym{TAG}(Val(:noargs), head, args...)
end

function Sym(::Val{:call}, x::Expr)
    @assert x.head === :call
    args = map(Sym, filter(!symignore, x.args))
    f, args = args[1], args[2:end]
    return convert(Sym, f(args...))
end
Sym{TAG}(::Val{:call}, x::Expr) where {TAG} = convert(Sym{TAG}, Sym(Val(:call), x))

function Sym(::Val{:macrocall}, x::Expr)
    @assert x.head === :macrocall
    return Sym{DEFAULT_TAG}(Val(:noargs), :macrocall, x.args...)
end
Sym{TAG}(::Val{:macrocall}, x::Expr) where {TAG} = convert(Sym{TAG}, Sym(Val(:macrocall), x))

symignore(x) = false
symignore(::LineNumberNode) = true

# Make Syms callable

(f::Sym)(args...) = apply(f, map(arg -> convert(Sym, arg), args)...)
(f::Sym)(args::Sym...) = apply(f, args...)

##################################################
# Convenience methods
##################################################

macro S_str(str::String, suffix::String = "")
    ex = esc_sym(Meta.parse(str))
    if suffix == "u"
        return :(unwrap($ex))
    end
    return ex
end

sym(x) = Sym{DEFAULT_TAG}(x)
sym(tag, x) = Sym{tag}(x)

macro sym(xs...)
    TAG = DEFAULT_TAG
    __return__ = Expr(:block)
    __return_tuple__ = Expr(:tuple)
    for x in xs
        if Meta.isexpr(x, :vect)
            if length(x.args) == 0
                TAG = DEFAULT_TAG
            elseif length(x.args) == 1
                TAG = esc(x.args[1])
            elseif length(x.args) > 1
                TAG = Expr(:curly, :Union, map(esc, x.args)...)
            else
                error("malformated tag: $x")
            end
            continue
        end
        @assert x isa Symbol
        x, symbol = esc(x), QuoteNode(x)
        push!(__return__.args, :($x = sym($TAG, $symbol)))
        push!(__return_tuple__.args, x)
    end
    push!(__return__.args, __return_tuple__)
    return __return__
end

##################################################
# Macro tools
##################################################

esc_sym(x; kwargs...) = Expr(:call, Sym, esc(x))
function esc_sym(x::Expr; interpolate=true, kwargs...)
    ex = interpolate ? Expr(:quote, unblock_interpolate(x; interpolate=interpolate, kwargs...)) : QuoteNode(x)
    return Expr(:call, Sym, ex)
end

unblock_interpolate(x; kwargs...) = x
unblock_interpolate(x::Symbol; interpolate=true, kwargs...) = interpolate ? Expr(:$, x) : x
function unblock_interpolate(x::QuoteNode; interpolate=true, kwargs...)
    x = Expr(:call, Sym, x)
    return interpolate ? Expr(:$, x) : x
end

function unblock_interpolate(x::Expr; kwargs...)
    return unblock_interpolate(Val(x.head), x; kwargs...)
end

function unblock_interpolate(::Val{head}, x::Expr; kwargs...) where {head}
    @assert x.head === head
    return unblock(Expr(head, map(arg -> unblock_interpolate(arg; kwargs...), x.args)...))
end

function unblock_interpolate(::Val{:call}, x::Expr; interpolate=true, kwargs...)
    @assert x.head === :call
    x.args = map(arg -> unblock_interpolate(arg; interpolate=false, kwargs...), x.args)
    ex = :(applicable($(x.args...)) ? $x : Core._expr(:call, $(x.args...)))
    return interpolate ? Expr(:$, ex) : ex
end

function unblock_interpolate(::Val{:(.)}, x::Expr; kwargs...)
    @assert x.head === :(.)
    return unblock(Expr(:(.), unblock_interpolate(x.args[1]; kwargs...), x.args[2:end]...))
end

function unblock_interpolate(::Val{:quote}, x::Expr; interpolate=true, kwargs...)
    @assert x.head === :quote
    return interpolate ? Expr(:$, x) : x
end

function unblock_interpolate(::Val{:macrocall}, x::Expr; interpolate=true, kwargs...)
    @assert x.head === :macrocall
    return interpolate ? Expr(:$, x) : x
end

function unblock_interpolate(::Val{:(::)}, x::Expr; interpolate=true, kwargs...)
    @assert x.head === :(::)
    if length(x.args) >= 2
        if Meta.isexpr(x.args[2], :(::)) && length(x.args[2].args) == 1
            x.args[1] = unblock_interpolate(x.args[1]; interpolate=false, kwargs...)
            x.args[2] = unblock_interpolate(x.args[2].args[1]; interpolate=false, kwargs...)
            ex = :(Core._expr(:(::), $(x.args...)))
        else
            ex = Expr(:call, sym, x.args[2], unblock_interpolate(x.args[1]; interpolate=false, kwargs...))
        end
    else
        x.args = map(arg -> unblock_interpolate(arg; interpolate=false, kwargs...), x.args)
        ex = :(Core._expr(:(::), $(x.args...)))
    end
    return interpolate ? Expr(:$, ex) : ex
end

##################################################
# Comparison
##################################################

Base.isequal(x::Sym, y::Sym) = false
function Base.isequal(x::Sym{TAG}, y::Sym{TAG}) where {TAG}
    return gethead(x) === gethead(y) && isequal(getargs(x), getargs(y))
end

##################################################
# Base extensions
##################################################

Base.hash(x::Sym, h::UInt) = hash((gethead(x), getargs(x)), h)
Base.copy(x::Sym{TAG}) where {TAG} = Sym{TAG}(gethead(x), copy(getargs(x))...)
Base.deepcopy(x::Sym{TAG}) where {TAG} = Sym{TAG}(gethead(x), deepcopy(getargs(x))...)
Base.iterate(x::Sym) = (x[1], 2)
Base.iterate(x::Sym, i) = length(x) + 1 === i ? nothing : (x[i], i + 1)

##################################################
# Utilities
##################################################

isatomic(x::Sym) = hashead(x, (:symbol, :fn, :function, :(->), :type, :object))
isfaithful(x::Sym{TAG}) where {TAG} = TAG isa Type && firstarg(x) isa TAG

@inline _typeof(x::Sym)::Sym{Type} = Sym(typeof)(x)

@inline tagof(x::Sym{TAG}) where {TAG} = TAG
@inline tagof(x::T) where {T} = T
@inline tagof(x::QuoteNode) = typeof(x.value)
@inline tagof(x::Symbol) = DEFAULT_TAG
@inline tagof(x::Expr) = tagof(Val(x.head), x)
@inline tagof(::Val, x::Expr) = DEFAULT_TAG
@inline tagof(::Val{:function}, x::Expr) = Function
@inline tagof(::Val{:(->)}, x::Expr) = Function

@inline istypetag(::Sym{TAG}) where {TAG} = TAG isa Type
@inline issubtag(a, b) = false
@inline issubtag(::Any, ::Type{Any}) = true
@inline issubtag(::Type, ::Type{Any}) = true
@inline issubtag(::Type{Any}, ::Type{Any}) = true
@inline issubtag(a::Symbol, b::Symbol) = a === b
@inline issubtag(a::Type, b::Type) = a <: b

@inline gethead(x) = getfield(x, :head)
@inline sethead!(x, val) = (setfield!(x, :head, val); x)
@inline hashead(x, head::Symbol) = gethead(x) === head
@inline hashead(x, heads::Union{Tuple,AbstractVector,AbstractSet}) = any(head -> gethead(x) === head, heads)
@inline hashead(heads) = (x::Sym) -> hashead(x, heads)

@inline getargs(x) = getfield(x, :args)
@inline getargs(x, idx) = getfield(x, :args)[idx]
@inline getargs(x, ::Nothing) = Any[]
@inline setargs!(x, val) = (setfield!(x, :args, val); x)
@inline firstarg(x) = first(getargs(x))
@inline tailargs(x) = view(getargs(x), 2:length(getargs(x)))


text:
using MuZero
using Test

@testset "MuZero.jl" begin
    @test true
end

include("../diagnostics/tac/game.jl")
include("../diagnostics/tac/inference_tests.jl")
dir = "$(@__DIR__)/output"
tests = [
    get_policyvalue_tests(Tac),
    get_reward_tests(Tac),
]
config = getconfig(Tac)
config.save_interval = 1
config.checkpoint_interval = 1
config.batch_size = 10
config.num_filters = 3
config.representation_resnet_layers = 1
config.dynamics_resnet_layers = 1
config.reward_fc_layers = 1
config.policy_fc_layers = 1
config.value_fc_layers = 1

@time @testset "Basic learner save/load" begin
    l1 = Learner(Tac; config=config, dir=dir, tests=tests)
    @test get_steps(l1) == 0
    save!(l1)
    l2 = Learner(Tac; config=config, uuid=l1.uuid, dir=dir, tests=tests)
    @test get_steps(l2) == 0
    @test l1.train.saved_at == l2.train.saved_at
    @test l1.training_started == l2.training_started
    @test l1.stats.losses == l2.stats.losses
    @test l1.stats.accuracies == l2.stats.accuracies
    @test l1.stats.step_count == l2.stats.step_count
end

@time @testset "Basic network training" begin
    l = Learner(Tac; config=config, dir=dir, tests=tests)
    selfplay!(l; limit=10)
    @test length(l.stats.losses) == 0
    @test length(l.stats.accuracies) == 0
    fit!(l, steps=1)
    @test length(l.stats.accuracies) == 2
    @test length(l.stats.losses) == 2
    @test get_steps(l) == 1
    fit!(l, steps=1)
    @test length(l.stats.accuracies) == 3
    @test length(l.stats.losses) == 3
    @test get_steps(l) == 2
end

@time @testset "Training/test networks save/load" begin
    l1 = Learner(Tac; config=config, dir=dir, tests=tests)
    selfplay!(l1; limit=10)
    fit!(l1, steps=1)
    batch_tensor = randtrainingbatch(l1) |> tensor
    train_net_initial_out1 = initial_inference(l1.train, batch_tensor.images) |> NetworkOutput
    test_net_initial_out1 = initial_inference(l1.test, batch_tensor.images) |> NetworkOutput
    loss1, accuracy1 = get_loss(l1), get_accuracy(l1)
    save!(l1)
    l2 = Learner(Tac; config=config, uuid=l1.uuid, dir=dir, tests=tests)
    train_net_initial_out2 = initial_inference(l2.train, batch_tensor.images) |> NetworkOutput
    test_net_initial_out2 = initial_inference(l2.test, batch_tensor.images) |> NetworkOutput
    loss2, accuracy2 = get_loss(l2), get_accuracy(l2)
    @test train_net_initial_out1 == train_net_initial_out2
    @test test_net_initial_out1 == test_net_initial_out2
    @test accuracy1 == accuracy2
    @test loss1 == loss2
end


text:
module SafeBuffers

export withbuffer

export BufferPool

"""
    BufferPool(n, factory)

A type with `n` lockable mutable buffers.

# Examples

```jldoctest
julia> b1 = BufferPool(2, () -> Ref{Int}())
BufferPool{Base.RefValue{Int64}}(🔒: 0 of 2, size = ())

julia> b2 = BufferPool(2, () -> Vector{Int}(undef, 2))
BufferPool{Vector{Int64}}(🔒: 0 of 2, size = (2,))

julia> b3 = BufferPool(Threads.nthreads(), () -> Matrix{Int}(undef, 2, 2))
BufferPool{Vector{Int64}}(🔒: 0 of 8, size = (2,2))
```
"""
struct BufferPool{T}
    N::Int
    condition::Threads.Condition
    buffers::Vector{T}
end

function BufferPool(n, factory)
    n > 0 || throw(ArgumentError("must have at least 1 buffer; got n = $n"))
    BufferPool(n, Threads.Condition(), [ factory() for _ in 1:n ])
end

function Base.show(io::IO, pool::BufferPool)
    print(io, typeof(pool), "(🔒: $(pool.N - length(pool.buffers)) of $(pool.N), size = $(size(first(pool.buffers))))")
end

"""
    withbuffer(f, pool)

Acquire a buffer in the `BufferPool` and call `f` on the acquired buffer.

`withbuffer` will wait as necessary for a buffer to become available. `f` must take one
argument, the buffer to be used/mutated.

# Examples

```jldoctest
julia> withbuffer(b2) do buf
           buf .= 1:2
           sum(abs2, buf)
       end
5
```
"""
function withbuffer(f, pool)
    buf = lock(pool.condition) do
        while isempty(pool.buffers)
            wait(pool.condition)
        end
        pop!(pool.buffers)
    end

    try
        f(buf)
    finally
        lock(pool.condition) do
            push!(pool.buffers, buf)
            notify(pool.condition)
        end
    end
end

end


text:
using Ansillary

struct State
    size
    issues
end

State() = State(Screen.raw(Screen.size), [])

abstract type Stage end

function draw(::Stage, _) end

handle(stage::Stage, _, _) = stage

tick(stage::Stage) = stage

struct Sleep <: Stage
    until
    next
end

tick(stage::Sleep) =
    if time() > stage.until
        stage.next
    else
        stage
    end

struct Introduction <: Stage end

function draw(::Introduction, _)
    Cursor.move!(Cursor.Coordinate(1, 1))

    println("Welcome to Ansillary's interactive test script!")
    println("At any point you can press Ctrl+c to exit the script.")
    println("Otherwise, please follow the instructions on screen.")
    println("Press any key to continue.")
end

function handle(stage::Introduction, _, _)
    Screen.clear!(Screen.All())
    Cursor.move!(Cursor.Coordinate(1, 1))

    SizeTest()
end

struct SizeTest <: Stage end

function draw(::SizeTest, state)
    rows = state.size.rows
    columns = state.size.columns
    println("We have calculated that your terminal is $rows rows by $columns columns.")
    println("Does this look correct? (y/n)")
end

function handle(stage::SizeTest, event, state)
    next = HideTestInitial()

    if event == Inputs.Character('y')
        Screen.clear!(Screen.All())
        Cursor.move!(Cursor.Coordinate(1, 1))

        next
    elseif event == Inputs.Character('n')
        Screen.clear!(Screen.All())
        Cursor.move!(Cursor.Coordinate(1, 1))

        push!(state.issues, "`Screen.size` did not correctly calculate the terminal size.")

        next
    else
        Cursor.move!(Cursor.Coordinate(3, 1))
        Screen.clear!(Screen.CurrentLine())
        println("You must press 'y' or 'n' not: $event")

        stage
    end
end

struct HideTestInitial <: Stage end

function draw(::HideTestInitial, _)
    println("This will test hiding the cursor (it should currently be visible).")
    println("Is the cursor currently visible on the screen? (y/n)")
end

handle(stage::HideTestInitial, event, state) =
    if event == Inputs.Character('y')
        Cursor.move!(Cursor.Coordinate(1, 1))
        Screen.clear!(Screen.CurrentLine())

        HideTestHidden()
    elseif event == Inputs.Character('n')
        Cursor.move!(Cursor.Coordinate(1, 1))
        Screen.clear!(Screen.CurrentLine())

        push!(state.issues, "Cursor was not visible at the beginning of the script.")

        MoveTest()
    else
        Cursor.move!(Cursor.Coordinate(3, 1))
        Screen.clear!(Screen.CurrentLine())
        println("You must press 'y' or 'n' not: $event")

        stage
    end

struct HideTestHidden <: Stage end

function draw(::HideTestHidden, _)
    Cursor.hide!()

    println("The cursor should now have disappeared.")
end

handle(stage::HideTestHidden, event, state) =
    if event == Inputs.Character('y')
        push!(state.issues, "`Cursor.hide!` did not hide the cursor.")

        # In case this was an accident.
        Cursor.show!()

        Cursor.move!(Cursor.Coordinate(1, 1))
        Screen.clear!(Screen.CurrentLine())

        MoveTest()
    elseif event == Inputs.Character('n')
        Cursor.move!(Cursor.Coordinate(1, 1))
        Screen.clear!(Screen.CurrentLine())

        HideTestFinal()
    else
        Cursor.move!(Cursor.Coordinate(3, 1))
        Screen.clear!(Screen.CurrentLine())
        println("You must press 'y' or 'n' not: $event")

        stage
    end

struct HideTestFinal <: Stage end

function draw(::HideTestFinal, _)
    Cursor.show!()

    println("The cursor should now have reappeared.")
end

function handle(stage::HideTestFinal, event, state)
    if event == Inputs.Character('y')
        Screen.clear!(Screen.All())
        Cursor.move!(Cursor.Coordinate(1, 1))

        MoveTest()
    elseif event == Inputs.Character('n')
        Screen.clear!(Screen.All())
        Cursor.move!(Cursor.Coordinate(1, 1))

        push!(state.issues, "`Cursor.show!` did not reveal the cursor.")

        MoveTest()
    else
        Cursor.move!(Cursor.Coordinate(3, 1))
        Screen.clear!(Screen.CurrentLine())
        println("You must press 'y' or 'n' not: $event")

        stage
    end
end

struct MoveTest <: Stage
    step
    tests
end

MoveTest() = MoveTest(1)

MoveTest(step) = MoveTest(
    step,
    [
        (Cursor.Coordinate(8, 4), 5),
        (Cursor.Up(2), 2),
        (Cursor.Left(2), 1),
        (Cursor.Right(4), 3),
        (Cursor.Down(4), 9),
        (Cursor.Column(2), 7),
    ],
)

function draw(stage::MoveTest, _)
    if stage.step == 1
        println("We will now test moving the cursor.")
        println("A grid of numbers will apear below.")
        println("When prompted you should press the number that the cursor is on.")
        println("If the cursor is not on a number then press any other key.")
        println()
        println(" 1 2 3")
        println()
        println(" 4 5 6")
        println()
        println(" 7 8 9")
        println()
        println("What number is the cursor on?")
        println("Step ")
    end

    Cursor.save() do
        Cursor.move!(Cursor.Coordinate(13, 6))
        print(stage.step)
    end

    movement, _ = stage.tests[stage.step]
    Cursor.move!(movement)
end

handle(stage::MoveTest, event, state) =
    if event in Inputs.Character.('0':'9')
        movement, expected = stage.tests[stage.step]

        if event != Inputs.Character('0' + expected)
            push!(
                state.issues,
                "`$movement` did not move the cursor onto the number $expected.",
            )
        end

        if stage.step >= length(stage.tests)
            Screen.clear!(Screen.All())
            Cursor.move!(Cursor.Coordinate(1, 1))

            ClearTestIntroduction()
        else
            MoveTest(stage.step + 1)
        end
    else
        Cursor.save() do
            Cursor.move!(Cursor.Coordinate(14, 1))
            Screen.clear!(Screen.CurrentLine())
            println("You must press a number key not: $event")
        end

        stage
    end

struct ClearTestIntroduction <: Stage end

function draw(::ClearTestIntroduction, _)
    println("Next we will test clearing the screen.")
    println("Your screen will fill with 'x's and 'o's, then there will be a brief pause.")
    println("Occasionally there will be a '|' to make it more obvious where the 'x's end.")
    println("After the pause the 'x's should disappear, leaving all of the 'o's and the '|'.")
    println("There will be another pause and you will be asked whether the 'x's disappeared.")
    println("If only the 'x's disappeared, press 'y'.")
    println("If any of the 'o's or '|' disappeared or there are some 'x's left, press 'n'.")
    println("Press any key to continue.")
end

function handle(::ClearTestIntroduction, _, _)
    Screen.clear!(Screen.All())
    Cursor.move!(Cursor.Coordinate(1, 1))

    ClearTestAll()
end

abstract type ClearTest <: Stage end

struct ClearTestRun <: Stage
    stage::ClearTest
end

struct ClearTestFinish <: Stage
    stage::ClearTest
end

function draw(stage::ClearTest, state)
    setup(stage, state)
end

tick(stage::ClearTest) = Sleep(time() + 2, ClearTestRun(stage))

draw(stage::ClearTestRun, _) = Screen.clear!(area(stage.stage))

tick(stage::ClearTestRun) = Sleep(time() + 2, ClearTestFinish(stage.stage))

function draw(stage::ClearTestFinish, _)
    Screen.clear!(Screen.All())
    Cursor.move!(Cursor.Coordinate(1, 1))

    println("Did the 'x's (and only the 'x's) disappear? (y/n)")
end

handle(stage::ClearTestFinish, event, state) =
    if event == Inputs.Character('y')
        Screen.clear!(Screen.All())
        Cursor.move!(Cursor.Coordinate(1, 1))

        next(stage.stage)
    elseif event == Inputs.Character('n')
        Screen.clear!(Screen.All())
        Cursor.move!(Cursor.Coordinate(1, 1))

        push!(
            state.issues,
            "`$(area(stage.stage))` did not clear the correct area of the screen.",
        )

        next(stage.stage)
    else
        Cursor.move!(Cursor.Coordinate(2, 1))
        Screen.clear!(Screen.CurrentLine())

        println("You must press 'y' or 'n' not: $event")

        stage
    end

row(character, state) = character^state.size.columns
line(character, state) = row(character, state) * '\n'

struct ClearTestAll <: ClearTest end

function setup(::ClearTestAll, state)
    print(line('x', state)^(state.size.rows - 1))
    print(row('x', state))
end

area(::ClearTestAll) = Screen.All()

next(::ClearTestAll) = ClearTestCurrentLine()

struct ClearTestCurrentLine <: ClearTest end

function setup(::ClearTestCurrentLine, state)
    print(line('o', state)^10)
    print(line('x', state))
    print(line('o', state)^(state.size.rows - 12))
    print(row('o', state))

    Cursor.move!(Cursor.Coordinate(11, 1))
end

area(::ClearTestCurrentLine) = Screen.CurrentLine()

next(::ClearTestCurrentLine) = ClearTestFromCursorBack()

struct ClearTestFromCursorBack <: ClearTest end

function setup(::ClearTestFromCursorBack, state)
    print(line('o', state)^15)
    println('x'^10, '|', 'o'^(state.size.columns - 11))
    print(line('o', state)^(state.size.rows - 17))
    print(row('o', state))

    Cursor.move!(Cursor.Coordinate(16, 10))
end

area(::ClearTestFromCursorBack) = Screen.FromCursorBack()

next(::ClearTestFromCursorBack) = ClearTestFromCursorForward()

struct ClearTestFromCursorForward <: ClearTest end

function setup(::ClearTestFromCursorForward, state)
    print(line('o', state)^20)
    println('o'^20, '|', 'x'^(state.size.columns - 21))
    print(line('o', state)^(state.size.rows - 22))
    print(row('o', state))

    Cursor.move!(Cursor.Coordinate(21, 22))
end

area(::ClearTestFromCursorForward) = Screen.FromCursorForward()

next(::ClearTestFromCursorForward) = ClearTestFromCursorUp()

struct ClearTestFromCursorUp <: ClearTest end

function setup(::ClearTestFromCursorUp, state)
    print(line('x', state)^7)
    println('x'^16, '|', 'o'^(state.size.columns - 17))
    print(line('o', state)^(state.size.rows - 9))
    print(row('o', state))

    Cursor.move!(Cursor.Coordinate(8, 16))
end

area(::ClearTestFromCursorUp) = Screen.FromCursorUp()

next(::ClearTestFromCursorUp) = ClearTestFromCursorDown()

struct ClearTestFromCursorDown <: ClearTest end

function setup(::ClearTestFromCursorDown, state)
    print(line('o', state)^17)
    println('o'^8, '|', 'x'^(state.size.columns - 9))
    print(line('x', state)^(state.size.rows - 19))
    print(row('x', state))

    Cursor.move!(Cursor.Coordinate(18, 10))
end

area(::ClearTestFromCursorDown) = Screen.FromCursorDown()

next(::ClearTestFromCursorDown) = PasteTest()

struct PasteTest <: Stage end

function draw(stage::PasteTest, _)
    println("Next we will test whether bracketed paste is working.")
    println("Please paste some (not too much) text into the terminal.")

    Inputs.paste!()
end

handle(::PasteTest, event, state) =
    if event == Inputs.PasteStart()
        PasteTestFindEnd()
    else
        push!(state.issues, "Terminal did not send the `PasteStart` event.")

        Inputs.nopaste!()

        PasteTestEatInput()
    end

struct PasteTestFindEnd <: Stage end

handle(stage::PasteTestFindEnd, event, _) =
    if event == Inputs.PasteEnd()
        Inputs.nopaste!()

        InputTest()
    else
        stage
    end

struct PasteTestEatInput <: Stage
    seen_tick
end

PasteTestEatInput() = PasteTestEatInput(false)

# Since pasted input gets sent all at once if we see two ticks then that's
# probably the end of input.
tick(stage::PasteTestEatInput) =
    if stage.seen_tick
        InputTest()
    else
        PasteTestEatInput(true)
    end

struct InputTest <: Stage
    step
    draw
    tests
end

InputTest(step, draw) = InputTest(
    step,
    draw,
    [
        Inputs.Character('g'),
        Inputs.Character('T'),
        Inputs.Character('5'),
        Inputs.Character('*'),
        Inputs.Insert(),
        Inputs.Delete(),
        Inputs.Home(),
        Inputs.End(),
        Inputs.PageUp(),
        Inputs.PageDown(),
        Inputs.Up(),
        Inputs.Down(),
        Inputs.Left(),
        Inputs.Right(),
        Inputs.Esc(),
        Inputs.Modified(Inputs.Character('d'), [Inputs.Ctrl()]),
        Inputs.Modified(Inputs.Character('a'), [Inputs.Ctrl(), Inputs.Alt()]),
        Inputs.Modified(Inputs.Insert(), [Inputs.Ctrl()]),
        Inputs.Modified(Inputs.PageUp(), [Inputs.Ctrl(), Inputs.Alt()]),
        Inputs.Modified(Inputs.Right(), [Inputs.Ctrl()]),
        Inputs.Modified(Inputs.Left(), [Inputs.Ctrl(), Inputs.Alt()]),
        Inputs.F(1),
        Inputs.Modified(Inputs.F(2), [Inputs.Ctrl()]),
    ],
)

InputTest() = InputTest(1, true)

function draw(stage::InputTest, _)
    if stage.draw
        Cursor.hide!()

        Cursor.move!(Cursor.Coordinate(1, 1))
        Screen.clear!(Screen.All())

        println("The remainder of the test will check input events.")
        println("This will not be an exhaustive test, for that use `examples/keylogger.jl`.")
        println("Please follow the instructions on screen.")

        print("Please press: ")
    end

    Cursor.move!(Cursor.Coordinate(4, 15))
    Screen.clear!(Screen.FromCursorForward())

    print(stage.tests[stage.step])
end

function handle(stage::InputTest, event, state)
    expected = stage.tests[stage.step]
    if event != expected
        push!(state.issues, "$expected received as $event.")
    end

    if stage.step >= length(stage.tests)
        Cursor.show!()
        nothing
    else
        InputTest(stage.step + 1, false)
    end
end

state = State()
if state.size.rows < 24 || state.size.columns < 80
    error("this script assumes a terminal size of at least 24x80")
end

Screen.alternative() do
    previous = nothing
    stage = Introduction()

    for event in Inputs.EventLoop(Inputs.Millisecond(100))
        if event == Inputs.CTRL_C
            break
        end

        if previous != stage
            draw(stage, state)
            previous = stage
        end

        stage = if event == Inputs.Tick()
            tick(stage)
        else
            handle(stage, event, state)
        end

        if stage === nothing
            break
        end
    end
end

if !isempty(state.issues)
    println("Test script found the following issues:")
    for issue in state.issues
        println('\t', issue)
    end
else
    println("Ansillary seems to handle your terminal well!")
end


text:
## Linear Algebra ##

function turing_chol(A::AbstractMatrix, check)
    chol = cholesky(A, check=check)
    (chol.factors, chol.info)
end
function ChainRules.rrule(::typeof(turing_chol), A::AbstractMatrix, check)
    factors, info = turing_chol(A, check)
    function turing_chol_pullback(Ȳ)
        f̄ = Ȳ[1]
        ∂A = ChainRules.chol_blocked_rev(f̄, factors, 25, true)
        return (ChainRules.NO_FIELDS, ∂A, ChainRules.DoesNotExist())
    end
    (factors,info), turing_chol_pullback
end
function turing_chol_back(A::AbstractMatrix, check)
    C, dC_pullback = rrule(turing_chol, A, check)
    function back(Δ)
        _, dC = dC_pullback(Δ)
        (dC, nothing)
    end
    C, back
end

function symm_turing_chol(A::AbstractMatrix, check, uplo)
    chol = cholesky(Symmetric(A, uplo), check=check)
    (chol.factors, chol.info)
end
function ChainRules.rrule(::typeof(symm_turing_chol), A::AbstractMatrix, check, uplo)
    factors, info = symm_turing_chol(A, check, uplo)
    function symm_turing_chol_pullback(Ȳ)
        f̄ = Ȳ[1]
        ∂A = ChainRules.chol_blocked_rev(f̄, factors, 25, true)
        return (ChainRules.NO_FIELDS, ∂A, ChainRules.DoesNotExist(), ChainRules.DoesNotExist())
    end
    return (factors,info), symm_turing_chol_pullback
end
function symm_turing_chol_back(A::AbstractMatrix, check, uplo)
    C, dC_pullback = rrule(symm_turing_chol, A, check, uplo)
    function back(Δ)
        _, dC = dC_pullback(Δ)
        (dC, nothing, nothing)
    end
    C, back
end

# Tracker's implementation of ldiv isn't good. We'll use Zygote's instead.
zygote_ldiv(A::AbstractMatrix, B::AbstractVecOrMat) = A \ B


text:
struct AirType <: AbstractGlass end
Base.show(io::IO, ::AirType) = print(io, "Air")
glassid(::AirType) = GlassID(AIR, 0)
glassname(::AirType) = "GlassCat.Air"

isair(::AirType) = true
isair(::AbstractGlass) = false
isair(a::GlassID) = a.type === AIR

function info(io::IO, ::AirType)
    println(io, "GlassCat.Air")
    println(io, "Material representing air, RI is always 1.0 at system temperature and pressure, absorption is always 0.0.")
end

"""
Special glass to represent air. Refractive index is defined to always be 1.0 for any temperature and pressure (other indices are relative to this).
"""
const Air = AirType()


text:
function convertToGensysOut(bb, phi, theF, cc, g0, g1, psi, ncpi)

    (nr, nc) = size(g1)
    (nrpsi,ncpsi) = size(psi)
    stateDim = size(bb, 2) - ncpi
    G1 = bb[1:nr, 1:nc]

    ststate = (g0 - g1) \ cc
    CC = (Compat.Matrix(I, nr, nr) - G1) * ststate

    thePsi = vcat(psi, zeros(ncpi, ncpsi))
    aa = phi * thePsi
    impact = aa[1:nr, :]

    # no unique way to represent these components
    (ywt, fmat, fwt) = smallF(theF, aa, stateDim)

    return CC, G1, impact, ywt, fmat, fwt

end # function


text:
using Images, ImageView, ImageFiltering

misty = load("misty.jpg")
#= You can find misty.jpg in the resources folder,
make sure to put it in the same directory you run the code =#

# let's blur misty
# calling reflect will convert correlation to a convolution
kernel = reflect(Kernel.gaussian(6))
# a kernel is a 2d matrix where the sum of the values are (usually) 1
blurred_misty = imfilter(misty, kernel)
imview(blurred_misty)
# imfilter function conlolves the kernel over the image

# we can see that misty is now blurred
# what this does: it takes the values from the kernel and multiplies them by corresponding RGB values
# to modify the pixel. it takes the sum of the multiplied pixel values to replace the center pixel value with it

# what happens is that for every pixel p, it multiplies the RGB values of the center pixel by the center value from the kernel,
# multiplies the pixel one to the left by the value one left to the center in the kernel and so on
# essentially think of it as a grid marching on the image.
# different kernels give out different results

# check out this awesome video if you want to learn more about how this works
# https://www.youtube.com/watch?v=8rrHTtUzyZA&t=960s

# that's actually it for blurring. very simple to implement, hard to understand, haha

max_stars_repo_name: gunalpms/imageprojl

max_issues_repo_name: gunalpms/imageprojl

max_forks_repo_name: gunalpms/imageprojl


text:
let
    gw = GW._test_gw()

    ptag = "HBKJSDBHK"
    pid = getpid()
    rfile = GW._reg_proc!(gw, pid, ptag)
    @test isfile(rfile)
    sdat1 = GW.get_sdat(gw, rfile)
    sdat2 = GW._find_proc_reg(gw, pid)
    sdat3 = GW._find_proc_reg(gw, ptag)
    empty!(gw.dat)
    GW._up_procs_from_disk(gw)
    sdat4 = GW.get_sdat(gw, rfile)

    for dat in [sdat1, sdat2, sdat3, sdat4]
        @test haskey(dat, "lstart")
    end

    @test GW._validate_proc(gw, pid)
    @test GW._validate_proc(gw, ptag)
    @test GW._validate_proc(gw, rfile)

    GW.worker_dir(gw) |> GW._rm
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers