text:
export ParamBox, inValOf, outValOf, inSymOf, outSymOf, dataOf, mapOf, outValCopy,
       inVarCopy, enableDiff!, disableDiff!, isDiffParam, toggleDiff!, gradHFenegy

using LinearAlgebra: eigen
using Symbolics: Num

"""

    ParamBox{T, V, F} <: DifferentiableParameter{ParamBox, T}

Parameter container that can enable parameter differentiations.

≡≡≡ Field(s) ≡≡≡

`data::Array{T<:Number, 0}`: The data (parameter) stored in a 0-D `Array` that can be
accessed by syntax `[]`.

`dataName::Symbol`: The name assigned to the data.

`map::Function`: The mathematical mapping of the data. The mapped result can be accessed by
syntax `()`.

`canDiff::Array{Bool, 0}`: Indicator that whether this container (mapped variable) is
marked as "differentiable".

`index::Array{<:Union{Int, Nothing}, 0}`: Additional index assigned to the parameter.

≡≡≡ Initialization Method(s) ≡≡≡

    ParamBox(data::Array{T, 0}, mapFunction::F, canDiff::Array{Bool, 0},
             index::Array{<:Union{Int, Nothing}, 0}, name::Symbol=:undef,
             dataName::Symbol=:undef) where {T<:Number, F<:Function} ->
    ParamBox{T, name, F}

    ParamBox(data::Array{T, 0},
             name::Symbol=:undef, mapFunction::F=itself, dataName::Symbol=:undef;
             canDiff::Bool=true,
             index::Union{Int, Nothing}=nothing) where {T<:Number, F<:Function} ->
    ParamBox{T, name, F}

    ParamBox(data::Number,
             name::Symbol=:undef, mapFunction::F=itself, dataName::Symbol=:undef;
             canDiff::Bool=true,
             index::Union{Int, Nothing}=nothing,
             paramType::Type{<:Number}=Float64) where {F<:Function} ->
    ParamBox{paramType, name, F}

`name` specifies the name of the (mapped) variable the `ParamBox` represents, which helps
with symbolic representation and automatic differentiation.

`mapFunction`: The (mathematical) mapping of the data, which will be stored in the field
`map`. It is for the case where the variable represented by the `ParamBox` is dependent on
another independent variable of which the value is the stored data in the container. After
initializing a `ParamBox`, e.g `pb1 = ParamBox(x, mapFunction=f)`, `pb[]` returns `x`, and
`pb()` returns `f(x)`. `mapFunction` is set to `$(itself)` in default, which is a dummy
function that maps the data to itself.

`canDiff` determines whether the mapped (math) variable is "marked" as differentiable (if
the mapping is a differentiable function), i.e, whether the mapped variable `ParamBox`
generates during the automatic differentiation procedure is treated as a dependent variable
or an independent variable regardless of the mapping relation.

`paramType` specifies the type of the stored parameter to avoid data type mutation.

≡≡≡ Example(s) ≡≡≡

```jldoctest; setup = :(push!(LOAD_PATH, "../../src/"); using Quiqbox)
julia> ParamBox(1.0)
ParamBox{Float64, :undef, :itself}(1.0)[∂][undef]

julia> ParamBox(1.0, :a)
ParamBox{Float64, :a, :itself}(1.0)[∂][a]

julia> ParamBox(1.0, :a, abs)
ParamBox{Float64, :a, :abs}(1.0)[∂][x_a]
```

**NOTE 1:** The rightmost "`[∂][IV]`" in the printing info indicates the differentiability
and the name of the represented independent variable `:IV`. When the `ParamBox` is marked
as "differentiable", "`[∂]`" is in color green (otherwise it's in grey).

**NOTE 2:** It's always the (mapped) variable `V` generated by a `ParamBox{<:Any, V}` that
is used to construct a basis, whereas the underlying independent variable is used to
differentiate the basis (in other words, only when `mapFunction = $(itself)` or
`canDiff = false` is the independent variable same as the mapped variable/represented
parameter).
"""
struct ParamBox{T, V, F} <: DifferentiableParameter{ParamBox, T}
    data::Array{T, 0}
    dataName::Symbol
    map::Function
    canDiff::Array{Bool, 0}
    index::Array{<:Union{Int, Nothing}, 0}

    function ParamBox(data::Array{T, 0}, mapFunction::F, canDiff, index,
                      name::Symbol=:undef, dataName::Symbol=:undef) where
                     {T<:Number, F<:Function}
        flag = (dataName == :undef)
        if typeof(mapFunction) === typeof(itself)
            dName = flag ? name : dataName
        else
            mapFuncStr = mapFunction |> nameOf |> string
            if startswith(mapFuncStr, '#')
                idx = parse(Int, mapFuncStr[2:end])
                fStr = "f_" * string(name) * numToSubs(idx)
                mapFunction = renameFunc(fStr, mapFunction)
            end
            dName = flag  ?  "x_" * string(name) |> Symbol  :  dataName
        end
        f = mapFunction
        new{T, name, nameOf(f)}(data, dName, f, canDiff, index)
    end
end

(pb::ParamBox)() = Base.invokelatest(pb.map, pb.data[])

function ParamBox(data::Array{T, 0}, name::Symbol=:undef, mapFunction::F=itself,
                  dataName::Symbol=:undef; canDiff::Bool=true,
                  index::Union{Int, Nothing}=nothing) where {T<:Number, F<:Function}
    idx = reshape(Union{Int, Nothing}[0], ()) |> collect
    idx[] = index
    ParamBox(data, mapFunction, fill(canDiff), idx, name, dataName)
end

ParamBox(data::Number, name::Symbol=:undef, mapFunction::F=itself, dataName::Symbol=:undef;
         canDiff::Bool=true, index::Union{Int, Nothing}=nothing,
         paramType::Type{<:Number}=Float64) where {F<:Function} =
ParamBox(fill(data |> paramType), name, mapFunction, dataName; canDiff, index)

"""

    inValOf(pb::ParamBox) -> Number

Return the value of stored data (independent variable) of the input `ParamBox`. Equivalent
to `pb[]`.
"""
inValOf(pb::ParamBox) = pb.data[]

"""

    outValOf(pb::ParamBox) -> Number

Return the value of mapped data (dependent variable) of the input `ParamBox`. Equivalent to
`pb()`.
"""
outValOf(pb::ParamBox) = pb()

"""

    inSymOf(pb::ParamBox) -> Symbol

Return the `Symbol` of the stored data (independent variable) of the input `ParamBox`.
"""
inSymOf(pb::ParamBox) = pb.dataName

"""

    outSymOf(pb::ParamBox) -> Symbol

Return the `Symbol` of the mapped data (dependent variable) of the input `ParamBox`.
"""
outSymOf(::ParamBox{<:Any, V}) where {V} = V

"""

    dataOf(pb::ParamBox{T}) where {T<:Number} -> Array{T, 0}

Return the 0-D `Array` of the data stored in the input `ParamBox`.
"""
dataOf(pb::ParamBox) = pb.data

"""

    mapOf(pb::ParamBox{<:Number, F}) where {F<:Function} -> F

Return the mapping function of the input `ParamBox`.
"""
mapOf(pb::ParamBox) = pb.map

"""

    outValCopy(pb::ParamBox{<:Number, V}) -> ParamBox{<:Number, V, :itself}

Return a new `ParamBox` of which the stored data is a **deep copy** of the mapped data from
the input `ParamBox`.
"""
outValCopy(pb::ParamBox) = ParamBox(pb(), outSymOf(pb), canDiff=pb.canDiff[])

"""

    inVarCopy(pb::ParamBox) -> ParamBox{<:Number, <:Any, :itself}

Return a new `ParamBox` of which the stored data is a **shallow copy** of the stored data
from the input `ParamBox`.

≡≡≡ Example(s) ≡≡≡

```jldoctest; setup = :(push!(LOAD_PATH, "../../src/"); using Quiqbox)
julia> pb1 = ParamBox(-2.0, :a, abs)
ParamBox{Float64, :a, :abs}(-2.0)[∂][x_a]

julia> pb2 = inVarCopy(pb1)
ParamBox{Float64, :x_a, :itself}(-2.0)[∂][x_a]

julia> pb1[] == pb2[] == -2.0
true

julia> pb1[] = 1.1
1.1

julia> pb2[]
1.1
```
"""
inVarCopy(pb::ParamBox) = ParamBox(pb.data, inSymOf(pb), canDiff=pb.canDiff[])

const NoDiffMark = superscriptSym['!']

"""

    enableDiff!(pb::ParamBox) -> ParamBox

Mark the input `ParamBox` as "differentiable" and return the marked `ParamBox`.
"""
function enableDiff!(pb::ParamBox)
    pb.canDiff[] = true
    pb
end

"""

    disableDiff!(pb::ParamBox) -> ParamBox

Mark the input `ParamBox` as "non-differentiable" and return the marked `ParamBox`.
"""
function disableDiff!(pb::ParamBox)
    pb.canDiff[] = false
    pb
end

"""

    isDiffParam(pb::ParamBox) -> Bool

Return the Boolean value of if the input `ParamBox` is "differentiable".
"""
isDiffParam(pb::ParamBox) = pb.canDiff[]

"""

    toggleDiff!(pb::ParamBox) -> Bool

Toggle the differentiability (`pb.canDiff[]`) of the input `ParamBox` and return the
altered result.
"""
toggleDiff!(pb::ParamBox) = begin pb.canDiff[] = !pb.canDiff[] end

function deriveBasisFunc(bf::CompositeGTBasisFuncs, par::ParamBox) where {N}
    varDict = getVarDictCore(bf)
    vr = getVar(par)
    info = diffInfo(bf, vr, varDict)
    diffInfoToBasisFunc(bf, info)
end

function oneBodyDerivativeCore(∂bfs::Vector{<:CompositeGTBasisFuncs},
                               bfs::Vector{<:CompositeGTBasisFuncs},
                               X::Matrix{Float64}, ∂X::Matrix{Float64},
                               ʃ::F, isGradient::Bool = false) where {F<:Function}
    dimOfʃ = 1+isGradient*2
    bsSize = ∂bfs |> length
    ∂ʃ = ones(bsSize, bsSize, dimOfʃ)
    ʃab = ones(bsSize, bsSize, dimOfʃ)
    ∂ʃab = ones(bsSize, bsSize, dimOfʃ)
    for i = 1:bsSize, j = 1:i
       ʃab[i,j,:] = ʃab[j,i,:] = ʃ(bfs[i], bfs[j])
    end
    for i = 1:bsSize, j = 1:i
        ∂ʃab[i,j,:] = ∂ʃab[j,i,:] = ʃ(∂bfs[i], bfs[j]) + ʃ(bfs[i], ∂bfs[j])
    end
    @views begin
        for e = 1:dimOfʃ
            for i=1:bsSize, j=1:i
                # X[i,j] == X[j,i]
                ∂ʃ[i,j,e] = ∂ʃ[j,i,e] =
                transpose( X[:,i]) * ∂ʃab[:,:,e] *  X[:,j] +
                transpose(∂X[:,i]) *  ʃab[:,:,e] *  X[:,j] +
                transpose( X[:,i]) *  ʃab[:,:,e] * ∂X[:,j]
            end
        end
    end
    ∂ʃ
end

function twoBodyDerivativeCore(∂bfs::Vector{<:CompositeGTBasisFuncs},
                               bfs::Vector{<:CompositeGTBasisFuncs},
                               X::Matrix{Float64}, ∂X::Matrix{Float64},
                               ʃ::F, isGradient::Bool = false) where {F<:Function}
    dimOfʃ = 1+isGradient*2
    bsSize = ∂bfs |> length
    ∂ʃ = ones(bsSize, bsSize, bsSize, bsSize, dimOfʃ)
    ʃabcd = ones(bsSize, bsSize, bsSize, bsSize, dimOfʃ)
    ʃ∂abcd = ones(bsSize, bsSize, bsSize, bsSize, dimOfʃ)
    for i = 1:bsSize, j = 1:i, k = 1:i, l = 1:(k==i ? j : k)
        ʃabcd[i,j,k,l,:] = ʃabcd[j,i,k,l,:] = ʃabcd[j,i,l,k,:] = ʃabcd[i,j,l,k,:] =
        ʃabcd[l,k,i,j,:] = ʃabcd[k,l,i,j,:] = ʃabcd[k,l,j,i,:] = ʃabcd[l,k,j,i,:] =
        ʃ(bfs[i],  bfs[j],  bfs[k],  bfs[l])
    end
    for i = 1:bsSize, j=1:bsSize, k=1:bsSize, l=1:k
        ʃ∂abcd[i,j,k,l,:] = ʃ∂abcd[i,j,l,k,:] = ʃ(∂bfs[i], bfs[j],  bfs[k],  bfs[l])
    end
    for e=1:dimOfʃ
        # [∂ʃ4[i,j,k,l] == ∂ʃ4[j,i,l,k] == ∂ʃ4[j,i,k,l] != ∂ʃ4[l,j,k,i]
        for i = 1:bsSize, j = 1:i, k = 1:i, l = 1:(k==i ? j : k)
            val = 0
            # ʃ∂abcd[i,j,k,l,:] == ʃ∂abcd[i,j,l,k,:] ==
            # ʃab∂cd[l,k,i,j,:] == ʃab∂cd[k,l,i,j,:]
            for a = 1:bsSize, b = 1:bsSize, c = 1:bsSize, d = 1:bsSize
                val += (  X[a,i]*X[b,j]* X[c,k]*X[d,l] +  X[a,j]*X[b,i]* X[c,k]*X[d,l] +
                          X[c,i]*X[d,j]* X[a,k]*X[b,l] +  X[c,i]*X[d,j]* X[a,l]*X[b,k]  ) *
                       ʃ∂abcd[a,b,c,d,e] +
                       ( ∂X[a,i]*X[b,j]* X[c,k]*X[d,l] + ∂X[a,j]*X[b,i]* X[c,k]*X[d,l] +
                          X[a,i]*X[b,j]*∂X[c,k]*X[d,l] +  X[a,i]*X[b,j]*∂X[c,l]*X[d,k]  ) *
                       ʃabcd[a,b,c,d,e]
            end
            ∂ʃ[i,j,k,l,e] = ∂ʃ[j,i,k,l,e] = ∂ʃ[j,i,l,k,e] = ∂ʃ[i,j,l,k,e] =
            ∂ʃ[l,k,i,j,e] = ∂ʃ[k,l,i,j,e] = ∂ʃ[k,l,j,i,e] = ∂ʃ[l,k,j,i,e] = val
        end
    end
    ∂ʃ
end

function derivativeCore(bs::Vector{<:CompositeGTBasisFuncs}, par::ParamBox,
                        S::Matrix{Float64}; oneBodyFunc::F1, twoBodyFunc::F2,
                        oneBodyGrad::Bool=false,
                        twoBodyGrad::Bool=false) where {F1<:Function, F2<:Function}
    # ijkl in chemists' notation of spatial bases (ij|kl).
    ∂bfs = deriveBasisFunc.(bs, Ref(par)) |> flatten
    bfs = decompose.(bs) |> flatten
    bsSize = basisSize(bs) |> sum
    ∂S = ones(bsSize, bsSize)
    ∂X = ones(bsSize, bsSize) # ∂X corresponds to the derivative of X = S^(-0.5)
    ∂X₀ = ones(bsSize, bsSize) # ∂X in its eigen basis
    for i=1:bsSize, j=1:i
        S∂ij = overlap(∂bfs[i], bfs[j])
        Si∂j = overlap(bfs[i], ∂bfs[j])
        ∂S[i,j] = ∂S[j,i] = S∂ij[] + Si∂j[]
    end
    X = S^(-0.5) |> Array
    λ, 𝑣 = eigen(S)
    ∂S2 = transpose(𝑣)*∂S*𝑣
    for i=1:bsSize, j=1:i
        ∂X₀[i,j] = ∂X₀[j,i] = (- ∂S2[i,j] * λ[i]^(-0.5) * λ[j]^(-0.5) *
                               (λ[i]^0.5 + λ[j]^0.5)^(-1))
    end
    for i=1:bsSize, j=1:bsSize
        ∂X[j,i] = [𝑣[j,k]*∂X₀[k,l]*𝑣[i,l] for k=1:bsSize, l=1:bsSize] |> sum
    end
    ∂ʃ2 = oneBodyDerivativeCore(∂bfs, bfs, X, ∂X, oneBodyFunc, oneBodyGrad)
    ∂ʃ4 = twoBodyDerivativeCore(∂bfs, bfs, X, ∂X, twoBodyFunc, twoBodyGrad)
    ∂ʃ2, ∂ʃ4
end

function ∂HFenergy(bs::Vector{<:CompositeGTBasisFuncs}, par::ParamBox,
                   C::Union{Matrix{Float64}, NTuple{2, Matrix{Float64}}},
                   S::Matrix{Float64}, mol::Vector{String},
                   nucCoords::Vector{<:AbstractArray},
                   nElectron::Union{Int, NTuple{2, Int}})
    Xinv = S^(0.5)
    Cₓ = (C isa Tuple) ? (Ref(Xinv) .* C) : (Xinv * C)
    ∂hij, ∂hijkl = derivativeCore(bs, par, S,
                                  oneBodyFunc=(i,j)->coreHijCore(i,j,mol,nucCoords),
                                  twoBodyFunc=eeInteractionCore)
    getEᵀ(dropdims(∂hij, dims=3), dropdims(∂hijkl, dims=5), Cₓ, nElectron)
end

function gradHFenegy(bs::Vector{<:CompositeGTBasisFuncs}, par::Vector{<:ParamBox},
                     C::Union{Matrix{Float64}, NTuple{2, Matrix{Float64}}},
                     S::Matrix{Float64}, mol::Vector{String},
                     nucCoords::Vector{<:AbstractArray};
                     nElectron::Union{Int, NTuple{2, Int}}=getCharge(mol))
    if length(C) == 2 && nElectron isa Int
        nElectron = (nElectron÷2, nElectron-nElectron÷2)
    end
    ∂HFenergy.(Ref(bs), par, Ref(C), Ref(S), Ref(mol), Ref(nucCoords), Ref(nElectron))
end

gradHFenegy(bs::Vector{<:CompositeGTBasisFuncs}, par::ParamBox,
            C::Union{Matrix{Float64}, NTuple{2, Matrix{Float64}}}, S::Matrix{Float64},
            mol::Vector{String}, nucCoords::Vector{<:AbstractArray};
            nElectron::Union{Int, NTuple{2, Int}}=getCharge(mol)) =
gradHFenegy(bs, [par], C, S, mol, nucCoords; nElectron)


text:
# TODO: still need to check the spin(+2) or spin(-2) sta
# TODO: summary methods for BlockDiagEquiRect{B} and Adjoint{T,BlockDiagEquiRect{B}}

# Type defs
# ================================================

struct ProjEquiRect{T} <: CartesianProj

    Ny          :: Int
    Nx          :: Int
    θspan       :: Tuple{Float64,Float64}
    φspan       :: Tuple{Float64,Float64}
    θ           :: Vector{Float64}
    φ           :: Vector{Float64}
    θedges      :: Vector{Float64}
    φedges      :: Vector{Float64}
    Ω           :: Vector{Float64}

    storage

end

struct BlockDiagEquiRect{B<:Basis, T<:Real, P<:ProjEquiRect{T}, A<:AbstractArray}  <: ImplicitOp{T}

    blocks :: A
    blocks_sqrt :: Ref{A}
    blocks_pinv :: Ref{A}
    logabsdet :: Ref{Tuple{T,Complex{T}}}
    proj :: P

end

struct AzFourier <: S0Basis end
const  QUAzFourier = Basis2Prod{    𝐐𝐔, AzFourier }
const IQUAzFourier = Basis3Prod{ 𝐈, 𝐐𝐔, AzFourier }
const AzBasis = Union{AzFourier, QUAzFourier, IQUAzFourier}

# Type Alias
# ================================================

make_field_aliases(
    "EquiRect",  ProjEquiRect,
    extra_aliases=OrderedDict(
        "AzFourier"    => AzFourier,
        "QUAzFourier"  => QUAzFourier,
        "IQUAzFourier" => IQUAzFourier,
    ),
)

typealias_def(::Type{<:ProjEquiRect{T}}) where {T} = "ProjEquiRect{$T}"

typealias_def(::Type{F}) where {B,M<:ProjEquiRect,T,A,F<:EquiRectField{B,M,T,A}} = "EquiRect$(typealias(B)){$(typealias(A))}"

# Proj
# ================================================

@memoize function ProjEquiRect(θ, φ, θedges, φedges, θspan, φspan, ::Type{T}, storage) where {T}
    Ny, Nx = length(θ), length(φ)
    Ω  = rem2pi(φedges[2] .- φedges[1], RoundDown) .* diff(.- cos.(θedges))
    ProjEquiRect{T}(Ny, Nx, θspan, φspan, θ, φ, θedges, φedges, Ω, storage)
end

"""
    ProjEquiRect(; Ny::Int, Nx::Int, θspan::Tuple, φspan::Tuple,         T=Float32, storage=Array)
    ProjEquiRect(; θ::Vector, φ::Vector, θedges::Vector, φedges::Vector, T=Float32, storage=Array)

Construct an EquiRect projection object. The projection can either be
specified by:

* The number of pixels `Ny` and `Nx` (corresponding to the `θ` and `φ`
  angular directions, respectively) and the span in radians of the
  field in these directions, `θspan` and `φspan`. The order in which
  the span tuples are given is irrelevant, either order will refer to
  the same field. Note, the spans correspond to the field size between
  outer pixel edges, not from pixel centers. If one wishes to call
  [`Cℓ_to_Cov`](@ref) with this projection, `φspan` must be an integer
  multiple of 2π, but other functionality will be available if this is
  not the case.
* A manual list of pixels centers and pixel edges, `θ`, `φ`, `θedges`,
  `φedges`.

"""
function ProjEquiRect(; T=Float32, storage=Array, kwargs...)

    arg_error() = error("Constructor takes either (θ, φ, θedges, φedges) or (Ny, Nx, θspan, φspan) keyword arguments.")

    if all(haskey.(Ref(kwargs), (:θ, :φ, :θedges, :φedges)))
        !any(haskey.(Ref(kwargs), (:Ny, :Nx, :θspan, :φspan))) || arg_error()
        @unpack (θ, φ, θedges, φedges) = kwargs
        θspan = (θedges[1], θedges[end])
        φspan = (φedges[1], φedges[end])
    elseif all(haskey.(Ref(kwargs), (:Ny, :Nx, :θspan, :φspan)))
        !all(haskey.(Ref(kwargs), (:θ, :φ, :θedges, :φedges))) || arg_error()
        @unpack (Nx, Ny, θspan, φspan) = kwargs
        θspan = (sort(collect(θspan))...,)
        φspan = (sort(collect(φspan))...,)
        φedges = rem2pi.(range(φspan..., length=Nx+1),           RoundDown)
        φ      = rem2pi.(range(φspan..., length=2Nx+1)[2:2:end], RoundDown)
        θedges = range(θspan..., length=Ny+1)
        θ      = range(θspan..., length=2Ny+1)[2:2:end]
    else
        arg_error()
    end

    ProjEquiRect(θ, φ, θedges, φedges, θspan, φspan, real_type(T), storage)

end

# Field Basis
# ================================================
"""
Jperm(ℓ::Int, n::Int) return the column number in the J matrix U^2
where U is unitary FFT. The J matrix looks like this:

|1   0|
|  / 1|
| / / |
|0 1  |

"""
function Jperm(ℓ::Int, n::Int)
    @assert 1 <= ℓ <= n
    ℓ==1 ? 1 : n - ℓ + 2
end

# AzFourier <-> Map
function AzFourier(f::EquiRectMap)
    nφ, T = f.Nx, real(f.T)
    EquiRectAzFourier(m_rfft(f.arr, 2) ./ T(√nφ), f.proj)
end

function Map(f::EquiRectAzFourier)
    nφ, T = f.Nx, real(f.T)
    EquiRectMap(m_irfft(f.arr, nφ, 2) .* T(√nφ), f.proj)
end

# QUAzFourier <-> QUMap
function QUAzFourier(f::EquiRectQUMap)
    nθ, nφ, T = f.Ny, f.Nx, real(f.T)
    qiumap = complex.(f.Qx, f.Ux)
    Uf = m_fft(qiumap, 2) ./ T(√nφ)
    arr = similar(Uf, 2nθ, nφ÷2+1)
    for ℓ = 1:nφ÷2+1
        if (ℓ==1) | ((ℓ==nφ÷2+1) & iseven(nφ))
            arr[1:nθ,     ℓ] .= Uf[:,ℓ]
            arr[nθ+1:2nθ, ℓ] .= conj.(Uf[:,ℓ])
        else
            arr[1:nθ,     ℓ] .= Uf[:,ℓ]
            arr[nθ+1:2nθ, ℓ] .= conj.(Uf[:,Jperm(ℓ,nφ)])
        end
    end
    EquiRectQUAzFourier(arr, f.proj)
end

function QUMap(f::EquiRectQUAzFourier)
    nθₓ2, nφ½₊1 = size(f.arr)
    nθ, nφ, T = f.Ny, f.Nx, real(f.T)
    @assert nφ½₊1 == nφ÷2+1
    @assert 2nθ   == nθₓ2

    pθk = similar(f.arr, nθ, nφ)
    for ℓ = 1:nφ½₊1
        if (ℓ==1) | ((ℓ==nφ½₊1) & iseven(nφ))
            pθk[:,ℓ] .= f.arr[1:nθ,ℓ]
        else
            pθk[:,ℓ]  .= f.arr[1:nθ,ℓ]
            pθk[:,Jperm(ℓ,nφ)] .= conj.(f.arr[nθ+1:2nθ,ℓ])
        end
    end
    qiumap = m_ifft(pθk, 2) .* T(√nφ)
    EquiRectQUMap(cat(real(qiumap), imag(qiumap), dims=3), f.proj)
end

function Base.getindex(f::EquiRectS0, k::Symbol)
    @match k begin
        :Ix => Map(f).arr
        :Il => AzFourier(f).arr
        _ => error("Invalid EquiRectS0 index $k")
    end
end
function Base.getindex(f::EquiRectS2, k::Symbol)
    @match k begin
        :Qx => QUMap(f).Qx
        :Ux => QUMap(f).Ux
        :Px => (qu=QUMap(f); complex.(qu.Qx,qu.Ux))
        :Pl => QUAzFourier(f).arr
        _ => error("Invalid EquiRectS2 index $k")
    end
end

function Base.summary(io::IO, f::EquiRectField)
    @unpack Ny,Nx,Nbatch = f
    print(io, "$(length(f))-element $Ny×$Nx$(Nbatch==1 ? "" : "(×$Nbatch)")-pixel ")
    Base.showarg(io, f, true)
end

# block-diagonal operator
# ================================================

# ## Constructors

function BlockDiagEquiRect{B}(block_matrix::A, proj::P) where {B<:AzBasis, T<:Real, P<:ProjEquiRect{T}, A<:AbstractArray}
    real(eltype(A)) == T || error("Mismatched eltype between $P and $A")
    BlockDiagEquiRect{B,T,P,A}(block_matrix, Ref{A}(), Ref{A}(), Ref{Tuple{T,Complex{T}}}((0,0)), proj)
end

# The following allows construction by a vector of blocks

function BlockDiagEquiRect{B}(vector_of_blocks::Vector{A}, proj::P) where {B<:AzBasis, T<:Real, P<:ProjEquiRect{T}, A<:AbstractMatrix}
    block_matrix = similar(vector_of_blocks[1], size(vector_of_blocks[1])..., length(vector_of_blocks))
    for b in eachindex(vector_of_blocks)
        block_matrix[:,:,b] .= vector_of_blocks[b]
    end
    BlockDiagEquiRect{B}(block_matrix, proj)
end

# ## Linear Algebra: tullio accelerated (operator, field)

# M * f

(*)(M::BlockDiagEquiRect{B}, f::EquiRectField) where {B<:Basis} = M * B(f)

function (*)(M::BlockDiagEquiRect{B}, f::F) where {B<:AzBasis, F<:EquiRectField{B}}
    promote_metadata_strict(M.proj, f.proj) # ensure same projection
    F(@tullio(Bf[p,iₘ] := M.blocks[p,q,iₘ] * f.arr[q,iₘ]), f.proj)
end

(*)(M::Adjoint{T,<:BlockDiagEquiRect{B}}, f::EquiRectField) where {T, B<:Basis} = M * B(f)

function (*)(M::Adjoint{T,<:BlockDiagEquiRect{B}}, f::F) where {T, B<:AzBasis, F<:EquiRectField{B}}
    promote_metadata_strict(M.parent.proj, f.proj) # ensure same projection
    F(@tullio(Bf[p,iₘ] := conj(M.parent.blocks[q,p,iₘ]) * f.arr[q,iₘ]), f.proj)
end

function rrule(::typeof(*), M::BlockDiagEquiRect{B}, f::EquiRectField{B′}) where {B<:Basis, B′<:Basis}
    function times_pullback(Δ)
        BΔ, Bf = B(Δ), B(f)
        Zygote.ChainRules.NoTangent(), @thunk(BlockDiagEquiRect{B}(@tullio(M̄[p,q,iₘ] := Bf.arr[p,iₘ] * conj(BΔ.arr[q,iₘ])), M.proj)'), B′(M' * BΔ)
    end
    M * f, times_pullback
end

# ## Linear Algebra: tullio accelerated (operator, operator)

# M₁ * M₂
function (*)(M₁::BlockDiagEquiRect{B}, M₂::BlockDiagEquiRect{B}) where {B<:AzBasis}
    promote_metadata_strict(M₁.proj, M₂.proj) # ensure same projection
    BlockDiagEquiRect{B}(@tullio(M₃[p,q,iₘ] := M₁.blocks[p,j,iₘ] * M₂.blocks[j,q,iₘ]), M₁.proj)
end

# M₁' * M₂
function (*)(M₁::Adjoint{T,<:BlockDiagEquiRect{B}}, M₂::BlockDiagEquiRect{B}) where {T, B<:AzBasis}
    promote_metadata_strict(M₁.parent.proj, M₂.proj) # ensure same projection
    BlockDiagEquiRect{B}(@tullio(M₃[p,q,iₘ] := conj(M₁.parent.blocks[j,p,iₘ]) * M₂.blocks[j,q,iₘ]), M₁.parent.proj)
end

# M₁ * M₂'
function (*)(M₁::BlockDiagEquiRect{B}, M₂::Adjoint{T,<:BlockDiagEquiRect{B}}) where {T, B<:AzBasis}
    promote_metadata_strict(M₁.proj, M₂.parent.proj) # ensure same projection
    BlockDiagEquiRect{B}(@tullio(M₃[p,q,iₘ] := M₁.blocks[p,j,iₘ] * conj(M₂.parent.blocks[q,j,iₘ])), M₁.proj)
end

# M₁ + M₂, M₁ - M₂, M₁ \ M₂, M₁ / M₂ ... also with mixed adjoints
# QUESTION: some of these may be sped up with @tullio

for op in (:+, :-, :/, :\)

    @eval function Base.$op(M₁::BlockDiagEquiRect{B}, M₂::BlockDiagEquiRect{B}) where {B<:AzBasis}
        promote_metadata_strict(M₁.proj, M₂.proj) # ensure same projection
        BlockDiagEquiRect{B}(
            map( $op, eachslice(M₁.blocks;dims=3), eachslice(M₂.blocks;dims=3) ),
            M₁.proj,
        )
    end

    @eval function Base.$op(M₁::Adjoint{T,<:BlockDiagEquiRect{B}}, M₂::BlockDiagEquiRect{B}) where {T, B<:AzBasis}
        promote_metadata_strict(M₁.parent.proj, M₂.proj) # ensure same projection
        BlockDiagEquiRect{B}(
            map( (m1,m2)->$op(m1',m2), eachslice(M₁.parent.blocks;dims=3), eachslice(M₂.blocks;dims=3) ),
            M₁.proj
        )
    end

    @eval function Base.$op(M₁::BlockDiagEquiRect{B}, M₂::Adjoint{T,<:BlockDiagEquiRect{B}}) where {T, B<:AzBasis}
        promote_metadata_strict(M₁.proj, M₂.parent.proj) # ensure same projection
        BlockDiagEquiRect{B}(
            map( (m1,m2)->$op(m1,m2'), eachslice(M₁.blocks;dims=3), eachslice(M₂.parent.blocks;dims=3) ),
            M₁.proj,
        )
    end

end

for op in (:*, :/)
    @eval Base.$op(a::Scalar, M::BlockDiagEquiRect{B}) where {B} = BlockDiagEquiRect{B}(broadcast($op, a, M.blocks), M.proj)
    @eval Base.$op(M::BlockDiagEquiRect{B}, a::Scalar) where {B} = BlockDiagEquiRect{B}(broadcast($op, a, M.blocks), M.proj)
    @eval Base.$op(a::Scalar, M::Adjoint{T,<:BlockDiagEquiRect{B}}) where {B,T} = (conj(a) * M.parent)'
    @eval Base.$op(M::Adjoint{T,<:BlockDiagEquiRect{B}}, a::Scalar) where {B,T} = (conj(a) * M.parent)'
end

# ## Linear Algebra: with arguments (operator, )

function LinearAlgebra.sqrt(M::BlockDiagEquiRect{B}) where {B<:AzBasis}
    if !isassigned(M.blocks_sqrt)
        M.blocks_sqrt[] = mapslices(sqrt, M.blocks, dims=(1,2))
    end
    BlockDiagEquiRect{B}(M.blocks_sqrt[], M.proj)
end

function LinearAlgebra.pinv(M::BlockDiagEquiRect{B}) where {B<:AzBasis}
    if !isassigned(M.blocks_pinv)
        M.blocks_pinv[] = mapslices(pinv, M.blocks, dims=(1,2))
    end
    BlockDiagEquiRect{B}(M.blocks_pinv[], M.proj)
end

# logdet and logabsdet

function LinearAlgebra.logdet(M::BlockDiagEquiRect{B}) where {B<:AzBasis}
    l, s = logabsdet(M)
    l + log(s)
end

function LinearAlgebra.logabsdet(M::BlockDiagEquiRect{B}) where {B<:AzBasis}
    if M.logabsdet[] == (0,0)
        M.logabsdet[] = mapreduce(logabsdet, ((l1,s1),(l2,s2))->(l1+l2,s1*s2), eachslice(M.blocks, dims=3))
    end
    M.logabsdet[]
end

@adjoint function LinearAlgebra.logabsdet(M::BlockDiagEquiRect)
    logabsdet(M), Δ -> (Δ[1] * pinv(M)',)
end

# dot products

LinearAlgebra.dot(a::EquiRectField, b::EquiRectField) = dot(Ł(a).arr, Ł(b).arr)

# needed by AD
function LinearAlgebra.dot(M₁::Adjoint{T,<:BlockDiagEquiRect{B}}, M₂::BlockDiagEquiRect{B}) where {T, B<:AzBasis}
    (@tullio a[] := conj(M₁.parent.blocks[q,p,iₘ]) * M₂.blocks[p,q,iₘ])[]
end

# mapblocks
# =====================================

function mapblocks(fun::Function, M::BlockDiagEquiRect{B}, f::EquiRectField) where {B<:AzBasis}
    mapblocks(fun, M, B(f))
end

function mapblocks(fun::Function, M::BlockDiagEquiRect{B}, f::F) where {B<:AzBasis, F<:EquiRectField{B}}
    promote_metadata_strict(M.proj, f.proj) # ensure same projection
    Mfarr = similar(f.arr)
    y_    = eachcol(Mfarr)
    x_    = eachcol(f.arr)
    Mb_   = eachslice(M.blocks; dims = 3)
    for (y, x, Mb) in zip(y_, x_, Mb_)
        y .= fun(Mb, x)
    end
    F(Mfarr, f.proj)
end

function mapblocks(fun::Function, Ms::BlockDiagEquiRect{B}...) where {B<:AzBasis}
    map(M->promote_metadata_strict(M.proj, Ms[1].proj), Ms)
    BlockDiagEquiRect{B}(
        map(
            i->fun(getindex.(getproperty.(Ms,:blocks),:,:,i)...), # This looks miserable:(
            axes(Ms[1].blocks,3),
        ),
        Ms[1].proj,
    )
end

# Other methods
# =========================================

# ## simulation

global_rng_for(::Type{<:BlockDiagEquiRect}) = Random.default_rng()

function white_noise(::Type{T}, pj::ProjEquiRect, rng::AbstractRNG=Random.default_rng()) where {T<:Real}
    EquiRectMap(randn(T, pj.Ny, pj.Nx), pj)
end

function white_noise(::Type{T}, pj::ProjEquiRect, rng::AbstractRNG=Random.default_rng()) where {T<:Complex}
    qiu = randn(T, pj.Ny, pj.Nx)
    EquiRectQUMap(real(qiu), imag(qiu), pj)
end

white_noise(ξ::EquiRectField{Map, T},         rng::AbstractRNG) where {T} = white_noise(T, ξ.proj, rng)
white_noise(ξ::EquiRectField{QUMap, T},       rng::AbstractRNG) where {T} = white_noise(Complex{T}, ξ.proj, rng)
white_noise(ξ::EquiRectField{AzFourier, T},   rng::AbstractRNG) where {T} = AzFourier(white_noise(T, ξ.proj, rng))
white_noise(ξ::EquiRectField{QUAzFourier, T}, rng::AbstractRNG) where {T} = QUAzFourier(white_noise(Complex{T}, ξ.proj, rng))

function simulate(rng::AbstractRNG, M::BlockDiagEquiRect{AzFourier,T}) where {T}
    sqrt(M) * white_noise(T, M.proj, rng)
end

function simulate(rng::AbstractRNG, M::BlockDiagEquiRect{QUAzFourier,T}) where {T}
    sqrt(M) * white_noise(Complex{T}, M.proj, rng)
end

# adapt_structure

function adapt_structure(storage, L::BlockDiagEquiRect{B}) where {B}
    BlockDiagEquiRect{B}(adapt(storage, L.blocks), adapt(storage, L.proj))
end

function Base.size(L::BlockDiagEquiRect{<:AzBasis})
    n,m,p = size(L.blocks)
    @assert n==m
    sz = n*p
    return (sz, sz)
end

# covariance and beam operators
# ================================================

function Cℓ_to_Cov(::Val, proj::ProjEquiRect, args...; kwargs...)
    error("Run `using CirculantCov` to use this function.")
end

@init @require CirculantCov="edf8e0bb-e88b-4581-a03e-dda99a63c493" begin

    function Cℓ_to_Cov(::Val{:I}, proj::ProjEquiRect{T}, CI::InterpolatedCℓs; units=1, ℓmax=10_000, progress=true) where {T}

        @unpack θ, φ, Ω = proj
        nθ, nφ  = length(θ), length(φ)
        ℓ       = 0:ℓmax

        CIℓ = nan2zero.(CI(ℓ))

        @assert real(T) == T
        blocks = zeros(T, nθ, nθ, nφ÷2+1)
        # TODO: do we want ngrid as an optional argmuent to Cℓ_to_Cov?
        Γ_I  = CirculantCov.Γθ₁θ₂φ₁φ⃗_Iso(ℓ, CIℓ; ngrid=50_000)
        # using full resolution ComplexF64 for internal construction
        ptmW    = FFTW.plan_fft(Vector{ComplexF64}(undef, nφ))

        pbar = Progress(nθ, progress ? 1 : Inf, "Cℓ_to_Cov: ")
        for k = 1:nθ
            for j = 1:nθ
                Iγⱼₖℓ⃗ = CirculantCov.γθ₁θ₂ℓ⃗(θ[j], θ[k], φ, Γ_I, ptmW)
                for ℓ = 1:nφ÷2+1
                    blocks[j,k,ℓ] = real(Iγⱼₖℓ⃗[ℓ])
                end
            end
            next!(pbar)
        end

        return BlockDiagEquiRect{AzFourier}(blocks, proj)

    end

    function Cℓ_to_Cov(::Val{:P}, proj::ProjEquiRect{T}, CEE::InterpolatedCℓs, CBB::InterpolatedCℓs; units=1, ℓmax=10_000, progress=true) where {T}

        @unpack θ, φ, Ω = proj
        nθ, nφ  = length(θ), length(φ)
        ℓ       = 0:ℓmax

        CBBℓ = nan2zero.(CBB(ℓ))
        CEEℓ = nan2zero.(CEE(ℓ))

        @assert real(T) == T
        blocks = zeros(Complex{T},2nθ,2nθ,nφ÷2+1)
        # TODO: do we want ngrid as an optional argmuent to Cℓ_to_Cov?
        ΓC_EB = CirculantCov.ΓCθ₁θ₂φ₁φ⃗_CMBpol(ℓ, CEEℓ, CBBℓ; ngrid=50_000)
        # using full resolution ComplexF64 for internal construction
        ptmW = FFTW.plan_fft(Vector{ComplexF64}(undef, nφ))

        pbar = Progress(nθ, progress ? 1 : Inf, "Cℓ_to_Cov: ")
        for k = 1:nθ
            for j = 1:nθ
                EBγⱼₖℓ⃗, EBξⱼₖℓ⃗ = CirculantCov.γθ₁θ₂ℓ⃗_ξθ₁θ₂ℓ⃗(θ[j], θ[k], φ, ΓC_EB..., ptmW)
                for ℓ = 1:nφ÷2+1
                    Jℓ = Jperm(ℓ, nφ)
                    blocks[j,    k,    ℓ] = EBγⱼₖℓ⃗[ℓ]
                    blocks[j,    k+nθ, ℓ] = EBξⱼₖℓ⃗[ℓ]
                    blocks[j+nθ, k,    ℓ] = conj(EBξⱼₖℓ⃗[Jℓ])
                    blocks[j+nθ, k+nθ, ℓ] = conj(EBγⱼₖℓ⃗[Jℓ])
                end
            end
            next!(pbar)
        end

        return BlockDiagEquiRect{QUAzFourier}(blocks, proj)

    end

end

function Cℓ_to_Beam(::Val{:I}, proj::ProjEquiRect{T}, CI::InterpolatedCℓs; units=1, ℓmax=10_000, progress=true) where {T}

    @unpack Ω = proj
    Ω′ = T.(Ω)

    Cov = Cℓ_to_Cov(:I, proj, CI; units, ℓmax, progress)
    @tullio Cov.blocks[j,k,iₘ] *= Ω′[k]

    return Cov
end

function Cℓ_to_Beam(::Val{:P}, proj::ProjEquiRect{T}, CI::InterpolatedCℓs; units=1, ℓmax=10_000, progress=true) where {T}

    @unpack θ, Ω = proj
    Ω′ = T.(Ω)

    Cov   = Cℓ_to_Cov(:I, proj, CI; units, ℓmax, progress)
    dcatΩ = Diagonal(vcat(Ω′, Ω′))
    zB    = zeros(T, length(θ), length(θ))

    Beam = BlockDiagEquiRect{QUAzFourier}(
        map(B->[B zB;zB B]*dcatΩ, eachslice(Cov.blocks; dims=3)),
        proj,
    )

    return Beam
end

Cℓ_to_Beam(pol::Symbol, args...; kwargs...) = Cℓ_to_Beam(Val(pol), args...; kwargs...)

# promotion
# ================================================

promote_basis_generic_rule(::Map, ::AzFourier) = Map()

promote_basis_generic_rule(::QUMap, ::QUAzFourier) = QUMap()

# used in broadcasting to decide the resulting metadata when
# broadcasting over two fields
function promote_metadata_strict(metadata₁::ProjEquiRect{T₁}, metadata₂::ProjEquiRect{T₂}) where {T₁,T₂}

    if (
        metadata₁.Ny    === metadata₂.Ny    &&
        metadata₁.Nx    === metadata₂.Nx    &&
        metadata₁.θspan === metadata₂.θspan &&
        metadata₁.φspan === metadata₂.φspan
    )

        # always returning the "wider" metadata even if T₁==T₂ helps
        # inference and is optimized away anyway
        promote_type(T₁,T₂) == T₁ ? metadata₁ : metadata₂

    else
        error("""Can't broadcast two fields with the following differing metadata:
        1: $(select(fields(metadata₁),(:Ny,:Nx,:θspan,:φspan)))
        2: $(select(fields(metadata₂),(:Ny,:Nx,:θspan,:φspan)))
        """)
    end

end

# used in non-broadcasted algebra to decide the resulting metadata
# when performing some operation across two fields. this is free to do
# more generic promotion than promote_metadata_strict (although this
# is currently not used, but in the future could include promoting
# resolution, etc...). the result should be a common metadata which we
# can convert both fields to then do a succesful broadcast
promote_metadata_generic(metadata₁::ProjEquiRect, metadata₂::ProjEquiRect) =
    promote_metadata_strict(metadata₁, metadata₂)

### preprocessing
# defines how ImplicitFields and BatchedReals behave when broadcasted
# with ProjEquiRect fields. these can return arrays, but can also
# return `Broadcasted` objects which are spliced into the final
# broadcast, thus avoiding allocating any temporary arrays.

function preprocess((_,proj)::Tuple{<:Any,<:ProjEquiRect}, r::Real)
    r isa BatchedReal ? adapt(proj.storage, reshape(r.vals, 1, 1, 1, :)) : r
end
# need custom adjoint here bc Δ can come back batched from the
# backward pass even though r was not batched on the forward pass
@adjoint function preprocess(m::Tuple{<:Any,<:ProjEquiRect}, r::Real)
    preprocess(m, r), Δ -> (nothing, Δ isa AbstractArray ? batch(real.(Δ[:])) : Δ)
end

### adapting

# dont adapt the fields in proj, instead re-call into the memoized
# ProjLambert so we always get back the singleton ProjEquiRect object
# for the given set of parameters (helps reduce memory usage and
# speed-up subsequent broadcasts which would otherwise not hit the
# "===" branch of the "promote_*" methods)
function adapt_structure(storage, proj::ProjEquiRect{T}) where {T}
    # TODO: make sure these are consistent with any arguments that
    # were added to the memoized constructor
    @unpack Ny, Nx, θspan, φspan = proj
    T′ = eltype(storage)
    ProjEquiRect(;Ny, Nx, T=(T′==Any ? T : real(T′)), θspan, φspan, storage)
end
adapt_structure(::Nothing, proj::ProjEquiRect{T}) where {T} = proj

# @adjoint function (::Type{F})(arr::A, proj::P) where {B<:SpatialBasis{AzFourier},P<:Proj,T,A<:AbstractArray{T},F<:BaseField{B}}
#     # F(arr, proj), Δ -> (Δ.arr .* adapt(Δ.storage, T.(rfft_degeneracy_fac(proj.Nx))'), nothing)
#     F(arr, proj), Δ -> (Δ.arr, nothing)
# end

ext: jl


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

@doc raw"""DeploymentSpec is the specification of the desired behavior of the Deployment.

    IoK8sApiAppsV1DeploymentSpec(;
        minReadySeconds=nothing,
        paused=nothing,
        progressDeadlineSeconds=nothing,
        replicas=nothing,
        revisionHistoryLimit=nothing,
        selector=nothing,
        strategy=nothing,
        template=nothing,
    )

    - minReadySeconds::Int32 : Minimum number of seconds for which a newly created pod should be ready without any of its container crashing, for it to be considered available. Defaults to 0 (pod will be considered available as soon as it is ready)
    - paused::Bool : Indicates that the deployment is paused.
    - progressDeadlineSeconds::Int32 : The maximum time in seconds for a deployment to make progress before it is considered to be failed. The deployment controller will continue to process failed deployments and a condition with a ProgressDeadlineExceeded reason will be surfaced in the deployment status. Note that progress will not be estimated during the time a deployment is paused. Defaults to 600s.
    - replicas::Int32 : Number of desired pods. This is a pointer to distinguish between explicit zero and not specified. Defaults to 1.
    - revisionHistoryLimit::Int32 : The number of old ReplicaSets to retain to allow rollback. This is a pointer to distinguish between explicit zero and not specified. Defaults to 10.
    - selector::IoK8sApimachineryPkgApisMetaV1LabelSelector : Label selector for pods. Existing ReplicaSets whose pods are selected by this will be the ones affected by this deployment. It must match the pod template's labels.
    - strategy::IoK8sApiAppsV1DeploymentStrategy : The deployment strategy to use to replace existing pods with new ones.
    - template::IoK8sApiCoreV1PodTemplateSpec : Template describes the pods that will be created.
"""
mutable struct IoK8sApiAppsV1DeploymentSpec <: SwaggerModel
    minReadySeconds::Any # spec type: Union{ Nothing, Int32 } # spec name: minReadySeconds
    paused::Any # spec type: Union{ Nothing, Bool } # spec name: paused
    progressDeadlineSeconds::Any # spec type: Union{ Nothing, Int32 } # spec name: progressDeadlineSeconds
    replicas::Any # spec type: Union{ Nothing, Int32 } # spec name: replicas
    revisionHistoryLimit::Any # spec type: Union{ Nothing, Int32 } # spec name: revisionHistoryLimit
    selector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: selector
    strategy::Any # spec type: Union{ Nothing, IoK8sApiAppsV1DeploymentStrategy } # spec name: strategy
    template::Any # spec type: Union{ Nothing, IoK8sApiCoreV1PodTemplateSpec } # spec name: template

    function IoK8sApiAppsV1DeploymentSpec(;minReadySeconds=nothing, paused=nothing, progressDeadlineSeconds=nothing, replicas=nothing, revisionHistoryLimit=nothing, selector=nothing, strategy=nothing, template=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1DeploymentSpec, Symbol("minReadySeconds"), minReadySeconds)
        setfield!(o, Symbol("minReadySeconds"), minReadySeconds)
        validate_property(IoK8sApiAppsV1DeploymentSpec, Symbol("paused"), paused)
        setfield!(o, Symbol("paused"), paused)
        validate_property(IoK8sApiAppsV1DeploymentSpec, Symbol("progressDeadlineSeconds"), progressDeadlineSeconds)
        setfield!(o, Symbol("progressDeadlineSeconds"), progressDeadlineSeconds)
        validate_property(IoK8sApiAppsV1DeploymentSpec, Symbol("replicas"), replicas)
        setfield!(o, Symbol("replicas"), replicas)
        validate_property(IoK8sApiAppsV1DeploymentSpec, Symbol("revisionHistoryLimit"), revisionHistoryLimit)
        setfield!(o, Symbol("revisionHistoryLimit"), revisionHistoryLimit)
        validate_property(IoK8sApiAppsV1DeploymentSpec, Symbol("selector"), selector)
        setfield!(o, Symbol("selector"), selector)
        validate_property(IoK8sApiAppsV1DeploymentSpec, Symbol("strategy"), strategy)
        setfield!(o, Symbol("strategy"), strategy)
        validate_property(IoK8sApiAppsV1DeploymentSpec, Symbol("template"), template)
        setfield!(o, Symbol("template"), template)
        o
    end
end # type IoK8sApiAppsV1DeploymentSpec

const _property_map_IoK8sApiAppsV1DeploymentSpec = Dict{Symbol,Symbol}(Symbol("minReadySeconds")=>Symbol("minReadySeconds"), Symbol("paused")=>Symbol("paused"), Symbol("progressDeadlineSeconds")=>Symbol("progressDeadlineSeconds"), Symbol("replicas")=>Symbol("replicas"), Symbol("revisionHistoryLimit")=>Symbol("revisionHistoryLimit"), Symbol("selector")=>Symbol("selector"), Symbol("strategy")=>Symbol("strategy"), Symbol("template")=>Symbol("template"))
const _property_types_IoK8sApiAppsV1DeploymentSpec = Dict{Symbol,String}(Symbol("minReadySeconds")=>"Int32", Symbol("paused")=>"Bool", Symbol("progressDeadlineSeconds")=>"Int32", Symbol("replicas")=>"Int32", Symbol("revisionHistoryLimit")=>"Int32", Symbol("selector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("strategy")=>"IoK8sApiAppsV1DeploymentStrategy", Symbol("template")=>"IoK8sApiCoreV1PodTemplateSpec")
Base.propertynames(::Type{ IoK8sApiAppsV1DeploymentSpec }) = collect(keys(_property_map_IoK8sApiAppsV1DeploymentSpec))
Swagger.property_type(::Type{ IoK8sApiAppsV1DeploymentSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAppsV1DeploymentSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1DeploymentSpec }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1DeploymentSpec[property_name]

function check_required(o::IoK8sApiAppsV1DeploymentSpec)
    (getproperty(o, Symbol("selector")) === nothing) && (return false)
    (getproperty(o, Symbol("template")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAppsV1DeploymentSpec }, name::Symbol, val)
end


text:
a = 2988348162058574136915891421498819466320163312926952423791023078876139
b = 2351399303373464486466122544523690094744975233415544072992656881240319
m = big(10) ^ 40
@show powermod(a, b, m)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    DigitalOutSet!(dev::DWFDevice, pin::Int, state::Bool)

Set the value of a digital output pin.
"""
function DigitalOutSet!(dev::DWFDevice, pin::Int, state::Bool)
    FDwfDigitalOutEnableSet(dev.hdwf[], pin, Int(state))
    return nothing
end


text:
@testset "CoNLL-U" begin

    trees = collect(test_treebank("english.conllu"))

    @test length(trees) == 2
    @test length.(trees) == [7, 19]

    for C in [ArcStandard(), ArcEager(), ListBasedNonProjective()], tree in trees
        tokens = [t.form for t in tree]
        oracle = Oracle(C, static_oracle)
    end

    for C in [ArcEager(), ArcHybrid()]
        oracle = Oracle(C, static_oracle, untyped)
        for (cfg, t) in Iterators.flatten(oracle.(trees))
            # @test DT.args(t) == ()
        end
    end

    # @test DependencyTrees.noval(CoNLLU).head == -1

    using DependencyTrees: from_conllu
    # make sure the errors get thrown correctly
    @test_throws MultiWordTokenError from_conllu("18-19	cannot	_	_	_	_	_	_	_	SpaceAfter=No")
    @test_throws EmptyTokenError from_conllu("0.1	nothing	_	_	_	_	_	_	_	_")

    c = from_conllu("1	Distribution	distribution	NOUN	S	Number=Sing	7	nsubj	_	_")
    @test c.feats == ["Number=Sing"]
    @test c.lemma == "distribution"

    @test_throws Exception from_conllu("1	2	3")

    sent = """
1	They	they	PRON	PRP	Case=Nom|Number=Plur	2	nsubj	2:nsubj|4:nsubj	_
2	buy	buy	VERB	VBP	Number=Plur|Person=3|Tense=Pres	0	root	0:root	_
3	and	and	CONJ	CC	_	4	cc	4:cc	_
4	sell	sell	VERB	VBP	Number=Plur|Person=3|Tense=Pres	2	conj	0:root|2:conj	_
5	books	book	NOUN	NNS	Number=Plur	2	obj	2:obj|4:obj	_
6	.	.	PUNCT	.	_	2	punct	2:punct	_"""

    graph = deptree(sent, from_conllu)
    for d in graph.tokens
        @test untyped(d) == ()
        @test typed(d) == (d.deprel,)
    end

    c = from_conllu("1	They	they	PRON	PRP	Case=Nom|Number=Plur	2	nsubj	2:nsubj|4:nsubj	_")
    @test length(c.deps) == 2

    sent = """
1	Sue	Sue	_	_	_	2	_	_	_
2	likes	like	_	_	_	0	_	_	_
3	coffee	coffee	_	_	_	2	_	_	_
4	and	and	_	_	_	2	_	_	_
5	Bill	Bill	_	_	_	2	_	_	_
5.1	likes	like	_	_	_	_	_	_	_
6	tea	tea	_	_	_	2	_	_	_
"""

    graph = deptree(sent)
    @test length(graph) == 6

    sent = """
1-2	vámonos	_	_	_	_	_	_	_	_
1	vamos	ir	_	_	_	0	_	_	_
2	nos	nosotros	_	_	_	1	_	_	_
3-4	al	_	_	_	_	_	_	_	_
3	a	a	_	_	_	5	_	_	_
4	el	el	_	_	_	5	_	_	_
5	mar	mar	_	_	_	1	_	_	_
"""
    graph = deptree(sent)
    @test length(graph) == 5
    @test_throws DependencyTrees.MultiWordTokenError from_conllu("1-2	vámonos	_	_	_	_	_	_	_	_")
end


text:
using LyceumDevTools
using LyceumDevTools.TestUtil
using Test

using LyceumCore


text:
@testset "read character matrix" begin
    mktemp() do path, _
        write(path, """
              1####2####
              ####3####4
              #5#6#7#8#9
              """)

        matrix = Advent.read_charmatrix(path)
        @test size(matrix) == (3, 10)
        @test matrix isa Matrix{Char}
        @test matrix[1,1] == '1'
        @test matrix[3,10] == '9'
    end
end


text:
# Choose parameters
c = 1.25
L0 = 27.0
L1 = 27.0

# Choose n points and distributions
m = 251
f0 = pdf(Beta(2.5, 3), linspace(0, 1, m))
f0 = f0 / sum(f0)
f1 = pdf(Beta(3, 2.5), linspace(0, 1, m))
f1 = f1 / sum(f1)  # Make sure sums to 1

# Create an instance of our WaldFriedman class
wf = WaldFriedman(c, L0, L1, f0, f1; m=m);

# Solve and simulate the solution
cdist, tdist = stopping_dist(wf; ndraws=5000)

a = plot([f0 f1],
    xlabel=L"$k$ Values",
    ylabel=L"Probability of $z_k$",
    labels=reshape([L"$f_0$"; L"$f_1$"], 1, 2),
    linewidth=2,
    title="Distributions over Outcomes")

b = plot(wf.pgrid, wf.sol.J,
    xlabel=L"$p_k$",
    ylabel="Value of Bellman",
    linewidth=2,
    title="Bellman Equation")
    plot!(fill(wf.sol.lb, 2), [minimum(wf.sol.J); maximum(wf.sol.J)],
    linewidth=2, color=:black, linestyle=:dash, label="", ann=(wf.sol.lb-0.05, 5., L"\beta"))
    plot!(fill(wf.sol.ub, 2), [minimum(wf.sol.J); maximum(wf.sol.J)],
    linewidth=2, color=:black, linestyle=:dash, label="", ann=(wf.sol.ub+0.02, 5., L"\alpha"),
    legend=:none)

counts = Array(Int64, maximum(tdist))
for i in 1:maximum(tdist)
    counts[i] = sum(tdist .== i)
end
c = bar(counts,
    xticks=0:1:maximum(tdist),
    xlabel="Time",
    ylabel="Frequency",
    title="Stopping Times",
    legend=:none)

counts = Array(Int64, 2)
for i in 1:2
    counts[i] = sum(cdist .== i-1)
end
d = bar([0; 1],
    counts,
    xticks=[0; 1],
    title="Correct Decisions",
    ann=(-.4, 0.6 * sum(cdist), "Percent Correct = $(sum(cdist)/length(cdist))"),
    legend=:none)

plot(a, b, c, d, layout=(2, 2), size=(1200, 800))


text:
module BenchmarkExtTest

for file in sort([file for file in readdir(@__DIR__) if
                  occursin(r"^test.*\.jl$", file)])
    include(file)
end

end # module


text:
@testset "Constrains" begin

   Unordered_con_c_e = UnOrdered_Constrain([Mod("c") , Mod("e")])
   Ordered_con_a_b = Ordered_Constrain([1 3],[Mod("a") Mod("b")])
   Ordered_con_c_b = Ordered_Constrain([2 3],[Mod("c") Mod("b")])

   @testset "single addition" begin
      @test isa(Unordered_con_c_e,UnOrdered_Constrain) == true
      @test isa(Ordered_con_a_b,Ordered_Constrain) == true
      @test isa((Ordered_con_c_b +  Ordered_con_a_b),Compose_Construct_Constrains) == true
      @test (Ordered_con_c_b +  Ordered_con_a_b).construct_con == [Ordered_con_c_b;Ordered_con_a_b]
      @test (Unordered_con_c_e +  Ordered_con_a_b).construct_con == [Unordered_con_c_e;Ordered_con_a_b]
   end

   @testset "mulitple_additions" begin
      @test  Set((Ordered_con_c_b +  Ordered_con_a_b + Unordered_con_c_e).construct_con) == Set([Ordered_con_c_b;Ordered_con_a_b;Unordered_con_c_e])
      @test  isa(sum([Ordered_con_c_b +  Ordered_con_a_b + Unordered_con_c_e]),Compose_Construct_Constrains) == true
   end

end


text:
using DifferentialEquations
using Plots
using Random
using BenchmarkTools

function infection_rate(u,p,t)
    (S,I,R) = u
    (β,c,τ) = p
    N = S+I+R
    β*c*I/N*S
end

function infection!(integrator)
    (β,c,τ) = integrator.p
    integrator.u[1] -= 1
    integrator.u[2] += 1

    # queue recovery callback
    add_tstop!(integrator, integrator.t + τ)
end

infection_jump = ConstantRateJump(infection_rate,infection!);

function recovery_condition(u,t,integrator)
    t == integrator.tstops[1]
end

function recovery!(integrator)
    if integrator.u[2] > 0
        integrator.u[2] -= 1
        integrator.u[3] += 1

        reset_aggregated_jumps!(integrator)
        popfirst!(integrator.tstops)
        integrator.tstops_idx -= 1
    end
end

recovery_callback = DiscreteCallback(recovery_condition, recovery!, save_positions = (false, false))

function affect_initial_recovery!(integrator)
    integrator.u[2] -= u0[2]
    integrator.u[3] += u0[2]

    reset_aggregated_jumps!(integrator)
end

cb_initial_recovery = DiscreteCallback((u,t,integrator) -> t == p[3], affect_initial_recovery!)

tmax = 40.0
tspan = (0.0,tmax);

δt = 0.1
t = 0:δt:tmax;

u0 = [990,10,0]; # S,I,R

p = [0.05,10.0,4.0]; # β,c,τ

Random.seed!(1234);

prob = DiscreteProblem(u0,tspan,p);

prob_jump = JumpProblem(prob, Direct(), infection_jump);

sol_jump = solve(prob_jump, SSAStepper(), callback = CallbackSet(cb_initial_recovery, recovery_callback), tstops = [p[3]]);

out_jump = sol_jump(t);

plot(
    out_jump,
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number"
)

integrator = init(prob_jump,SSAStepper(), callback = recovery_callback);
for i in 1:10
	add_tstop!(integrator, integrator.t + p[3])
end
solve!(integrator)
sol_jump2 = integrator.sol

@benchmark solve(prob_jump, SSAStepper(), callback = CallbackSet(cb_initial_recovery, recovery_callback), tstops = [p[3]]);


text:
## Common utility functions

@static if VERSION < v"1.8-"
    # copy-pasted from the implementation of @lazy_str in Base
    macro lazy_str(text)
        parts = Any[]
        lastidx = idx = 1
        while (idx = findnext('$', text, idx)) !== nothing
            lastidx < idx && push!(parts, text[lastidx:idx-1])
            idx += 1
            expr, idx = Meta.parseatom(text, idx; filename=string(__source__.file))
            push!(parts, esc(expr))
            lastidx = idx
        end
        lastidx <= lastindex(text) && push!(parts, text[lastidx:end])
        :(string($(parts...)))
    end
end

"""
    @toggleassert expression

Internal macro used to assert and expression conditionally on a build-time constant.
To toggle on or off these assertions, the constant has to be modified in the source code
and the module rebuilt afterwards.
"""
macro toggleassert(ex...)
    # Change the following boolean and reload CrystalNets to enable/disable assertions
    ENABLE_ASSERT = true
    if ENABLE_ASSERT
        return esc(quote
            @assert $(ex...)
        end)
    end
    return nothing
end

function no_metadata_metafmt(level::Logging.LogLevel, _module, group, id, file, line)
    @nospecialize
    color = Logging.default_logcolor(level)
    prefix = string(level == Logging.Warn ? "Warning" : string(level), ':')
    return color, prefix, ""
end
const minimal_logger = Logging.ConsoleLogger(; meta_formatter=no_metadata_metafmt)

macro ifwarn(ex)
    return quote
        if (DOWARN[]::Bool)
            @static if VERSION < v"1.7-"
                $(esc(ex))
            else
                with_logger(minimal_logger) do
                    $(esc(ex))
                end
            end
        end
    end
end

function recursive_readdir!(stored, prefix, path)
    for f in readdir(path; sort=false, join=false)
        name = joinpath(prefix, f)
        file = joinpath(path, f)
        if isdir(file)
            recursive_readdir!(stored, name, file)
        elseif !isempty(splitext(f)[2]) # otherwise, not a crystallographic file
            push!(stored, name)
        end
    end
    nothing
end
function recursive_readdir(path)
    ret = String[]
    while isdirpath(path)
        path = path[1:end-1]
    end
    sizehint!(ret, length(readdir(path; sort=false, join=false)))
    recursive_readdir!(ret, "", path)
    return ret
end

function tmpexportname(path, pre, name, ext)
    if name isa Nothing
        return tempname(path; cleanup=false)*ext
    end
    i = 0
    pre = pre*(name::String)*'_'
    x = pre*string(i)*ext
    paths = Set{String}(readdir(path; sort=false, join=false))
    while x ∈ paths
        i += 1
        x = pre*string(i)*ext
    end
    return joinpath(path, x)
end

"""
    export_default(c::Union{PeriodicGraph,CrystalNet,Crystal}, obj=nothing, name=nothing, path=tempdir(); repeats=nothing)

Export a VTF representation of an object at the given `path`.

`obj` is a `String` describing the nature of the object, such as "net", "clusters" or
"subnet" for example. Default is `string(typeof(c))`.

`name` is a `String` inserted in the exported file name. Default is a [`tempname`](https://docs.julialang.org/en/v1/base/file/#Base.Filesystem.tempname).

`repeats` is the maximum distance between a represented atom out of the unit cell and one
inside. Default is between 2 and 6, depending on `obj` and the size of the graph.
"""
export_default(g::PeriodicGraph, args...; kwargs...) = export_default(CrystalNet(g), args...; kwargs...)
function export_default(c, obj=nothing, name=nothing, path=tempdir(); repeats=nothing)
    _repeats = repeats isa Integer ? repeats : begin
        if obj isa AbstractString && (obj == "net" || startswith(obj, "clusters") || startswith(obj, "subnet"))
            2
        else
            nv(c.pge.g) == 0 ? 1 : clamp(fld(600, nv(c.pge.g)), 2, 6)
        end
    end
    if !isempty(path)
        _name = tmpexportname(path, (obj isa Nothing ? string(typeof(c)) : obj)*'_', name, ".vtf")
        truepath = replace(joinpath(path, _name), ('\\' => "/"))
        if obj isa Nothing
            println(lazy"Saving file (representing $(typeof(c))) at $truepath")
        else
            println(lazy"Export of $obj is enabled: saving file at $truepath")
        end
        try
            export_vtf(truepath, c, _repeats)
        catch e
            if e isa SystemError
                @error lazy"Failed to export because of the following error: $e"
            else
                rethrow()
            end
        end
    end
    nothing
end

function db_options(; kwargs...)
    restore_warns = false
    if DOWARN[] && !get(kwargs, :force_warn, false)
        restore_warns = true
        @ifwarn begin
            @warn "Warnings may critically alter performance for this operation and were turned off."
            @info "Use CrystalNets.toggle_warning(false) or the --no-warn option to remove this warning"
            @info "Use the force_warn option to force printing warnings"
        end
        DOWARN[] = false
    end
    # if kwargs explicitly set one of the two, it will take precedence
    return Options(; export_input="", export_attributions="", export_clusters="",
                     export_net="", export_subnets="", kwargs...),
           restore_warns
end

const atomic_numbers = IdDict{Symbol, Int}(
    :H => 1,
    :D => 1,
    :He => 2,
    :Li => 3,
    :Be => 4,
    :B => 5,
    :C => 6,
    :N => 7,
    :O => 8,
    :F => 9,
    :Ne => 10,
    :Na => 11,
    :Mg => 12,
    :Al => 13,
    :Si => 14,
    :P => 15,
    :Pc => 15, # fictitious atom representing an organic P
    :S => 16,
    :Ss => 16, # fictitious atom representing an organic S
    :Cl => 17,
    :Ar => 18,
    :K => 19,
    :Ca => 20,
    :Sc => 21,
    :Ti => 22,
    :V => 23,
    :Cr => 24,
    :Mn => 25,
    :Fe => 26,
    :Co => 27,
    :Ni => 28,
    :Cu => 29,
    :Zn => 30,
    :Ga => 31,
    :Ge => 32,
    :As => 33,
    :Se => 34,
    :Br => 35,
    :Kr => 36,
    :Rb => 37,
    :Sr => 38,
    :Y => 39,
    :Zr => 40,
    :Nb => 41,
    :Mo => 42,
    :Tc => 43,
    :Ru => 44,
    :Rh => 45,
    :Pd => 46,
    :Ag => 47,
    :Cd => 48,
    :In => 49,
    :Sn => 50,
    :Sb => 51,
    :Te => 52,
    :I => 53,
    :Xe => 54,
    :Cs => 55,
    :Ba => 56,
    :La => 57,
    :Ce => 58,
    :Pr => 59,
    :Nd => 60,
    :Pm => 61,
    :Sm => 62,
    :Eu => 63,
    :Gd => 64,
    :Tb => 65,
    :Dy => 66,
    :Ho => 67,
    :Er => 68,
    :Tm => 69,
    :Yb => 70,
    :Lu => 71,
    :Hf => 72,
    :Ta => 73,
    :W => 74,
    :Re => 75,
    :Os => 76,
    :Ir => 77,
    :Pt => 78,
    :Au => 79,
    :Hg => 80,
    :Tl => 81,
    :Pb => 82,
    :Bi => 83,
    :Po => 84,
    :At => 85,
    :Rn => 86,
    :Fr => 87,
    :Ra => 88,
    :Ac => 89,
    :Th => 90,
    :Pa => 91,
    :U => 92,
    :Np => 93,
    :Pu => 94,
    :Am => 95,
    :Cm => 96,
    :Bk => 97,
    :Cf => 98,
    :Es => 99,
    :Fm => 100,
    :Md => 101,
    :No => 102,
    :Lr => 103,
    :Rf => 104,
    :Db => 105,
    :Sg => 106,
    :Bh => 107,
    :Hs => 108,
    :Mt => 109,
    :Ds => 110,
    :Rg => 111,
    :Cn => 112,
    :Nh => 113,
    :Fl => 114,
    :Mc => 115,
    :Lv => 116,
    :Ts => 117,
    :Og => 118,
    :Uue => 119
)

const element_categories = Symbol[ # populated using PeriodicTable.jl
    :nonmetal, :noble, :metal, :metal, :metalloid, :nonmetal, :nonmetal,
    :nonmetal, :halogen, :noble, :metal, :metal, :metal, :metalloid, :nonmetal,
    :nonmetal, :halogen, :noble, :metal, :metal, :metal, :metal, :metal, :metal,
    :metal, :metal, :metal, :metal, :metal, :metal, :metal, :metalloid,
    :metalloid, :nonmetal, :halogen, :noble, :metal, :metal, :metal, :metal,
    :metal, :metal, :metal, :metal, :metal, :metal, :metal, :metal, :metal,
    :metal, :metalloid, :metalloid, :halogen, :noble, :metal, :metal,
    :lanthanide, :lanthanide, :lanthanide, :lanthanide, :lanthanide,
    :lanthanide, :lanthanide, :lanthanide, :lanthanide, :lanthanide,
    :lanthanide, :lanthanide, :lanthanide, :lanthanide, :lanthanide, :metal,
    :metal, :metal, :metal, :metal, :metal, :metal, :metal, :metal, :metal,
    :metal, :metal, :metal, :halogen, :noble, :metal, :metal, :actinide,
    :actinide, :actinide, :actinide, :actinide, :actinide, :actinide, :actinide,
    :actinide, :actinide, :actinide, :actinide, :actinide, :actinide, :actinide,
    :metal, :metal, :metal, :metal, :metal, :metal, :metal, :metal, :metal,
    :metal, :metal, :metal, :metal, :halogen, :noble, :metal]

function string_atomtype(t::Symbol)
    @static if VERSION < v"1.7-"
        replace(replace(string(t), "Pc" => 'P'), "Ss" => 'S')
    else
        replace(string(t), "Pc" => 'P', "Ss" => 'S')
    end
end

function representative_atom(t::Symbol, default::Int=0)
    at = get(atomic_numbers, t, 0)
    at == 0 || return t, at
    styp = string(t)
    @toggleassert length(styp) ≥ 2
    t = islowercase(styp[2]) ? Symbol(styp[1:2]) : Symbol(styp[1])
    return t, get(atomic_numbers, t, default)
end

"""
    issingular(x::SMatrix{N,N,T}) where {N,T<:Rational}

Test whether a NxN matrix is singular.
"""
function issingular(x::SMatrix{N,N,T}) where {N,T<:Rational}
    try
        return iszero(det(x))
    catch e
        e isa OverflowError || rethrow()
    end
    return iszero(det(SMatrix{N,N,widen(T)}(x)))
end

function issingular(x::SMatrix{3,3,T,9})::Bool where T<:Rational
@inbounds begin
    (i, j, k) = iszero(x[1,1]) ? (iszero(x[1,2]) ? (3,1,2)  : (2,1,3)) : (1,2,3)
    iszero(x[1,i]) && return true
    x1i = x[1,i]
    factj = (x[1,j] // x1i)
    factk = (x[1,k] // x1i)
    y11 = x[2,j] - factj * x[2,i]
    y12 = x[3,j] - factj * x[3,i]
    y21 = x[2,k] - factk * x[2,i]
    y22 = x[3,k] - factk * x[3,i]
    return widemul(y11, y22) == widemul(y12, y21)
    # This can overflow so the input matrix should already have a wide enough type
end
end

function issingular(x::SMatrix{2,2,T,4})::Bool where T<:Rational
@inbounds begin
    if iszero(x[1,1])
        (iszero(x[1,2]) || iszero(x[2,1])) && return true
        return false
    end
    x12 = x[1,2] // x[1,1]
    return x[2,2] == widemul(x[2,1], x12)
    # This can overflow so the input matrix should already have a wide enough type
end
end

function issingular(x::SMatrix{1,1,T})::Bool where T<:Rational
    iszero(@inbounds x[1,1])
end

function isrank3(x::Matrix{T}) where T<:Rational
    _n, n = size(x)
    @toggleassert _n == 3
    n < 3 && return false
    cols = collect(eachcol(x))
    u1 = popfirst!(cols)
    u2 = u1
    while iszero(u1) && !isempty(cols)
        u1 = popfirst!(cols)
    end
    n = length(cols)
    k = iszero(@inbounds u1[1]) ? iszero(@inbounds u1[2]) ? 3 : 2 : 1
    i = 1
    @inbounds while i < n
        u2 = cols[i]
        (u2[k]//u1[k])*u1 == u2 || break
        i+=1
    end
    i == n && return false
    @inbounds for j in i+1:n
        a = SMatrix{3,3,T,9}(hcat(u1, u2, cols[j]))
        issingular(a) || return true
    end
    return false
end

function isrank2(x::Matrix{T}) where T<:Rational
    _n, n = size(x)
    @toggleassert _n == 2
    n < 2 && return false
    cols = collect(eachcol(x))
    u1 = popfirst!(cols)
    u2 = u1
    while iszero(u1) && !isempty(cols)
        u1 = popfirst!(cols)
    end
    n = length(cols)
    k = iszero(@inbounds u1[1]) ? 2 : 1
    i = 1
    @inbounds while i < n
        u2 = cols[i]
        (u2[k]//u1[k])*u1 == u2 || break
        i+=1
    end
    i == n && return false
    return true
end

function isrank1(x::AbstractMatrix{T}) where T<:Rational
    _n, n = size(x)
    @toggleassert _n == 1
    @inbounds for y in eachcol(x)
        iszero(y[1]) || return true
    end
    return false
end

"""
    back_to_unit(r::Rational)

Return `x - floor(Int, x)`
"""
function back_to_unit(r::Rational)
    return Base.unsafe_rational(mod(numerator(r), denominator(r)), denominator(r))
end

"""
    nextword(l::String, k::Int)

Return the triplet of indices `(i, j, x)` such that `l[i:j]` is the next word in the
string `l` after position `k`.
Use `k = x` to get the following word, and so forth.

`(0, 0, 0)` is returned if there is no next word.
"""
function nextword(l::String, i::Int)
    n = lastindex(l)
    i == n && return (0, 0, 0)

    i::Int = nextind(l, i)
    start::Int = 0
    while i <= n
        if isspace(l[i])
            i = nextind(l, i)
        else
            if l[i] == '#'
                i = findnext(isequal('\n'), l, i)
                i === nothing && return (0, 0, 0)
            else
                start = i
                break
            end
        end
    end
    start == 0 && return (0, 0, 0)

    inquote = l[start] == '\'' || l[start] == '"'
    quotesymb = l[start]
    inmultiline = l[start] == ';' && l[prevind(l,start)] == '\n'
    instatus = inmultiline | inquote
    instatus && (start = nextind(l, start); i = start)
    while i <= n
        c = l[i]
        if !isspace(c)
            if instatus
                if inmultiline && c == ';' && l[prevind(l,i)] == '\n'
                    return (start, prevind(l, i-1), i)
                elseif inquote && c == quotesymb && i != n
                    nextl = l[nextind(l, i)]
                    if isspace(nextl) || nextl == '#'
                        return (start, prevind(l, i), i)
                    end
                end
            elseif c == '#' && isspace(l[prevind(l,i)])
                return (start, prevind(l, i), prevind(l, i))
            end
            i = nextind(l, i)
        elseif !instatus
            return (start, prevind(l, i), i)
        else
            i = nextind(l, i)
        end
    end
    if instatus
        if inmultiline
            error(lazy"Invalid syntax: opening multiline field at position $start is not closed")
        end
        if inquote
            error(lazy"Invalid syntax: opening quote $quotesymb at position $start is not closed")
        end
    end
    return (start, n, n)
end

function angle(p1, p2)
    return acosd(clamp(dot(p1, p2)/(norm(p1)*norm(p2)), -1.0, 1.0))
end
function dihedral(p1, p2, p3)
    β = angle(cross(p1, p2), cross(p2, p3))
    return isnan(β) ? 0.0 : β
end

function isinterrupt(@nospecialize(e))
    return e isa InterruptException || (e isa TaskFailedException && e.task.result isa InterruptException)
end
function isoverfloworinexact(@nospecialize(e))
    (e isa OverflowError || e isa InexactError) && return true
    if e isa TaskFailedException
        result = e.task.result
        return result isa OverflowError || result isa InexactError
    end
    return false
end


text:
ADCME.options.newton_raphson.verbose = true
Δt = 0.01
mmesh = Mesh(30, 30, 1/30)
bdnode = bcnode(mmesh)
bdnode = [bdnode; bdnode .+ mmesh.ndof]
nu = constant(0.0001ones(get_ngauss(mmesh)))
nodes = fem_nodes(mmesh)
u = @. sin(2π * nodes[:,1])
v = @. cos(2π * nodes[:,2])
u0 = [u;v]
u0[bdnode] .= 0.0
us = solve_burgers(u0, 10)

sess = Session(); init(sess)
U = run(sess, us)
figure(figsize=(10,5))
close("all")
subplot(121)
visualize_scalar_on_fem_points(U[end,1:mmesh.nnode], mmesh)
subplot(122)
visualize_scalar_on_fem_points(U[end,mmesh.ndof + 1:mmesh.ndof + mmesh.nnode], mmesh)
savefig("test_mfem.png")

close("all")
visualize_vector_on_fem_points(U[end,1:mmesh.nnode], U[end,mmesh.ndof + 1:mmesh.ndof + mmesh.nnode], mmesh)
savefig("mfem_quiver.png")
# function ff(x)
#     unext = constant(x[1:mmesh.ndof])
#     vnext = constant(x[mmesh.ndof+1:end])
#     r, J = calc_residual_and_jacobian(unext, vnext, u, v, mmesh)
#     run(sess, r), run(sess, J)
# end

# ff(rand(2mmesh.ndof))
# sess = Session(); init(sess)
# test Jacobian
# test_jacobian(ff, rand(2mmesh.ndof))


text:
module ConeProgramDiff

using LinearAlgebra, SparseArrays                       #stdlib
using JuMP, MathOptInterface, MathOptSetDistances       #Opt utils
using IterativeSolvers, Roots                           #Eq solvers
using SCS, ECOS, Hypatia                                #Opt solvers
using Random

const MOI = MathOptInterface
const MOSD = MathOptSetDistances

include("utils.jl")
include("cones.jl")
include("cone_solve.jl")
include("generate_program.jl")

export solve_and_diff
export project_onto_cone, d_project_onto_cone
export random_cone_program, l1_minimization_program

end

# TODO: Check out https://github.com/Jutho/LinearMaps.jl for projections


text:
repo = joinpath(homedir(), "ReferenceImages", "gallery")
recordings = joinpath(@__DIR__, "test_recordings")
cp(repo, recordings, force = true)
?


text:
# ---------------------------- #
# bitboard stuff on primitives #
# ---------------------------- #

encode(r::Integer, c::Integer) = UInt64(1) << (8*(r-1) + (c-1))
encode(rc::Tuple{Integer,Integer}) = encode(rc[1], rc[2])
const CENTERS_BITS = sum(encode(i) for i in [(4, 4), (5, 5), (5, 4), (4, 5)])

function _one_directon_right_shift(g::UInt64, p::UInt64, n::Integer)
    g |= p & (g >> (n * 1))
    p &= (p >> (n * 1))
    g |= p & (g >> (n * 2))
    p &= (p >> (n * 2))
    g |= p & (g >> (n * 4))
    g
end

function _one_directon_left_shift(g::UInt64, p::UInt64, n::Integer)
    g |= p & (g << (n * 1))
    p &= (p << (n * 1))
    g |= p & (g << (n * 2))
    p &= (p << (n * 2))
    g |= p & (g << (n * 4))
    g
end

function _get_g_directions(p1::UInt64, p2::UInt64)
    mask_a = 0xfefefefefefefefe
    mask_h = 0x7f7f7f7f7f7f7f7f

    gN = _one_directon_right_shift(p1, p2, 8)
    gS = _one_directon_left_shift(p1, p2, 8)
    gE = _one_directon_left_shift(p1, p2 & mask_a, 1)
    gW = _one_directon_right_shift(p1, p2 & mask_h, 1)
    gNE = _one_directon_right_shift(p1, p2 & mask_a, 7)
    gNW = _one_directon_right_shift(p1, p2 & mask_h, 9)
    gSE = _one_directon_left_shift(p1, p2 & mask_a, 9)
    gSW = _one_directon_left_shift(p1, p2 & mask_h, 7)

    return gN, gS, gE, gW, gNE, gNW, gSE, gSW
end

function legal_actions_bits(mover::UInt64, opp::UInt64)
    occupied = mover | opp
    if (CENTERS_BITS & occupied) != CENTERS_BITS
        return CENTERS_BITS - occupied
    end

    mask_a = 0xfefefefefefefefe
    mask_h = 0x7f7f7f7f7f7f7f7f
    empty = 0xffffffffffffffff ⊻ occupied

    gN, gS, gE, gW, gNE, gNW, gSE, gSW = _get_g_directions(mover, opp)

    legal = UInt64(0)
    legal |= ((gN & ~mover) >> 8) & empty
    legal |= ((gS & ~mover) << 8) & empty
    legal |= ((gE & ~mover & mask_h) << 1) & empty
    legal |= ((gW & ~mover & mask_a) >> 1) & empty
    legal |= ((gNE & ~mover & mask_h) >> 7) & empty
    legal |= ((gNW & ~mover & mask_a) >> 9) & empty
    legal |= ((gSE & ~mover & mask_h) << 9) & empty
    legal |= ((gSW & ~mover & mask_a) << 7) & empty
    legal
end

function next_state_bits(action::UInt64, mover::UInt64, opponent::UInt64)
    gN, gS, gE, gW, gNE, gNW, gSE, gSW = _get_g_directions(
        action, opponent
    )

    mask_a = 0xfefefefefefefefe
    mask_h = 0x7f7f7f7f7f7f7f7f
    mover += action

    flips = UInt64(0)
    ((gN >> 8) & mover > 0)           && (flips |= gN)
    ((gS << 8) & mover > 0)           && (flips |= gS)
    ((gE << 1) & mask_a & mover > 0)  && (flips |= gE)
    ((gW >> 1) & mask_h & mover > 0)  && (flips |= gW)
    ((gNE >> 7) & mask_a & mover > 0) && (flips |= gNE)
    ((gNW >> 9) & mask_h & mover > 0) && (flips |= gNW)
    ((gSE << 9) & mask_a & mover > 0) && (flips |= gSE)
    ((gSW << 7) & mask_h & mover > 0) && (flips |= gSW)

    new_mover = mover | flips
    new_opponent = opponent & (~flips)

    new_mover, new_opponent
end

score(x::UInt64) = count_ones(x)

function decode(x::UInt64)
    out = Array{Tuple{Int,Int}}(uninitialized, count_ones(x))
    i = 0
    v = UInt64(1)  # equal to encode(1, 1)
    for r in 1:8
        for c in 1:8
            if v & x > 0
                out[i += 1] = (r, c)
            end
            v <<= 1  # now equal to encode(r, c+1) if (c <= 7)
                     #           or encode(r+1, 1) if (c == 8) and (r < 8)
        end
    end
    out
end

legal_actions(p1::UInt64, p2::UInt64) = decode(legal_actions_bits(p1, p2))

function game_over(p1::UInt64, p2::UInt64)::Tuple{Bool,Int,Int}
    if p1 == 0
        return true, 0, score(p2)
    end

    if p2 == 0
        return true, score(p1), 0
    end

    occupied = p1 | p2
    if occupied == (UInt64(1) << (64)) - 1
        return true, score(p1), score(p2)
    end

    if isempty(legal_actions(p1, p2)) && isempty(legal_actions(p2, p1))
        return true, score(p1), score(p2)
    end

    false, 0, 0
end

# ---------- #
# Board type #
# ---------- #

struct Board
    p1_placed::UInt64
    p2_placed::UInt64
end

function Base.show(io::IO, ::MIME"text/plain", b::Board)
    # corners
    top_left = "\u250C"
    top_right = "\u2510"
    bottom_left = "\u2514"
    bottom_right = "\u2518"

    # lines
    hori = "\u2500"
    vert = "\u2502"

    # line + center piece
    hori_down = "\u252C"
    hori_up = "\u2534"

    print(io, " ", join(map(i->string("  ", i, " "), 1:8)), "\n")
    print(io, " ", top_left, hori, repeat(string(hori^2, hori_down, hori), 7), hori^2, top_right)

    for r in UInt64(1):8
        println(io)
        print(io, r)
        print(io, vert)
        for c in UInt64(1):8
            v = encode(r, c)
            if v & b.p1_placed > 0
                print(io, MARKERS[1])
            elseif v & b.p2_placed > 0
                print(io, MARKERS[2])
            else
                print(io, MARKERS[0])
            end
            print(io, vert)
        end
    end

    println(io, "\n ", bottom_left, hori, repeat(string(hori^2, hori_up, hori), 7), hori^2, bottom_right)
end
Base.show(io::IO, b::Board) = show(io, MIME"text/plain"(), b)


text:
module MaxHelpingHandKevinBarrier

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/KevinBarrier" KevinBarrier(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight,
    color::String="62222b", particleColor::String="ffffff", flashOnHit::Bool=true, invisible::Bool=false)

const placements = Ahorn.PlacementDict(
    "Kevin Barrier (max480's Helping Hand)" => Ahorn.EntityPlacement(
        KevinBarrier,
        "rectangle"
    )
)

Ahorn.minimumSize(entity::KevinBarrier) = 8, 8
Ahorn.resizable(entity::KevinBarrier) = true, true

function Ahorn.selection(entity::KevinBarrier)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return Ahorn.Rectangle(x, y, width, height)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::KevinBarrier, room::Maple.Room)
    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    color = Ahorn.argb32ToRGBATuple(parse(Int, get(entity.data, "color", "62222b"), base=16))[1:3] ./ 255
    Ahorn.drawRectangle(ctx, 0, 0, width, height, color, (0.0, 0.0, 0.0, 0.0))
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
mutable struct Node{T}
	m::JuMP.Model
	z::Float64
	x::Vector{Float64}
	zparent::Float64
	nbcover::Int
	f0::Set{Int}
	f1::Set{Int}
	bound1::Float64
	bound2::Float64
	cstr1::T
	cstr2::T
end

function Node(m, b1, b2, c1::T, c2::T) where T
	Node(m, NaN, zeros(m.numCols), Inf, 0, Set{Int}(), Set{Int}(), Float64(b1), Float64(b2), c1, c2)
end
# Base.copy(n::Node) = begin warn("deepcopy of a node") ; return Node(copy(n.m), NaN, n.x, n.z, n.nbcover, copy(n.f0), copy(n.f1), n.bound1, n.bound2, n.cstr1, n.cstr2) end
shallowcopy(n::Node) = Node(n.m, NaN, n.x, n.z, n.nbcover, copy(n.f0), copy(n.f1), n.bound1, n.bound2, n.cstr1, n.cstr2)

mutable struct NodeParragh{T}
	m::JuMP.Model
	x::Vector{Vector{Float64}}
	z_avg::Float64
	nbcover::Int
	f0::Set{Int}
	f1::Set{Int}
	bound1::Float64
	bound2::Float64
	cstr1::T
	cstr2::T
end

function NodeParragh(m, b1, b2, c1::T, c2::T) where T
	NodeParragh(m, Vector{Float64}[], NaN, 0, Set{Int}(), Set{Int}(), Float64(b1), Float64(b2), c1, c2)
end
# Base.copy(n::NodeParragh) = NodeParragh(copy(n.m), Vector{Float64}[], n.nbcover, copy(n.f0), copy(n.f1))
shallowcopy(n::NodeParragh) = NodeParragh(n.m, Vector{Float64}[], n.z_avg, n.nbcover, copy(n.f0), copy(n.f1), n.bound1, n.bound2, n.cstr1, n.cstr2)

Base.show(io::IO, n::NodeParragh) = print(io, "Node[f0=$(n.f0), f1=$(n.f1), <$(n.bound1)-$(n.bound2)>]")


text:
function convert{DA <: Array, FSA <: FixedArray}(::Type{DA}, b::FSA)
    result = Array(eltype(b), size(b)...)
    @inbounds for i=1:length(b)
        result[i] = b[i]
    end
    result
end

convert{R, C, R2, C2, T}(a::Type{Mat{R, C, T}}, b::Mat{R2, C2, T}) =
    Mat(ntuple(c -> ntuple(r -> b[r,c], Val{R}), Val{C}))

#conversion
convert{T <: Tuple}(::Type{T}, x::Real)                   = ntuple(ConstFunctor(eltype(T)(x)), Val{length(T.parameters),})
convert{T <: Tuple, FSA <: FixedArray}(::Type{T}, x::FSA) = map(eltype(T), x.(1)[1:length(T.parameters)])
convert{T <: FixedArray}(t::Type{T}, f::T)                = f
convert{FSA1 <: FixedArray}(t::Type{FSA1}, f::FixedArray) =
    map(ConversionIndexFunctor(f, eltype_or(FSA1, eltype(typeof(f)))), FSA1)

convert{T <: FixedArray}(::Type{T}, x) = T(x)
convert{T <: FixedArray}(::Type{T}, x...) = T(x)


text:
using MixedModeBroadcastAD: broadcast_gradients!, Wrt
using CUDAdrv, CUDAnative

sigm(x)      = @fastmath 1 / (1 + exp(-x))
cuda_sigm(x) = @fastmath 1 / (1 + CUDAnative.exp_fast(-x))
cuda_tanh(x) = CUDAnative.tanh(x)

###################
# kernel selector #
###################

@noinline broadcast_wrapper(f::F) where {F} = (inputs, derivs, buffers) -> broadcast_gradients!(f, inputs, derivs)

function initialize_inputs(::Type{A}, dims::Int, uniform::Bool) where {A<:AbstractArray}
    # set up control variables
    function random_control(dims)
        control = (round.(rand(Float32, dims)), round.(rand(Float32, dims)))

        # ensure that we hit all three cases in the HMLSTM update algorithm at least once
        @assert dims >= 3
        control[1][1] = 1.0f0 # FLUSH case
        control[2][1] = 0.0f0
        control[1][2] = 0.0f0 # COPY case
        control[2][2] = 0.0f0
        control[1][3] = 0.0f0 # UPDATE case
        control[2][3] = 1.0f0

        return control
    end
    function nondivergent_control(dims)
        # repeat sequences of warp uniform control values to avoid divergent execution
        warpsize = 32
        @assert rem(dims, warpsize) == 0
        subdim = Int(dims / warpsize)

        control = random_control(subdim)
        Tuple(collect(Iterators.flatten(fill(x, warpsize) for x in bools)) for bools in control)
    end
    control = (uniform ? nondivergent_control : random_control)(dims)

    control = (convert(A, control[1]), convert(A, control[2]))
    return (control..., (convert(A, rand(Float32, dims, dims)) for _ in 1:4)...,)
end

function get_hmlstm_kernel(tfstyle::Bool, usegpu::Bool, dims::Int = 2048, uniform::Bool = false)
    if usegpu
        scalar_kernel = gpu_hmlstm_update_c_scalar
        A = CuArray
    else
        scalar_kernel = cpu_hmlstm_update_c_scalar
        A = Array
    end
    inputs = initialize_inputs(A, dims, uniform)
    if tfstyle
        kernel! = tf_hmlstm_update_c_gradients!
        derivs = similar.(inputs[3:end])
        buffers = similar.(derivs[2:end])
    else
        kernel! = broadcast_wrapper(scalar_kernel)
        derivs = similar.(inputs[3:end])
        inputs = (inputs[1], inputs[2], Wrt.(inputs[3:end])...)
        buffers = ()
    end
    return kernel!, inputs, derivs, buffers
end

function get_arity_scaling_kernel(usegpu::Bool, dims::Int = 1024, arity::Int = 2)
    A = usegpu ? CuArray : Array
    kernel! = broadcast_wrapper(usegpu ? gpu_arity_scaling : cpu_arity_scaling)
    inputs = ((convert(A, rand(Float32, dims, dims)) for _ in 1:arity)...,)
    derivs = similar.(inputs)
    buffers = ()
    return kernel!, Wrt.(inputs), derivs, buffers
end

###############################################
# idiomatic Julia forward-pass scalar kernels #
###############################################

@fastmath function cpu_hmlstm_update_c_scalar(z, zb, c, f, i, g)
    if z == 1.0f0 # FLUSH
        return sigm(i) * tanh(g)
    elseif zb == 0.0f0 # COPY
        return c
    else # UPDATE
        return sigm(f) * c + sigm(i) * tanh(g)
    end
end

@fastmath function gpu_hmlstm_update_c_scalar(z, zb, c, f, i, g)
    if z == 1.0f0 # FLUSH
        return cuda_sigm(i) * cuda_tanh(g)
    elseif zb == 0.0f0 # COPY
        return c
    else # UPDATE
        return cuda_sigm(f) * c + cuda_sigm(i) * cuda_tanh(g)
    end
end

@generated function cpu_arity_scaling(args::Vararg{Any,N}) where {N}
    mapped = [:(args[$i] > 0.5f0 ? args[$i] : -args[$i]) for i in 1:N]
    ex = :(tanh($(pop!(mapped))))
    while !isempty(mapped)
        ex = Expr(:call, :*, ex, :(tanh($(pop!(mapped)))))
    end
    return ex
end

@generated function gpu_arity_scaling(args::Vararg{Any,N}) where {N}
    mapped = [:(args[$i] > 0.5f0 ? args[$i] : -args[$i]) for i in 1:N]
    ex = :(cuda_tanh($(pop!(mapped))))
    while !isempty(mapped)
        ex = Expr(:call, :*, ex, :(cuda_tanh($(pop!(mapped)))))
    end
    return ex
end

#########################################
# TF-style HMLSTM gradient calculations #
#########################################
# This code implements the computation described by the HLO graph and profile images found
# in this directory. The former was generated by running the TensorFlow code in `kernels.py`
# with the flag `TF_XLA_FLAGS=--xla_generate_hlo_graph=.*`, while the latter was generated
# by profiling the executed kernels using `nvprof`. Where reasonable, variable names used in
# these kernels match those used in the HLO graph.

function tf_hmlstm_update_c_gradients!(inputs::NTuple{6,AbstractArray},
                                       derivs::NTuple{4,AbstractArray},
                                       buffers::NTuple{3,AbstractArray})

    z, zb, c, f, i, g = inputs
    ∇c, ∇f, ∇i, ∇g = derivs
    P0, P1, P2, P3, P4, P5 = c, z, zb, f, g, i
    _tanh_func = ifelse(isa(first(inputs), CuArray), cuda_tanh, tanh)
    # these are executed in the same order as shown in `tf_hmlstm_profile.png`
    fusion2 = tf_fusion_2_or_5!(_tanh_func, buffers[1], P5) # sigm.(i)
    tanh1   = broadcast!(_tanh_func, buffers[2], P4) # tanh.(g)
    fusion5 = tf_fusion_2_or_5!(_tanh_func, buffers[3], P3) # sigm.(f)
    fusion  = tf_fusion!(∇g, fusion2, tanh1, P1, P2)
    fusion1 = tf_fusion1!(∇i, fusion2, tanh1, P1, P2)
    fusion3 = tf_fusion3!(∇f, fusion5, P0, P1, P2)
    fusion4 = tf_fusion4!(∇c, fusion5, P1, P2)
    return nothing
end

@fastmath function tf_fusion!(∇g, fusion2, tanh1, P1, P2)
    P5 = P1
    P4 = P2
    P3 = 0.0f0
    P2 = 1.0f0
    P1 = fusion2
    P0 = tanh1
    return broadcast!(∇g, P0, P1, P2, P3, P4, P5) do p0, p1, p2, p3, p4, p5
        equalto7 = p4 == p3
        equalto13 = p5 == p2
        select7 = ifelse(equalto13, p3, p2)
        select5 = ifelse(equalto13, p2, p3)
        select6 = ifelse(equalto7, p3, select7)
        multiply17 = p0 * p0
        multiply18 = select6 * p1
        multiply19 = select5 * p1
        subtract3 = p2 - multiply17
        add5 = multiply19 + multiply18
        return add5 * subtract3
    end
end

@fastmath function tf_fusion1!(∇i, fusion2, tanh1, P1, P2)
    P5 = P1
    P4 = P2
    P3 = 0.0f0
    P2 = 1.0f0
    P1 = tanh1
    P0 = fusion2
    return broadcast!(∇i, P0, P1, P2, P3, P4, P5) do p0, p1, p2, p3, p4, p5
        equalto9 = p3 == p4
        equalto15 = p2 == p5
        select8 = ifelse(equalto15, p2, p3)
        select10 = ifelse(equalto15, p3, p2)
        select9 = ifelse(equalto9, p3, select10)
        multiply22 = select9 * p1
        multiply23 = select8 * p1
        add6 = multiply22 + multiply23
        multiply21 = add6 * p0
        subtract4 = p2 - p0
        return multiply21 * subtract4
    end
end

@fastmath function tf_fusion3!(∇f, fusion5, P0, P1, P2)
    P5 = P1
    P4 = P2
    P3 = 0.0f0
    P2 = 1.0f0
    P1 = P0
    P0 = fusion5
    return broadcast!(∇f, P0, P1, P2, P3, P4, P5) do p0, p1, p2, p3, p4, p5
        equalto11 = p3 == p4
        equalto17 = p2 == p5
        select12 = ifelse(equalto17, p3, p2)
        select11 = ifelse(equalto11, p3, select12)
        multiply28 = select11 * p1
        multiply27 = multiply28 * p0
        subtract5 = p2 - p0
        return multiply27 * subtract5
    end
end

@fastmath function tf_fusion4!(∇c, fusion5, P1, P2)
    P4 = P1
    P3 = P2
    P2 = 0.0f0
    P1 = 1.0f0
    P0 = fusion5
    return broadcast!(∇c, P0, P1, P2, P3, P4) do p0, p1, p2, p3, p4
        equalto5 = p3 == p2
        equalto19 = p4 == p1
        select14 = ifelse(equalto19, p2, p1)
        select4 = ifelse(equalto5, select14, p2)
        select13 = ifelse(equalto5, p2, select14)
        multiply29 = select13 * p0
        return select4 + multiply29
    end
end

# fusion.2 and fusion.5 are exactly the same,
# so we just use this method for both kernels
@fastmath function tf_fusion_2_or_5!(_tanh_func, output, P_3_or_5)
    P1 = 0.5f0
    P0 = P_3_or_5
    return broadcast!(output, P0, P1) do p0, p1
        return p1 + (p1 * _tanh_func(p1 * p0))
    end
end


text:
using LinearAlgebra

function RLS(input_trn, output_trn, input_tst, output_tst, ϕ; M1=1, M2=1, M3=1, N=1, λ=1.00)

    # Maximum delay
    maxM = maximum([M1,M2,M3])

    # Zero-padding of signals
    input_trn = [zeros(maxM,); input_trn]
    input_tst = [zeros(maxM,); input_tst]
    output_trn = [zeros(maxM,); output_trn]
    output_tst = [zeros(maxM,); output_tst]

    # Signal splits
    T_trn = length(input_trn)
    T_tst = length(input_tst)

    # Parameters
    P = λ.*Matrix{Float64}(I,N,N)
    w_k = zeros(N,)

    # Preallocate prediction array
    predictions = zeros(T_trn,)
    errors = zeros(T_trn,)

    for (ii,k) in enumerate(maxM+1:T_trn)

        # Update history vectors
        u_kmin1 = input_trn[k-1:-1:k-M1]
        y_kmin1 = output_trn[k-1:-1:k-M2]
        e_kmin1 = errors[k-1:-1:k-M3]

        # Update data vector
        ϕx = ϕ([input_trn[k]; u_kmin1; y_kmin1; e_kmin1])

        # Update weights
        α = output_trn[k] - w_k'*ϕx
        g = P*ϕx*inv(λ + ϕx'*P*ϕx)
        P = inv(λ)*P - g*ϕx'*inv(λ)*P
        w_k = w_k + α*g

        # Prediction
        predictions[k] = w_k'*ϕx
        errors[k] = output_trn[k] - predictions[k]

    end

    "1-step ahead prediction"

    # Prepare array
    predictions = zeros(T_tst,)
    errors = zeros(T_tst,)

    for k in maxM+1:T_tst

        # Update history vectors
        u_kmin1 = input_tst[k-1:-1:k-M1]
        y_kmin1 = output_tst[k-1:-1:k-M2]
        e_kmin1 = errors[k-1:-1:k-M3]

        # Posterior predictive
        ϕx = ϕ([input_tst[k]; u_kmin1; y_kmin1; e_kmin1])
        predictions[k] = w_k'*ϕx

        # Update error
        errors[k] = output_tst[k] - predictions[k]

    end

    "Simulation"

    # Prepare array
    simulations = zeros(T_tst,)

    for k in maxM+1:T_tst

        # Update history vectors
        u_kmin1 = input_tst[k-1:-1:k-M1]
        y_kmin1 = simulations[k-1:-1:k-M2]
        e_kmin1 = zeros(M3,)

        # Posterior predictive
        ϕx = ϕ([input_tst[k]; u_kmin1; y_kmin1; e_kmin1])
        simulations[k] = w_k'*ϕx

    end

    "Evaluation"

    # Compute root mean square error
    RMS_prd = sqrt(mean((predictions[maxM+1:end] - output_tst[maxM+1:end]).^2))
    RMS_sim = sqrt(mean((simulations[maxM+1:end] - output_tst[maxM+1:end]).^2))

    return RMS_sim, RMS_prd
end

function run_system(input_tst, output_tst, ϕ, θ; M1=1, M2=1, M3=1, N=1)

    # Maximum delay
    maxM = maximum([M1,M2,M3])

    # Zero-padding of signals
    input_tst = [zeros(maxM,); input_tst]
    output_tst = [zeros(maxM,); output_tst]

    # Length of validation signal
    T_tst = length(input_tst)

    "1-step ahead prediction"

    # Prepare array
    predictions = zeros(T_tst,)
    errors = zeros(T_tst,)

    for k in maxM+1:T_tst

        # Update history vectors
        u_kmin1 = input_tst[k-1:-1:k-M1]
        y_kmin1 = output_tst[k-1:-1:k-M2]
        e_kmin1 = errors[k-1:-1:k-M3]

        # Posterior predictive
        ϕx = ϕ([input_tst[k]; u_kmin1; y_kmin1; e_kmin1])
        predictions[k] = θ'*ϕx

        # Update error
        errors[k] = output_tst[k] - predictions[k]

    end

    "Simulation"

    # Prepare array
    simulations = zeros(T_tst,)

    for k in maxM+1:T_tst

        # Update history vectors
        u_kmin1 = input_tst[k-1:-1:k-M1]
        y_kmin1 = simulations[k-1:-1:k-M2]
        e_kmin1 = zeros(M3,)

        # Posterior predictive
        ϕx = ϕ([input_tst[k]; u_kmin1; y_kmin1; e_kmin1])
        simulations[k] = θ'*ϕx

    end

    "Evaluation"

    # Compute root mean square error
    RMS_prd = sqrt(mean((predictions[maxM+1:end] - output_tst[maxM+1:end]).^2))
    RMS_sim = sqrt(mean((simulations[maxM+1:end] - output_tst[maxM+1:end]).^2))

    return RMS_sim, RMS_prd
end


text:
using OREnvironment
using ORInterface
using Test

@testset "Building Status" begin
    numConstraints = 5;
    status = OREnvironment.constructStatus(numConstraints);
    # conditions after creation
    @test OREnvironment.is_feasible(status) == false;
    @test OREnvironment.is_optimal(status) == false;
    @test OREnvironment.get_objfunction(status) == 0.0;
    for i in 1:numConstraints
        @test OREnvironment.get_constraint_consumption(status, i) == 0.0;
    end
    # conditions after modification
    OREnvironment.set_feasible!(status, true);
    OREnvironment.set_optimal!(status, true);
    OREnvironment.set_objfunction!(status, 16.0);
    OREnvironment.set_constraint_consumption!(status, 1, 15.0);
    @test OREnvironment.is_feasible(status) == true;
    @test OREnvironment.is_optimal(status) == true;
    @test OREnvironment.get_objfunction(status) == 16.0;
    @test OREnvironment.get_constraint_consumption(status, 1) == 15.0;
end

@testset "is_first_status_better" begin
    numConstraints = 5;
    s1 = OREnvironment.constructStatus(numConstraints);
    s2 = OREnvironment.constructStatus(numConstraints);
    OREnvironment.set_objfunction!(s1, 34.0);
    OREnvironment.set_objfunction!(s2, 3.0);

    @testset "feasibility not required" begin
        feasibility = false;
        @test OREnvironment.is_first_status_better(s1,s2,:max, feasibility) == true;
        @test OREnvironment.is_first_status_better(s1,s2,:min, feasibility) == false;
        OREnvironment.set_objfunction!(s2, 34.0);
        @test OREnvironment.is_first_status_better(s1,s2,:max, feasibility) == false;
        @test OREnvironment.is_first_status_better(s2,s1,:max, feasibility) == false;
        @test OREnvironment.is_first_status_better(s1,s2,:min, feasibility) == false;
        @test OREnvironment.is_first_status_better(s2,s1,:min, feasibility) == false;

        OREnvironment.set_objfunction!(s1, 34.0);
        OREnvironment.set_objfunction!(s2, 3.0);
        OREnvironment.set_feasible!(s1, true);
        OREnvironment.set_feasible!(s2, false);
        @test OREnvironment.is_first_status_better(s1,s2,:max, feasibility) == true;
        @test OREnvironment.is_first_status_better(s1,s2,:min, feasibility) == false;
        OREnvironment.set_objfunction!(s2, 34.0);
        @test OREnvironment.is_first_status_better(s1,s2,:max, feasibility) == false;
        @test OREnvironment.is_first_status_better(s2,s1,:max, feasibility) == false;
        @test OREnvironment.is_first_status_better(s1,s2,:min, feasibility) == false;
        @test OREnvironment.is_first_status_better(s2,s1,:min, feasibility) == false;

        OREnvironment.set_objfunction!(s1, 34.0);
        OREnvironment.set_objfunction!(s2, 3.0);
        OREnvironment.set_feasible!(s1, false);
        OREnvironment.set_feasible!(s2, true);
        @test OREnvironment.is_first_status_better(s1,s2,:max, feasibility) == true;
        @test OREnvironment.is_first_status_better(s1,s2,:min, feasibility) == false;
        OREnvironment.set_objfunction!(s2, 34.0);
        @test OREnvironment.is_first_status_better(s1,s2,:max, feasibility) == false;
        @test OREnvironment.is_first_status_better(s2,s1,:max, feasibility) == false;
        @test OREnvironment.is_first_status_better(s1,s2,:min, feasibility) == false;
        @test OREnvironment.is_first_status_better(s2,s1,:min, feasibility) == false;

        OREnvironment.set_objfunction!(s1, 34.0);
        OREnvironment.set_objfunction!(s2, 3.0);
        OREnvironment.set_feasible!(s1, true);
        OREnvironment.set_feasible!(s2, true);
        @test OREnvironment.is_first_status_better(s1,s2,:max, feasibility) == true;
        @test OREnvironment.is_first_status_better(s1,s2,:min, feasibility) == false;
        OREnvironment.set_objfunction!(s2, 34.0);
        @test OREnvironment.is_first_status_better(s1,s2,:max, feasibility) == false;
        @test OREnvironment.is_first_status_better(s2,s1,:max, feasibility) == false;
        @test OREnvironment.is_first_status_better(s1,s2,:min, feasibility) == false;
        @test OREnvironment.is_first_status_better(s2,s1,:min, feasibility) == false;
    end

    @testset "feasibility required" begin
        feasibility = true;
        OREnvironment.set_objfunction!(s1, 34.0);
        OREnvironment.set_objfunction!(s2, 3.0)
        OREnvironment.set_feasible!(s1, false);
        OREnvironment.set_feasible!(s2, false);
        @test OREnvironment.is_first_status_better(s1,s2,:max, feasibility) == true;
        @test OREnvironment.is_first_status_better(s1,s2,:min, feasibility) == false;
        OREnvironment.set_objfunction!(s2, 34.0);
        @test OREnvironment.is_first_status_better(s1,s2,:max, feasibility) == false;
        @test OREnvironment.is_first_status_better(s2,s1,:max, feasibility) == false;
        @test OREnvironment.is_first_status_better(s1,s2,:min, feasibility) == false;
        @test OREnvironment.is_first_status_better(s2,s1,:min, feasibility) == false;

        OREnvironment.set_objfunction!(s1, 34.0);
        OREnvironment.set_objfunction!(s2, 3.0);
        OREnvironment.set_feasible!(s1, true);
        OREnvironment.set_feasible!(s2, false);
        @test OREnvironment.is_first_status_better(s1,s2,:max, feasibility) == true;
        @test OREnvironment.is_first_status_better(s1,s2,:min, feasibility) == true;
        OREnvironment.set_objfunction!(s2, 34.0);
        @test OREnvironment.is_first_status_better(s1,s2,:max, feasibility) == true;
        @test OREnvironment.is_first_status_better(s2,s1,:max, feasibility) == false;
        @test OREnvironment.is_first_status_better(s1,s2,:min, feasibility) == true;
        @test OREnvironment.is_first_status_better(s2,s1,:min, feasibility) == false;

        OREnvironment.set_objfunction!(s1, 34.0);
        OREnvironment.set_objfunction!(s2, 3.0);
        OREnvironment.set_feasible!(s1, false);
        OREnvironment.set_feasible!(s2, true);
        @test OREnvironment.is_first_status_better(s1,s2,:max, feasibility) == false;
        @test OREnvironment.is_first_status_better(s1,s2,:min, feasibility) == false;
        OREnvironment.set_objfunction!(s2, 34.0);
        @test OREnvironment.is_first_status_better(s1,s2,:max, feasibility) == false;
        @test OREnvironment.is_first_status_better(s2,s1,:max, feasibility) == true;
        @test OREnvironment.is_first_status_better(s1,s2,:min, feasibility) == false;
        @test OREnvironment.is_first_status_better(s2,s1,:min, feasibility) == true;

        OREnvironment.set_objfunction!(s1, 34.0);
        OREnvironment.set_objfunction!(s2, 3.0);
        OREnvironment.set_feasible!(s1, true);
        OREnvironment.set_feasible!(s2, true);
        @test OREnvironment.is_first_status_better(s1,s2,:max, feasibility) == true;
        @test OREnvironment.is_first_status_better(s1,s2,:min, feasibility) == false;
        OREnvironment.set_objfunction!(s2, 34.0);
        @test OREnvironment.is_first_status_better(s1,s2,:max, feasibility) == false;
        @test OREnvironment.is_first_status_better(s2,s1,:max, feasibility) == false;
        @test OREnvironment.is_first_status_better(s1,s2,:min, feasibility) == false;
        @test OREnvironment.is_first_status_better(s2,s1,:min, feasibility) == false;
    end
end

@testset "Worst Value" begin
    @test OREnvironment.worst_value(:max) == -Inf;
    @test OREnvironment.worst_value(:min) == Inf;
end

@testset "update_status" begin
    numConstraints = 5;
    # first status
    s1 = OREnvironment.constructStatus(numConstraints);
    OREnvironment.set_objfunction!(s1, 34.0);
    OREnvironment.set_feasible!(s1, true);
    OREnvironment.set_optimal!(s1, true);
    OREnvironment.set_constraint_consumption!(s1, 1, 1.2);
    OREnvironment.set_constraint_consumption!(s1, 2, 2.2);
    OREnvironment.set_constraint_consumption!(s1, 3, 3.2);
    OREnvironment.set_constraint_consumption!(s1, 4, 4.2);
    OREnvironment.set_constraint_consumption!(s1, 5, 5.2);

    @test OREnvironment.is_feasible(s1) == true;
    @test OREnvironment.is_optimal(s1) == true;
    @test OREnvironment.get_objfunction(s1) == 34.0;
    @test OREnvironment.get_constraint_consumption(s1, 1) == 1.2;
    @test OREnvironment.get_constraint_consumption(s1, 2) == 2.2;
    @test OREnvironment.get_constraint_consumption(s1, 3) == 3.2;
    @test OREnvironment.get_constraint_consumption(s1, 4) == 4.2;
    @test OREnvironment.get_constraint_consumption(s1, 5) == 5.2;

    # second status
    s2 = OREnvironment.constructStatus(numConstraints);
    @test OREnvironment.is_feasible(s2) == false;
    @test OREnvironment.is_optimal(s2) == false;
    @test OREnvironment.get_objfunction(s2) == 0.0;
    for i in 1:numConstraints
        @test OREnvironment.get_constraint_consumption(s2, i) == 0.0;
    end

    # update
    OREnvironment.update_status!(s2, s1);

    @test OREnvironment.is_feasible(s1) == true;
    @test OREnvironment.is_optimal(s1) == true;
    @test OREnvironment.get_objfunction(s1) == 34.0;
    @test OREnvironment.get_constraint_consumption(s1, 1) == 1.2;
    @test OREnvironment.get_constraint_consumption(s1, 2) == 2.2;
    @test OREnvironment.get_constraint_consumption(s1, 3) == 3.2;
    @test OREnvironment.get_constraint_consumption(s1, 4) == 4.2;
    @test OREnvironment.get_constraint_consumption(s1, 5) == 5.2;

    @test OREnvironment.is_feasible(s2) == true;
    @test OREnvironment.is_optimal(s2) == true;
    @test OREnvironment.get_objfunction(s2) == 34.0;
    @test OREnvironment.get_constraint_consumption(s2, 1) == 1.2;
    @test OREnvironment.get_constraint_consumption(s2, 2) == 2.2;
    @test OREnvironment.get_constraint_consumption(s2, 3) == 3.2;
    @test OREnvironment.get_constraint_consumption(s2, 4) == 4.2;
    @test OREnvironment.get_constraint_consumption(s2, 5) == 5.2;
end


text:
struct ArbitraryMotion
    # e.g. ptpJointSpace!
    func::Function
    inputArg1       # e.g. positions, waitingPeriod, robotOrDepot
    inputArg2       # e.g. nothing, movablePart
    inputArg3       # e.g. nothing, waitingPeriod
end

mutable struct ReferencePath
    names::Vector{AbstractString}     # names are set by the user
    initPos::Vector{Float64}          # initial positions
    position::Vector{Float64}         # current reference position
    velocity::Vector{Float64}         # current reference velocity
    acceleration::Vector{Float64}     # current reference acceleration
    ptpPath::Union{PTP_path,Nothing}  # current PTP_path definition
    nextEventTime::Float64            # time at which the next event is triggered
    v_max::Vector{Float64}            # max velocity is set by the user
    a_max::Vector{Float64}            # max acceleration is set by the user
    refMotion::Vector{ArbitraryMotion} # stores explicit robot functions

    # both are updated in calculateRobotMovement
    instantiatedModel::ModiaLang.SimulationModel     # SimulationModel of model
    # scene::Composition.Scene                       # scene

    function ReferencePath(names::AbstractVector;
                position = zeros(size(names,1)),
                v_max    = ones(size(names,1)),
                a_max    = ones(size(names,2)))
        new(names, position, position, zeros(size(names,1)), zeros(size(names,1)), nothing, 0.001, v_max, a_max, ArbitraryMotion[])
    end
end

ReferencePath(; names, position, v_max, a_max) = ReferencePath(names, position=position, v_max=v_max, a_max=a_max)

getRefPathPosition(referencePath, index) = referencePath.position[index]
getRefPathPosition(; referencePath, index) = referencePath.position[index]

getRefPathInitPosition(referencePath, index) = referencePath.initPos[index]
getRefPathInitPosition(; referencePath, index) = referencePath.initPos[index]

## -------------------- utility functions ------------------------
function setAttributesReferencePath!(ref::ReferencePath, waitingPeriod::Float64)
    ref.nextEventTime = waitingPeriod
    ModiaLang.setNextEvent!(ref.instantiatedModel, ref.nextEventTime)
end

# getComponent(referencePath::ReferencePath, path::String) =
#    getPathComponent(referencePath.simulationState.model.model.assembly, path)

#--------------  Functions that can be executed at an event instant -------
### ------------------ ptpJointSpace! -------------------------------------
function ptpJointSpace!(ref::ReferencePath, positions::Matrix{Float64}, dummy1, dummy2)
    pos1 = size(transpose(ref.position))
    pos2 = size(positions)
    if pos1[2] != pos2[2]
        error("... from ptpJointSpace: Please check out dimensions of input $positions (= $pos2). Its dimensions should be (n, $(pos1[2])), n > 0.")
    end
    ref.ptpPath = PTP_path(ref.names,
    positions = [transpose(ref.position); positions],
    startTime = ModiaLang.getTime(ref.instantiatedModel),
    v_max     = ref.v_max,
    a_max     = ref.a_max)
    setAttributesReferencePath!(ref, pathEndTime(ref.ptpPath))
end

ptpJointSpace!(ref::ReferencePath, positions::Vector{Float64}, dummy1, dummy2) =
    ptpJointSpace!(ref, reshape(positions, :, 1), dummy1, dummy2)

ptpJointSpace!(ref::ReferencePath, positions::Float64, dummy1, dummy2) =
    ptpJointSpace!(ref, reshape([positions], :, 1), dummy1, dummy2)