text:
function phase_run(j)
  info("initiating run phase")
  r = [@spawnat w runfun(j) for w in workers()]
  map(wait, r)
  info("run phase finished")
end

function runfun(j)
  info("runfun worker $(myid()) initiated")
  length(j.inputs) > 2 && throw(ArgumentError("Cannot operate on more than two inputs."))
  length(j.inputs) == 1 && runfunsinglearg(j)
  length(j.inputs) == 2 && runfunmultiplearg(j)
  info("runfun worker $(myid()) finished")
end

function runfunsinglearg(j)
  info("runfunsinglearg $(myid()) initiated")
  input_filename = j.inputs[1]
  inputfile = joinpath(splitdir, string(myid()), input_filename)
  io_input = open(inputfile)
  io_output = write_sink(j)
  for (linenum, line) in enumerate(eachline(io_input))
    try
      v = j.fun(line)
      v == "" && continue
      write(io_output, json(v) * "\n")
    catch e
      if typeof(e) in [MapperException, UDFException]
        warn("[JOB $(j.jobid)] mapper encountered an exception of type $(typeof(e)) on line $(linenum): $(e.msg)")
      else
        rethrow(e)
      end
    end
  end

  if length(ctx) > 0
    for (k, v) in ctx
      write(io_output, json((k, v)) * "\n")
    end
  end
  close(io_input)
  close(io_output)
  info("runfunsinglearg $(myid()) finished")
end

function runfunmultiplearg(j)
  info("runfunmultiplearg $(myid()) initiated")
  input_filenamex = j.inputs[1]
  input_filenamey = j.inputs[2]
  io_output = write_sink(j)
  inputfilex = joinpath(splitdir, string(myid()), input_filenamex)
  for worker in workers()
    inputfiley = joinpath(splitdir, string(worker), input_filenamey)
    io_inputy = open(inputfiley)
    io_inputx = open(inputfilex)
    for (lineynum, liney) in enumerate(eachline(io_inputy))
      for (linexnum, linex) in enumerate(eachline(io_inputx))
        try
          v = j.fun(linex, liney)
          v == "" && continue
          write(io_output, json(v) * "\n")
        catch e
          if typeof(e) in [MapperException, UDFException]
            warn("[JOB $(j.jobid)] mapper encountered an exception of type $(typeof(e)) on line $(linenum): $(e.msg)")
          else
            rethrow(e)
          end
        end
      end
    end
    close(io_inputy)
    close(io_inputx)
  end
  if length(ctx) > 0
    for (k, v) in ctx
      write(io_output, json((k, v)) * "\n")
    end
  end
  close(io_output)
  info("runfunmultiplearg $(myid()) finished")
end

function runcombiner(filename, fun_combiner)
  r = [@spawnat w readlines(joinpath(splitdir, string(w), filename)) for w in workers()]
  combined = []
  map(wait, r)
  for p in r
    append!(combined, fetch(p))
  end
  fun_combiner(combined)
end

max_stars_repo_name: deyandyankov/nmr

max_issues_repo_name: deyandyankov/nmr

max_forks_repo_name: deyandyankov/nmr


text:
@inline function check_index(dict::CardinalDict{V,K}, idx::K) where {K<:Signed, V}
    return one(K) <= idx <= dict.idxmax
end
@inline function check_index(dict::CardinalDict{V,K}, idx::J) where {J<:Signed,K,V}
    return check_index(dict, idx%K)
end

@inline function guarded_index(dict::CardinalDict{V,K}, idx::K) where {K<:Signed, V}
    check_index(dict, idx) || throw(ErrorException("invalid index ($(idx))"))
    return nothing
end
@inline function guarded_index(dict::CardinalDict{V,K}, idx::J) where {J<:Signed,K,V}
    return guarded_index(dict, idx%K)
end

@inline function check_present(dict::CardinalDict{V,K}, idx::K) where {K<:Signed, V}
    return getindex(dict.present, idx)
end
@inline function check_present(dict::CardinalDict{V,K}, idx::J) where {J<:Signed,K,V}
    return check_present(dict, idx%K)
end

@inline function guarded_present(dict::CardinalDict{V,K}, idx::K) where {K<:Signed, V}
    check_present(dict, idx) || throw(ErrorException("there is no content at index ($(idx))"))
    return nothing
end
@inline function guarded_present(dict::CardinalDict{V,K}, idx::J) where {J<:Signed,K,V}
    return guarded_present(dict, idx%K)
end

@inline function keyval(dict::CardinalDict{V,K}, idx::K) where {K<:Signed, V}
    guarded_index(dict, idx)
    return idx
end
@inline function keyval(dict::CardinalDict{V,K}, idx::J) where {J<:Signed,K,V}
    return keyval(dict, idx%K)
end

@inline function unsafe_keyval(dict::CardinalDict{V,K}, idx::K) where {K<:Signed, V}
    return idx
end
@inline function unsafe_keyval(dict::CardinalDict{V,K}, idx::J) where {J<:Signed,K,V}
    return idx%K
end

# some essential Julia associative, dict API functions

function Base.haskey(dict::CardinalDict{V,K}, idx::K) where {K<:Signed, V}
    guarded_index(dict, idx)
    return check_present(dict, idx)
end
@inline function Base.haskey(dict::CardinalDict{V,K}, idx::J) where {J<:Signed, K<:Signed, V}
    return haskey(dict, idx%K)
end

# this is quite unsafe
function unsafe_haskey(dict::CardinalDict{V,K}, idx::K) where {K<:Signed, V}
    return check_present(dict, idx)
end
@inline function unsafe_haskey(dict::CardinalDict{V,K}, idx::J) where {J<:Signed, K<:Signed, V}
    return unsafe_haskey(dict, idx%K)
end

function Base.getindex(dict::CardinalDict{V,K}, idx::K) where {K<:Signed, V}
    guarded_index(dict, idx)
    guarded_present(dict, idx)
    return getindex(dict.content, idx)
end
@inline function Base.getindex(dict::CardinalDict{V,K}, idx::J) where {J<:Signed, K<:Signed, V}
    return getindex(dict, idx%K)
end

function unsafe_getindex(dict::CardinalDict{V,K}, idx::K) where {K<:Signed, V}
    guarded_present(dict, idx)
    return getindex(dict.content, idx)
end
@inline function unsafe_getindex(dict::CardinalDict{V,K}, idx::J) where {J<:Signed, K<:Signed, V}
    return unsafe_getindex(dict, idx%K)
end

function Base.get(dict::CardinalDict{V,K}, idx::K, default::V) where {K<:Signed, V}
    guarded_index(dict, idx)
    return check_present(dict, idx) ? getindex(dict.content, idx) : default
end
@inline function Base.get(dict::CardinalDict{V,K}, idx::J, default::V) where {J<:Signed, K<:Signed, V}
    return get(dict, idx%K, default)
end

function Base.get!(dict::CardinalDict{V,K}, idx::K, default::V) where {K<:Signed, V}
    guarded_index(dict, idx)
    result = default
    if check_present(dict, idx)
       result = getindex(dict, idx)
    else
       setindex!(dict, default, idx)
    end
    return result
end
@inline function Base.get!(dict::CardinalDict{V,K}, idx::J, default::V) where {J<:Signed, K<:Signed, V}
    return get!(dict, idx%K, default)
end

function Base.setindex!(dict::CardinalDict{V,K}, value::V, idx::K) where {K<:Signed, V}
     guarded_index(dict, idx)
     setindex!(dict.present, true, idx)
     setindex!(dict.content, value, idx)
     return value
end
@inline function Base.setindex!(dict::CardinalDict{V,K}, value::V, idx::J) where {J<:Signed, K<:Signed, V}
    return setindex!(dict, value, idx%K)
end

function unsafe_setindex!(dict::CardinalDict{V,K}, value::V, idx::K) where {K<:Signed, V}
     setindex!(dict.present, true, idx)
     setindex!(dict.content, value, idx)
     return value
end
@inline function unsafe_setindex!(dict::CardinalDict{V,K}, value::V, idx::J) where {J<:Signed, K<:Signed, V}
    return unsafe_setindex!(dict, value, idx%K)
end

ext: jl


text:
struct Barrier

    # surface defined as points where func(x,y) = val
    # rmap reflects points across the surface
    func::Function
    val::Float64
    rmap::Function
    orientation::Int64

end

function smoothstep(x, y, delta, bar::Barrier)

    u = bar.orientation*(bar.func(x,y) - bar.val) / delta + 0.5

    if u < 0
        return 0.0
    elseif 0 < u < 1
        return 3*u^2 - 2*u^3
    else
        return 1.0
    end
end

function reflection(x, y, bar::Barrier)

    if smoothstep(x, y, 0.0001, bar) < 0.5
        return bar.rmap(x, y)
    else
        return x, y
    end
end

function check_if_inside(x, y, deltas::Vector{Float64}, bars::Vector{Barrier})

    xy_is_inside = true

    for l=1:length(bars)
        if smoothstep(x, y, deltas[l], bars[l]) == 0
            xy_is_inside = false
            break
        end
    end

    return xy_is_inside
end

#convenience Barriers that will include all points or exclude all points
function Everywhere()

    func(x, y) = 0
    val = -Inf
    rmap(x, y) = x, y
    return Barrier(func, val, rmap, 1)
end

function Nowhere()

    func(x, y) = 0
    val = Inf
    rmap(x, y) = x, y
    return Barrier(func, val, rmap, 1)
end

# functions for geometric operations with polygons
function xy_to_uv(x, y, r1, r2)

    xs = x - r1[1]
    ys = y - r1[2]
    t = atan(r2[2] - r1[2], r2[1] - r1[1])

    u = xs * cos(t) + ys * sin(t)
    v = -xs * sin(t) + ys * cos(t)

    return Float64[u, v]
end

function uv_to_xy(u, v, r1, r2)

    t = atan(r2[2] - r1[2], r2[1] - r1[1])
    x = u * cos(t) - v * sin(t)
    y = u * sin(t) + v * cos(t)

    return r1 + Float64[x, y]
end

function poly_dist(x, y, vert)

    N = size(vert, 2) - 1
    h = zeros(Float64, N)

    for i=1:N
        r1 = vert[:,i]
        r2 = vert[:,i+1]
        u, v = xy_to_uv(x, y, r1, r2)
        if u < 0
            h[i] = norm([x,y] - r1) * sign(v)
        elseif 0 <= u < norm(r2 - r1)
            h[i] = v
        else
            h[i] = norm([x,y] - r2) * sign(v)
        end
    end

    h2, k = findmin(h .^ 2)
    return h[k], k
end

function ray_reflection(x, y, r1, r2)

    u, v = xy_to_uv(x, y, r1, r2)
    return uv_to_xy(u, -v, r1, r2)
end

# Barrier from polygon of given vertices
function PolyBarrier(vert)

    func(x,y) = poly_dist(x, y, vert)[1]

    function rmap(x, y)

        h, ns = poly_dist(x, y, vert)
        r1 = vert[:,ns]
        r2 = vert[:,ns+1]
        return ray_reflection(x, y, r1, r2)
    end

    return Barrier(func, 0.0, rmap, 1)
end

# constructors for cases where the "do-block" syntax might make the code more readable
Barrier(f::Function) = Barrier(f()...)
PolyBarrier(f::Function) = PolyBarrier(f())

# might be convenient to generate a polygon from parameterized curve
PolyBarrier(f::Function, N::Int64) = PolyBarrier(hcat(f.(LinRange(0, 1, N))))


text:
function LinearAlgebra.mul!(
    C::Matrix{T},
    blkA::BlockedSparse{T},
    adjB::Adjoint{T,<:BlockedSparse{T}},
    α::Number,
    β::Number,
) where {T}
    A = blkA.cscmat
    B = adjB.parent.cscmat
    B.m == size(C, 2) && A.m == size(C, 1) && A.n == B.n || throw(DimensionMismatch(""))
    anz = nonzeros(A)
    arv = rowvals(A)
    bnz = nonzeros(B)
    brv = rowvals(B)
    isone(β) || rmul!(C, β)
    @inbounds for j = 1:A.n
        for ib in nzrange(B, j)
            αbnz = α * bnz[ib]
            jj = brv[ib]
            for ia in nzrange(A, j)
                C[arv[ia], jj] += anz[ia] * αbnz
            end
        end
    end
    C
end

LinearAlgebra.mul!(
    C::Matrix{T},
    A::BlockedSparse{T},
    adjB::Adjoint{T,<:Matrix{T}},
    α::Number,
    β::Number,
) where {T} = mul!(C, A.cscmat, adjB, α, β)

LinearAlgebra.mul!(
    C::BlockedSparse{T},
    A::BlockedSparse{T},
    adjB::Adjoint{T,<:BlockedSparse{T}},
    α::Number,
    β::Number,
) where {T} = mul!(C.cscmat, A.cscmat, adjB.parent.cscmat', α, β)

LinearAlgebra.mul!(
    C::StridedVecOrMat{T},
    A::StridedVecOrMat{T},
    adjB::Adjoint{T,<:BlockedSparse{T}},
    α::Number,
    β::Number,
) where {T} = mul!(C, A, adjB.parent.cscmat', α, β)

LinearAlgebra.mul!(
    C::StridedVector{T},
    adjA::Adjoint{T,<:BlockedSparse{T}},
    B::StridedVector{T},
    α::Number,
    β::Number,
) where {T} = mul!(C, adjA.parent.cscmat', B, α, β)

function LinearAlgebra.ldiv!(
    adjA::Adjoint{T,<:LowerTriangular{T,UniformBlockDiagonal{T}}},
    B::StridedVector{T},
) where {T}
    A = adjA.parent
    length(B) == size(A, 2) || throw(DimensionMismatch(""))
    m, n, k = size(A.data.data)
    fv = A.data.facevec
    bb = reshape(B, (n, k))
    for j = 1:k
        ldiv!(adjoint(LowerTriangular(fv[j])), view(bb, :, j))
    end
    B
end

"""
    _rdivlowertriblk!(A::AbstractMatrix, Bd::Array, b, S)

Divide in-place a block of `S` columns of `A` starting at `(b-1)*s+1` on the right by the lower-triangular `b`th face of Bd.
"""
function _rdivlowertriblk!(A::AbstractMatrix, Bd::Array, b, S)
    coloffset = (b - 1) * S
    @inbounds @simd for i = 1:size(A, 1)
        for j = 1:S
            Aij = A[i, j + coloffset]
            for k = 1:j-1
                Aij -= A[i, k + coloffset] * Bd[j, k, b]'
            end
            A[i, j + coloffset] = Aij / Bd[j, j, b]'
        end
    end
end

function LinearAlgebra.rdiv!(
    A::Matrix{T},
    adjB::Adjoint{T,<:LowerTriangular{T,UniformBlockDiagonal{T}}},
) where {T}
    m, n = size(A)
    Bd = adjB.parent.data
    r, s, blk = size(Bd.data)
    n == size(Bd, 1) && r == s || throw(DimensionMismatch())
    for b = 1:blk
        _rdivlowertriblk!(A, Bd.data, b, s)
    end
    A
end

function LinearAlgebra.rdiv!(
    A::BlockedSparse{T,S,P},
    B::Adjoint{T,<:LowerTriangular{T,UniformBlockDiagonal{T}}},
) where {T,S,P}
    Bpd = B.parent.data
    j, k, l = size(Bpd.data)
    cbpt = A.colblkptr
    nzv = A.cscmat.nzval
    P == j == k && length(cbpt) == (l + 1) || throw(DimensionMismatch(""))
    for (j, f) in enumerate(Bpd.facevec)
        rdiv!(reshape(view(nzv, cbpt[j]:(cbpt[j+1]-1)), :, P), adjoint(LowerTriangular(f)))
    end
    A
end


text:
using Surrogates
using Test
using SafeTestsets

@time @safetestset "AbstractGP.jl" begin include("AbstractGP.jl") end
@time @safetestset  "Radials.jl" begin include("Radials.jl") end
@time @safetestset  "Kriging.jl" begin include("Kriging.jl") end
@time @safetestset  "Sampling" begin include("sampling.jl") end
@time @safetestset  "Optimization" begin include("optimization.jl") end
@time @safetestset  "LinearSurrogate" begin include("linearSurrogate.jl") end
@time @safetestset  "Lobachevsky" begin include("lobachevsky.jl") end
@time @safetestset  "RandomForestSurrogate" begin include("random_forest.jl") end
#@time @safetestset  "SVMSurrogate" begin include("SVMSurrogate.jl") end
@time @safetestset  "NeuralSurrogate" begin include("neuralSurrogate.jl") end
@time @safetestset  "InverseDistanceSurrogate" begin include("inverseDistanceSurrogate.jl") end
@time @safetestset  "SecondOrderPolynomialSurrogate" begin include("secondOrderPolynomialSurrogate.jl") end
@time @safetestset  "AD_Compatibility" begin include("AD_compatibility.jl") end
@time @safetestset  "Wendland" begin include("Wendland.jl") end
#@time @safetestset  "MOE" begin include("MOE.jl") end write em algorithm to get rid fo deps
@time @safetestset  "VariableFidelity" begin include("VariableFidelity.jl") end
@time @safetestset  "Earth" begin include("earth.jl") end
@time @safetestset  "Gradient Enhanced Kriging" begin include("GEK.jl") end
@time @safetestset  "Section Samplers" begin include("SectionSampleTests.jl") end


text:
testpath(f) = joinpath(dirname(@__FILE__), f)

for test_file in [
  "simulations.jl"
  ]
  include(testpath(test_file))
end

max_stars_repo_name: zhangliye/SimHPC

max_issues_repo_name: zhangliye/SimHPC

max_forks_repo_name: zhangliye/SimHPC


text:
@doc raw"""

Implements equations for Lagrangian Lotka-Volterra models in 4d of the form
```math
L(q, \dot{q}) = \frac{1}{2} (\log q)^T A \frac{\dot{q}}{q} + q^T B \dot{q} - H (q) ,
```
with hamiltonian
```math
H(q) = a_1 q^1 + a_2 q^2 + a_3 q^3 + a_4 q^4 + b_1 \log q^1 + b_2 \log q^2 + b_3 \log q^3 + b_4 \log q^4 + b_5 \log q^5 .
```

"""
module LotkaVolterra4dLagrangian

    using GeometricBase
    using GeometricEquations
    using ModelingToolkit
    using Parameters
    using RuntimeGeneratedFunctions
    using Symbolics

    RuntimeGeneratedFunctions.init(@__MODULE__)

    # export hamiltonian, ϑ, ϑ₁, ϑ₂, ω

    export lotka_volterra_4d_ode,
           lotka_volterra_4d_iode, lotka_volterra_4d_idae,
           lotka_volterra_4d_lode
        #    lotka_volterra_4d_ldae,
        #    lotka_volterra_4d_dg

    # include("lotka_volterra_4d_plots.jl")

    Δt = 0.1
    nt = 1000

    const q₀ = [ 2.0,  1.0,  1.0,  1.0]
    const a₀ = [ 1.0,  1.0,  1.0,  1.0]
    const b₀ = [-1.0, -2.0, -1.0, -1.0]

    const reference_solution = [1.6390462434739954, 1.3764800055785835, 0.37903204434372284, 1.4399236281802124]

    const A_antisym = 1//2 .* [
                0 -1 +1 -1
               +1  0 -1 +1
               -1 +1  0 -1
               +1 -1 +1  0]

    const A_positive = [
                0  0  1  0
                1  0  0  1
                0  1  0  0
                1  0  1  0]

    const A_upper = [
                0 -1 +1 -1
                0  0 -1 +1
                0  0  0 -1
                0  0  0  0]

    const A_lower = [
                0  0  0  0
               +1  0  0  0
               -1 +1  0  0
               +1 -1 +1  0]

    const A_quasicanonical_antisym = 1//2 .* [
                0 -1 +1 -1
               +1  0 -1  0
               -1 +1  0 -1
               +1  0 +1  0]

    const A_quasicanonical_reduced = 1//2 .* [
                0  0 +1  0
               +2  0 -2  0
               -1  0  0  0
               +2  0 +2  0]

    const B = [ 0  1  1  1
                1  0  1  1
                1  1  0  1
                1  1  1  0]

    H(x, a, b) = sum(a .* x) + sum(b .* log.(x))
    K(x, v, A, B) = sum(log.(x) .* (A * (v ./ x))) + sum(x .* (B * v))
    L(x, v, A, B, a, b) = K(x, v, A, B) - H(x, a, b)

    function substitute_ẋ_with_v!(eqs, ẋ, v)
        for i in eachindex(eqs)
            eqs[i] = substitute(eqs[i], [ẋᵢ=>vᵢ for (ẋᵢ,vᵢ) in zip(ẋ,v)])
        end
    end

    function substitute_variables!(eqs, x, v, X, V)
        for i in eachindex(eqs)
            eqs[i] = substitute(eqs[i], [z=>Z for (z,Z) in zip([x..., v...], [X..., V...])])
        end
    end

    function get_functions(A,B,a,b)
        @parameters t

        @variables x₁(t), x₂(t), x₃(t), x₄(t)
        @variables v₁(t), v₂(t), v₃(t), v₄(t)
        @variables X[1:4]
        @variables V[1:4]
        @variables P[1:4]
        @variables F[1:4]

        x = [x₁, x₂, x₃, x₄]
        v = [v₁, v₂, v₃, v₄]

        Dt = Differential(t)
        Dx = Differential.(x)
        Dv = Differential.(v)
        DX = Differential.(X)
        DV = Differential.(V)

        let L = simplify(L(x, v, A, B, a, b)), K = simplify(K(x, v, A, B)), H = simplify(H(x, a, b))
            EL = [expand_derivatives(Dx[i](L) - Dt(Dv[i](L))) for i in eachindex(Dx,Dv)]
            ∇H = [expand_derivatives(dx(H)) for dx in Dx]
            f  = [expand_derivatives(dx(L)) for dx in Dx]
            f̄  = [expand_derivatives(dx(K)) for dx in Dx]
            g  = [expand_derivatives(Dt(dv(L))) for dv in Dv]
            ϑ  = [expand_derivatives(dv(L)) for dv in Dv]
            ω  = [expand_derivatives(simplify(Dx[i](ϑ[j]) - Dx[j](ϑ[i]))) for i in eachindex(Dx,ϑ), j in eachindex(Dx,ϑ)]
            Σ  = simplify.(inv(ω))
            ẋ  = simplify.(Σ * ∇H)

            for eq in (EL, ∇H, f, f̄, g, ϑ, ω, Σ, ẋ)
                substitute_ẋ_with_v!(eq, Dt.(x), v)
                substitute_variables!(eq, x, v, X, V)
            end

            H = substitute(H, [z=>Z for (z,Z) in zip([x..., v...], [X..., V...])])
            L = substitute(L, [z=>Z for (z,Z) in zip([x..., v...], [X..., V...])])
            ϕ = [P[i] - ϑ[i] for i in eachindex(P,ϑ)]
            ψ = [F[i] - g[i] for i in eachindex(F,g)]

            # component wise code generation
            # code_EL = [build_function(eq, t, X, V)[2] for eq in EL]
            # code_f  = [build_function(eq, t, X, V)[2] for eq in f ]
            # code_f̄  = [build_function(eq, t, X, V)[2] for eq in f̄ ]
            # code_g  = [build_function(eq, t, X, V)[2] for eq in g ]
            # code_∇H = [build_function(eq, t, X)[2]    for eq in ∇H]
            # code_p  = [build_function(eq, t, X)[1]    for eq in ϑ ]
            # code_ϑ  = [build_function(eq, t, X)[2]    for eq in ϑ ]
            # code_ω  = [build_function(eq, t, X)[2]    for eq in ω ]
            # code_P  = [build_function(eq, t, X)[2]    for eq in Σ ]
            # code_ẋ  = [build_function(eq, t, X)[2]    for eq in ẋ ]
            # code_ϕ  = [build_function(eq, t, X, P)[2]        for eq in ϕ ]
            # code_ψ  = [build_function(eq, t, X, V, P, F)[2]  for eq in ψ ]
            # code_H  =  build_function(H,  t, X)

            code_EL = build_function(EL, t, X, V)[2]
            code_f  = build_function(f,  t, X, V)[2]
            code_f̄  = build_function(f̄,  t, X, V)[2]
            code_g  = build_function(g,  t, X, V)[2]
            code_∇H = build_function(∇H, t, X)[2]
            code_p  = build_function(ϑ,  t, X)[1]
            code_ϑ  = build_function(ϑ,  t, X)[2]
            code_ω  = build_function(ω,  t, X)[2]
            code_P  = build_function(Σ,  t, X)[2]
            code_ẋ  = build_function(ẋ,  t, X)[2]
            code_ϕ  = build_function(ϕ,  t, X, P)[2]
            code_ψ  = build_function(ψ,  t, X, V, P, F)[2]
            code_H  = build_function(H,  t, X)
            code_L  = build_function(L,  t, X, V)

            return (
                EL = @RuntimeGeneratedFunction(Symbolics.inject_registered_module_functions(code_EL)),
                ∇H = @RuntimeGeneratedFunction(Symbolics.inject_registered_module_functions(code_∇H)),
                f  = @RuntimeGeneratedFunction(Symbolics.inject_registered_module_functions(code_f)),
                f̄  = @RuntimeGeneratedFunction(Symbolics.inject_registered_module_functions(code_f̄)),
                g  = @RuntimeGeneratedFunction(Symbolics.inject_registered_module_functions(code_g)),
                p  = @RuntimeGeneratedFunction(Symbolics.inject_registered_module_functions(code_p)),
                ϑ  = @RuntimeGeneratedFunction(Symbolics.inject_registered_module_functions(code_ϑ)),
                ω  = @RuntimeGeneratedFunction(Symbolics.inject_registered_module_functions(code_ω)),
                P  = @RuntimeGeneratedFunction(Symbolics.inject_registered_module_functions(code_P)),
                ẋ  = @RuntimeGeneratedFunction(Symbolics.inject_registered_module_functions(code_ẋ)),
                ϕ  = @RuntimeGeneratedFunction(Symbolics.inject_registered_module_functions(code_ϕ)),
                ψ  = @RuntimeGeneratedFunction(Symbolics.inject_registered_module_functions(code_ψ)),
                H  = @RuntimeGeneratedFunction(Symbolics.inject_registered_module_functions(code_H)),
                L  = @RuntimeGeneratedFunction(Symbolics.inject_registered_module_functions(code_L)),
            )
        end
    end

    function lotka_volterra_4d_ode(q₀=q₀, A=A_antisym, B=zeros(Int, 4, 4), a=a₀, b=b₀)
        funcs = get_functions(A,B,a,b)
        ODE((t,x,ẋ) -> funcs[:ẋ](ẋ,t,x), q₀;
            invariants = (h = funcs[:H],))
    end

    function lotka_volterra_4d_iode(q₀=q₀, A=A_antisym, B=zeros(Int, 4, 4), a=a₀, b=b₀)
        funcs = get_functions(A,B,a,b)
        IODE((t,x,v,ϑ) -> funcs[:ϑ](ϑ,t,x),
             (t,x,v,f) -> funcs[:f](f,t,x,v),
             (t,x,v,f̄) -> funcs[:f̄](f̄,t,x,v),
             q₀, funcs[:p](0, q₀);
             v̄ = (t,x,v) -> funcs[:ẋ](v,t,x),
             invariants = (h = (t,x,v) -> funcs[:H](t,x),))
    end

    function lotka_volterra_4d_lode(q₀=q₀, A=A_antisym, B=zeros(Int, 4, 4), a=a₀, b=b₀)
        funcs = get_functions(A,B,a,b)
        LODE((t,x,v,ϑ) -> funcs[:ϑ](ϑ,t,x),
             (t,x,v,f) -> funcs[:f](f,t,x,v),
             (t,x,v,f̄) -> funcs[:f̄](f̄,t,x,v),
             (t,x,v)   -> funcs[:L](t,x,v),
             (t,x,v,ω) -> funcs[:ω](ω,t,x),
             q₀, funcs[:p](0, q₀);
             v̄ = (t,x,v) -> funcs[:ẋ](v,t,x),
             invariants = (h = (t,x,v) -> funcs[:H](t,x),))
    end

    function lotka_volterra_4d_idae(q₀=q₀, A=A_antisym, B=zeros(Int, 4, 4), a=a₀, b=b₀)
        funcs = get_functions(A,B,a,b)
        IDAE((t,x,v,ϑ) -> funcs[:ϑ](ϑ,t,x,v),
             (t,x,v,f) -> funcs[:f](f,t,x,v),
             (t,x,p,v,u) -> u .= v,
             (t,x,p,v,f̄) -> funcs[:f̄](f̄,t,x,v),
             (t,x,p,ϕ) -> funcs[:ϕ](ϕ,t,x,p),
             q₀, funcs[:p](0, q₀), zero(q₀);
             v̄ = (t,x,v) -> funcs[:ẋ](v,t,x),
             invariants = (h = (t,x,v) -> funcs[:H](t,x),))
    end

    # function lotka_volterra_4d_ldae(q₀=q₀, p₀=ϑ(0, q₀), λ₀=zero(q₀), params=p)
    #     LDAE(lotka_volterra_4d_ϑ, lotka_volterra_4d_f_ham,
    #          lotka_volterra_4d_g, lotka_volterra_4d_g̅,
    #          lotka_volterra_4d_ϕ, lotka_volterra_4d_ψ,
    #          q₀, p₀, λ₀; parameters=params, h=hamiltonian, v=lotka_volterra_4d_v)
    # end

    # function lotka_volterra_4d_dg(q₀=q₀, p₀=ϑ(0, q₀), params=p)
    #     IODE(lotka_volterra_4d_ϑ, lotka_volterra_4d_f,
    #          lotka_volterra_4d_g, q₀, p₀;
    #          parameters=params, h=hamiltonian, v=lotka_volterra_4d_v)
    # end

end


text:
#### Soil model

export SoilModel, SoilParamFunctions

"""
    AbstractSoilParameterFunctions{FT <: AbstractFloat}
"""
abstract type AbstractSoilParameterFunctions{FT <: AbstractFloat} end

"""
    SoilParamFunctions{FT} <: AbstractSoilParameterFunctions{FT}

Necessary parameters for the soil model. These will eventually be prescribed
functions of space (and time).
# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct SoilParamFunctions{FT} <: AbstractSoilParameterFunctions{FT}
    "Aggregate porosity of the soil"
    porosity::FT = FT(NaN)
    "Hydraulic conductivity at saturation. Units of m s-1."
    Ksat::FT = FT(NaN)
    "Specific storage of the soil. Units of m s-1."
    S_s::FT = FT(NaN)
    "Residual Water Fraction - default is zero; unitless."
    θ_r::FT = FT(0.0)
    "Volume fraction of gravels, relative to soil solids only; unitless."
    ν_ss_gravel::FT = FT(NaN)
    "Volume fraction of SOM, relative to soil solids only; unitless."
    ν_ss_om::FT = FT(NaN)
    "Volume fraction of quartz, relative to soil solids only; unitless."
    ν_ss_quartz::FT = FT(NaN)
    "Bulk volumetric heat capacity of dry soil. Units of J m-3 K-1."
    ρc_ds::FT = FT(NaN)
    "Particle density for soil solids. Units of kg m-3"
    ρp::FT = FT(NaN)
    "Thermal conductivity of the soil solids. Units of W m-1 K-1."
    κ_solid::FT = FT(NaN)
    "Saturated thermal conductivity for unfrozen soil. Units of W m-1 K-1."
    κ_sat_unfrozen::FT = FT(NaN)
    "Saturated thermal conductivity for frozen soil. Units of W m-1 K-1."
    κ_sat_frozen::FT = FT(NaN)
    "Adjustable scale parameter for determining Kersten number in the Balland and Arp formulation; unitless."
    a::FT = FT(0.24)
    "Adjustable scale parameter for determining Kersten number in the Balland and Arp formulation; unitless."
    b::FT = FT(18.1)
    "Parameter used in the Balland and Arp formulation for κ_dry; unitless"
    κ_dry_parameter::FT = FT(0.053)

end

"""
    SoilModel{PF, W, H} <: BalanceLaw

A BalanceLaw for soil modeling.
Users may over-ride prescribed default values for each field.

# Usage

    SoilModel(
        param_functions,
        water,
        heat,
    )

# Fields
$(DocStringExtensions.FIELDS)
"""
struct SoilModel{PF, W, H} <: BalanceLaw
    "Soil Parameter Functions"
    param_functions::PF
    "Water model"
    water::W
    "Heat model"
    heat::H
end

function vars_state(soil::SoilModel, st::Prognostic, FT)
    @vars begin
        water::vars_state(soil.water, st, FT)
        heat::vars_state(soil.heat, st, FT)
    end
end

function vars_state(soil::SoilModel, st::Auxiliary, FT)
    @vars begin
        water::vars_state(soil.water, st, FT)
        heat::vars_state(soil.heat, st, FT)
    end
end

function vars_state(soil::SoilModel, st::Gradient, FT)
    @vars begin
        water::vars_state(soil.water, st, FT)
        heat::vars_state(soil.heat, st, FT)
    end
end

function vars_state(soil::SoilModel, st::GradientFlux, FT)
    @vars begin
        water::vars_state(soil.water, st, FT)
        heat::vars_state(soil.heat, st, FT)
    end
end

function compute_gradient_argument!(
    land::LandModel,
    soil::SoilModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    compute_gradient_argument!(land, soil, soil.heat, transform, state, aux, t)
    compute_gradient_argument!(land, soil, soil.water, transform, state, aux, t)
end

function compute_gradient_flux!(
    land::LandModel,
    soil::SoilModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)

    compute_gradient_flux!(
        land,
        soil,
        soil.water,
        diffusive,
        ∇transform,
        state,
        aux,
        t,
    )
    compute_gradient_flux!(
        land,
        soil,
        soil.heat,
        diffusive,
        ∇transform,
        state,
        aux,
        t,
    )

end

function land_nodal_update_auxiliary_state!(
    land::LandModel,
    soil::SoilModel,
    state::Vars,
    aux::Vars,
    t::Real,
)
    land_nodal_update_auxiliary_state!(land, soil, soil.water, state, aux, t)
    land_nodal_update_auxiliary_state!(land, soil, soil.heat, state, aux, t)
end

function land_init_aux!(
    land::LandModel,
    soil::SoilModel,
    aux::Vars,
    geom::LocalGeometry,
)
    soil_init_aux!(land, soil, soil.water, aux, geom)
    soil_init_aux!(land, soil, soil.heat, aux, geom)
end

"""
    abstract type AbstractSoilComponentModel <: BalanceLaw
"""
abstract type AbstractSoilComponentModel <: BalanceLaw end

## When PrescribedModels are chosen, all balance law functions and boundary
## condition functions should do nothing. Since these models are of super
## type AbstractSoilComponentModel, use AbstractSoilComponentModel in
## argument type. The more specific Water and Heat models will have
## different methods (see soil_water.jl and soil_heat.jl).

vars_state(m::AbstractSoilComponentModel, st::AbstractStateType, FT) = @vars()

function soil_init_aux!(
    land::LandModel,
    soil::SoilModel,
    m::AbstractSoilComponentModel,
    aux::Vars,
    geom::LocalGeometry,
) end

function land_nodal_update_auxiliary_state!(
    land::LandModel,
    soil::SoilModel,
    m::AbstractSoilComponentModel,
    state::Vars,
    aux::Vars,
    t::Real,
)

end

function compute_gradient_argument!(
    land::LandModel,
    soil::SoilModel,
    m::AbstractSoilComponentModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)

end

function compute_gradient_flux!(
    land::LandModel,
    soil::SoilModel,
    m::AbstractSoilComponentModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)

end


text:
module CompilerUtils

using Core: OpaqueClosure
using Core.Compiler:
    CodeInfo,
    CodeInstance,
    NativeInterpreter,
    OptimizationParams,
    OptimizationState,
    compact!,
    convert_to_ircode,
    copy,
    copy_exprargs,
    slot2reg,
    specialize_method

function methodinstance_from_expr_new_opaque_closure(inst::Expr)
    if !Meta.isexpr(inst, :new_opaque_closure)
        error("not an `Expr(new_opaque_closure, ...)`")
    end

    # Ref: jl_new_opaque_closure
    argt, _isva, _rt_lb, rt_ub, m, = inst.args
    argt::Type{<:Tuple}
    rt_ub::Type
    m::Method
    sig = Tuple{OpaqueClosure{argt,rt_ub},argt.parameters...}
    mi = specialize_method(m, sig, Core.svec())
    if isdefined(mi, :cache)
        return mi
    end
    if !Base.uncompressed_ir(m).inferred
        return mi
    end

    # If the `CodeInfo` is already inferred (likely from Tapir) and not cached,
    # insert a `CodeInstance` to avoid confusing `code_llvm` etc.
    inferred_result = m.source
    codeinst = CodeInstance(
        mi,
        Any,              # rettype
        nothing,          # inferred_const
        inferred_result,
        Int32(0x00),      # const_flags
        m.primary_world,  # min_world
        typemax(UInt),    # max_world
    )
    ccall(:jl_mi_cache_insert, Cvoid, (Any, Any), mi, codeinst)

    return mi
end

function ircode_from_codeinfo(ci::CodeInfo)
    # Making a copy here, as, e.g., `convert_to_ircode` mutates `ci`:
    ci = copy(ci)

    for (i, t) in pairs(ci.ssavaluetypes)
        if is_intermediate_value_type(t)
            ci.ssavaluetypes[i] = Any
        end
    end

    linfo = ci.parent  # MethodInstance
    interp = NativeInterpreter()
    params = OptimizationParams(interp)
    opt = OptimizationState(linfo, ci, params, interp)
    nargs = Int(opt.nargs) - 1
    preserve_coverage = false
    code = copy_exprargs(ci.code)
    ir = convert_to_ircode(ci, code, preserve_coverage, nargs, opt)
    ir = slot2reg(ir, ci, nargs, opt)
    ir = compact!(ir)
    return ir
end

end  # module CompilerUtils


text:
## Determine the value of R_HOME from the environment or through evaluating `R RHOME`
const rh = haskey(ENV,"R_HOME") ? ENV["R_HOME"] : (ENV["R_HOME"] = chomp(readall(`R RHOME`)))

## Create and test the libR path
const libR = joinpath(rh,"lib",string("libR.",Base.Sys.shlib_ext))
dlopen_e(libR) == C_NULL && error("Unable to load $libR\n\nPlease re-run Pkg.build(package), and restart Julia.")

## Check for the values of other environment variables expected by R
envset = true
const envkeys = ["R_HOME","R_DOC_DIR","R_INCLUDE_DIR","R_SHARE_DIR","LD_LIBRARY_PATH"]
for kk in envkeys
    if !haskey(ENV,kk)
        envset = false
        break
    end
end

if !envset                              # extract the values through Rscript
    vals = split(readall(`Rscript -e 'for (nm in c("R_HOME","R_DOC_DIR","R_INCLUDE_DIR","R_SHARE_DIR","LD_LIBRARY_PATH")) print(Sys.getenv(nm),quote=FALSE)'`),'\n')
    for i in 1:length(envkeys)
        ENV[envkeys[i]] = vals[i][5:end]  # first part of each value is "[1] "
    end
end

## Write the deps.jl file
open("./deps.jl","w") do io
    println(io,"# This is an auto-generated file; do not edit\n")
    println(io,"const libR=\"",libR,"\"")
    println(io,"ENV[\"R_HOME\"]=\"",ENV["R_HOME"],"\"")
    println(io,"ENV[\"R_DOC_DIR\"]=\"",ENV["R_DOC_DIR"],"\"")
    println(io,"ENV[\"R_INCLUDE_DIR\"]=\"",ENV["R_INCLUDE_DIR"],"\"")
    println(io,"ENV[\"R_SHARE_DIR\"]=\"",ENV["R_SHARE_DIR"],"\"")
    println(io,"ENV[\"LD_LIBRARY_PATH\"]=\"",ENV["LD_LIBRARY_PATH"],"\"")
end


text:
using Test, Random, Distributed, Statistics
Random.seed!(0)
using Hyperopt, Plots
using Optim

@testset "BOHB" begin
    f(a;c=10) = sum(@. 100 + (a-3)^2 + (c-100)^2)
    @test_nowarn Hyperband(50)
    let bohb = @hyperopt for i=18, sampler=Hyperband(R=50, η=3, inner=BOHB(dims=[Hyperopt.Continuous(), Hyperopt.Continuous()])), a = LinRange(1,5,800), c = exp10.(LinRange(-1,3,1800))
            # println(i, "\t", a, "\t", b, "\t", c)
            if !(state === nothing)
                a,c = state
            end
            res = Optim.optimize(x->f(x[1],c=x[2]), [a,c], NelderMead(), Optim.Options(f_calls_limit=i))
            Optim.minimum(res), Optim.minimizer(res)
        end
        @test length(bohb.history) == 69
        @test length(bohb.results) == 69
        @test minimum(bohb) < 300
        @hyperopt for i=1, ho=bohb, sampler=Hyperband(R=50, η=3, inner=BOHB(dims=[Hyperopt.Continuous(), Hyperopt.Continuous()])), a = LinRange(1,5,800), c = exp10.(LinRange(-1,3,1800))
            # println(i, "\t", a, "\t", b, "\t", c)
            if !(state === nothing)
                a,c = state
            end
            res = Optim.optimize(x->f(x[1],c=x[2]), [a,c], NelderMead(), Optim.Options(f_calls_limit=i))
            Optim.minimum(res), Optim.minimizer(res)
        end
        @test length(bohb.history) == 138
        @test length(bohb.results) == 138
    end

    # extra robust option
    @macroexpand @hyperopt for i=18, sampler=Hyperband(R=50, η=3, inner=BOHB(dims=[Hyperopt.UnorderedCategorical(5), Hyperopt.Continuous(), Hyperopt.Continuous()])),
        algorithm = [SimulatedAnnealing(), ParticleSwarm(), NelderMead(), BFGS(), NewtonTrustRegion()],
        a = LinRange(1,5,1800),
        c = exp10.(LinRange(-1,3,1800))
        if !(state === nothing)
            a,c,algorithm = state
        end
        println(i, " algorithm: ", typeof(algorithm))
        res = Optim.optimize(x->f(x[1],c=x[2]), [a,c], algorithm, Optim.Options(time_limit=2i+2, show_trace=false, show_every=5))
        Optim.minimum(res), (Optim.minimizer(res)..., algorithm)
    end

    bohb = @hyperopt for i=18, sampler=Hyperband(R=50, η=3, inner=RandomSampler()),
        algorithm = [SimulatedAnnealing(), ParticleSwarm(), NelderMead(), BFGS(), NewtonTrustRegion()],
        a = LinRange(1,5,1800),
        c = exp10.(LinRange(-1,3,1800))
        if !(state === nothing)
            a,c,algorithm = state
        end
        println(i, " algorithm: ", typeof(algorithm))
        res = Optim.optimize(x->f(x[1],c=x[2]), [a,c], algorithm, Optim.Options(time_limit=2i+2, show_trace=false, show_every=5))
        Optim.minimum(res), (Optim.minimizer(res)..., algorithm)
    end

    # ho = @hyperopt for i=18, sampler=Hyperband(R=50, η=3, inner=BOHB(dims=[Hyperopt.UnorderedCategoricalDim(5), Hyperopt.ContinuousDim(), Hyperopt.ContinuousDim()])),
    ho = @hyperopt for i=18, sampler=Hyperband(R=50, η=3, inner=BOHB(dims=[Hyperopt.UnorderedCategorical(5), Hyperopt.Continuous(), Hyperopt.Continuous()])),
        algorithm = [SimulatedAnnealing(), ParticleSwarm(), NelderMead(), BFGS(), NewtonTrustRegion()],
        a = LinRange(1,5,1800),
        c = exp10.(LinRange(-1,3,1800))
        if !(state === nothing)
            a,c,algorithm = state
        end
        println(i, " algorithm: ", typeof(algorithm))
        res = Optim.optimize(x->f(x[1],c=x[2]), [a,c], algorithm, Optim.Options(time_limit=2i+2, show_trace=false, show_every=5))
        Optim.minimum(res), (Optim.minimizer(res)..., algorithm)
    end

end


text:
module CoolTensors

export Tensor, @t_str, raise, lower, TCartesianIndex
export TScalar, TVector, TCovector, TLinearMap, TBilinear

struct Tensor{T,N,ipos,A<:AbstractArray{T,N}} <: AbstractArray{T,N}
    parent::A
end

Base.size(t::Tensor) = size(t.parent)
function Base.getindex(t::Tensor, inds::Int...)
    @boundscheck checkbounds(t.parent, inds...)
    @inbounds t.parent[inds...]
end

struct IndexPos{N}
    # TODO: Should this be UInt64 instead?
    x::UInt
end

index_pos(::Tensor{<:Any,<:Any,ipos}) where {ipos} = ipos
function (ipos::IndexPos{N})(a::AbstractArray{T,N}) where {T,N}
    Tensor{T,N,ipos,typeof(a)}(a)
end
(::IndexPos{0})(x::Number) = x

Base.getindex(ipos::IndexPos, i) = (ipos.x >> (i-1)) % Bool
Base.lastindex(::IndexPos{N}) where {N} = N
function deleteat(ipos::IndexPos{N}, i::Int) where {N}
    mask = (UInt(1) << i) - UInt(1)
    IndexPos{N-1}((ipos.x & (mask >> 1)) | ((ipos.x & ~mask) >> 1))
end
function Base.show(io::IO, ipos::IndexPos{N}) where {N}
    print(io, "t\"")
    for i in 1:N
        print(io, ipos[i] ? ''' : ',')
    end
    print(io, '\"')
    nothing
end

macro t_str(s)
    x = mapreduce(|, enumerate(s), init=UInt(0)) do (i, c)
        if c === '''
            UInt(1) << (i-1)
        elseif c === ','
            UInt(0)
        else
            throw(ArgumentError("Index positions can only be ''' and ','"))
        end
    end
    :(IndexPos{$(length(s))}($x))
end

for (typed_f, f) in [
    (:getindex, :vect), (:typed_hcat, :hcat), (:typed_vcat, :vcat), (:typed_hvcat, :hvcat),
]
    @eval Base.$typed_f(ipos::IndexPos, args...) = ipos(Base.$f(args...))
end
Base.getindex(::IndexPos{0}, x) = t""(fill(x))

const TScalar = Tensor{T,0,t"",<:AbstractArray{T,0}} where {T}
const TVector = Tensor{T,1,t"'",<:AbstractVector{T}} where {T}
const TCovector = Tensor{T,1,t",",<:AbstractVector{T}} where {T}
const TLinearMap = Tensor{T,2,t"',",<:AbstractMatrix{T}} where {T}
const TBilinear = Tensor{T,2,t",,",<:AbstractMatrix{T}} where {T}

function (t::Tensor{<:Any,N})(v::Union{TVector,TCovector}; dims=N) where {N}
    N == 0 && throw(ArgumentError("Cannot contract Scalar with $(typeof(v))"))
    index_pos(t)[dims] !== index_pos(v)[1] ||
        throw(ArgumentError("Tensor indices don't match"))

    # TODO: use an external package for tensor reductions
    parent = mapreduce(
        (x, y) -> x .* y, (x, y) -> x .+ y, eachslice(t.parent, dims=dims), v;
        init=zeros(promote_type(eltype(t), eltype(v)), size(t)[setdiff(1:end, dims)]),
    )
    deleteat(index_pos(t), dims)(parent)
end
(t::Tensor)(::Colon; dims=0) = t
function (t::Tensor{<:Any,N})(args...) where {N}
    foldl((t, (i,v)) -> t(v, dims=N-i+1), enumerate(args), init=t)
end

function raise(t::Tensor{<:Any,N,ipos}, i::Int) where {N,ipos}
    IndexPos{N}(ipos.x | (UInt(1) << (i-1)))(t.parent)
end
raise(t::Tensor, inds::Vararg{Int}) = foldl((t, i) -> raise(t, i), inds, init=t)

function lower(t::Tensor{<:Any,N,ipos}, i::Int) where {N,ipos}
    IndexPos{N}(ipos.x & ~(UInt(1) << (i-1)))(t.parent)
end
lower(t::Tensor, inds::Vararg{Int}) = foldl((t, i) -> lower(t, i), inds, init=t)

struct TCartesianIndex{N,ipos}
    I::CartesianIndex{N}
end

function TCartesianIndex{N,ipos}(t::NTuple{N,Int}) where {N,ipos}
    TCartesianIndex{N,ipos}(CartesianIndex(t))
end
TCartesianIndex{N,ipos}(t::Vararg{Int,N}) where {N,ipos} = TCartesianIndex{N,ipos}(t)

function Base.getindex(t::Tensor{<:Any,N,ipos}, cI::TCartesianIndex{N,ipos}) where {N,ipos}
    t.parent[cI.I]
end

Base.typed_vcat(t::Tensor, i...) = Base.typed_hvcat(t, ntuple(one, length(i)), i...)
Base.typed_hcat(t::Tensor, i...) = Base.typed_hvcat(t, (length(i),), i...)

function index_pos_from_alternations(alt::Tuple, current_ipos=true)
    i = 0
    x = UInt(0)
    while !isempty(alt)
        if first(alt) > 0
            x |= UInt(current_ipos) << i
            i += 1
            alt = (first(alt) - 1, Base.tail(alt)...)
        else
            current_ipos = !current_ipos
            alt = Base.tail(alt)
        end
    end
    IndexPos{i}(x)
end

function Base.typed_hvcat(t::Tensor{<:Any,N,ipos}, alt::Tuple, i...) where {N,ipos}
    i, _ipos = if !isempty(i) && first(i) === (\)
        Base.tail(i), index_pos_from_alternations(Base.tail(alt), false)
    else
        i, index_pos_from_alternations(alt)
    end
    _ipos === ipos ||
        throw(ArgumentError("Index positions $_ipos don't match indices of Tensor $ipos"))
    I = TCartesianIndex{N,ipos}(reduce((x, y) -> (x..., y...), i))
    t[I]
end

end


text:
abstract type AbstractCoordinatePoint{T, S} <: StaticArrays.FieldVector{3, T} end
abstract type AbstractCoordinateVector{T, S} <: StaticArrays.FieldVector{3, T} end

include("Points.jl")
include("Vectors.jl")


text:
using LinearAlgebra

@testset "lsa" begin
    crps = Corpus([StringDocument("this is a string document"), TokenDocument("this is a token document")])
    F1 = lsa(crps)

    update_lexicon!(crps)
    m = DocumentTermMatrix(crps)

    # lsa on crps
    @test typeof(F1) <: LinearAlgebra.SVD
    @test F1.U * LinearAlgebra.Diagonal(F1.S) * F1.Vt == Matrix(tf_idf(m))

    # lsa on dtm
    F2 = lsa(m)
    @test typeof(F2) <: LinearAlgebra.SVD
    @test F2.U * LinearAlgebra.Diagonal(F2.S) * F2.Vt == Matrix(tf_idf(m))

    @test F1 == F2
end


text:
# Linear Program Contractor
function poorLP_contractor(f,g,X::Array{Interval},UBD,nx::Int64)
    # gets prior McCormick relaxation setting
    mu_temp = EAGOSmoothMcCormickGrad.MC_param.mu
    set_diff_relax!(0)

    # sets variable bounds
    l = [X[i].lo for i=1:nx]
    u = [X[i].lo for i=1:nx]
    l_out = copy(l)
    u_out = copy(u)

    # computes relaxation of g and f
    x_mc = [SMCg{nx,Float64}(mid.(X[i]),mid.(X[i]),seed(nx,i),seed(nx,i),X[i],false,X,mid.(X))]
    f_mc = f(x)
    g_mc = g(x)

    ### constructs linear problem
    A = zeros(length(g_mc)+1,nx)
    b = zeros(length(g_mc)+1)
    A[1,:] = [f.cv_grad[i] for i=1:nx]
    b[1] = UBD - f.cv + sum([f.cv_grad[i]*mid.(X[i]) for i=1:nx])
    A[2:length(g_mc)+1,:] = [g_mc[i].cv_grad[j] for i=1:length(g_mc),j=1:nx]
    b[2:length(g_mc)+1] = [sum([g_mc[i].cv_grad[j]*mid.(X[i]) for j=1:nx])-g_mc[i].cv for i=1:length(g_mc)]

    for hind=1:nx
        for i=1:(length(g_mc)+1)
            if (A[i,hind] > 0.0)
                upper_bnd = (1.0/A[i,hind])*(b[i]-sum([ hind == j ? 0.0 : min(A[i,j]*u_out[j],A[i,j]*l_out[j]) for j=1:nx]))
                if (u_out[hind] > upper_bnd)
                    u_out[hind] = upper_bnd
                end
            elseif (A[i,hind] < 0.0)
                lower_bnd = (1.0/A[i,hind])*(b[i]-sum([ hind == j ? 0.0 : min(A[i,j]*u_out[j],A[i,j]*l_out[j]) for j=1:nx]))
                if (l_out[hind] < lower_bnd)
                    u_out[hind] = upper_bnd
                end
            end
        end
        if (lout[hind]>uout[hind])
            return X, false
        end
    end

    # reset McCormick relaxation to original value
    set_diff_relax!(mu_temp)

    temp_arr = [Interval(l_out[i],u_out[i]) for i=1:nx]
    return temp_arr,true
end


text:
export rotate

function rotate(x::Array{Float64,2}; angle=0.0, point=[0.0, 0.0, 0.0], vector_=[1.0, 0.0, 0.0])
    t = angle/180*pi
    x = x .- vec(point)
    vector_ = vec(vector_)
    unit_vector = vector_/(sqrt(sum(vector_.^2)))
    l, m, n = unit_vector
    c, s = cos(t), sin(t)
    rot = [l*l*(1-c)+c     m*l*(1-c)-n*s   n*l*(1-c)+m*s;
           l*m*(1-c)+n*s   m*m*(1-c)+c     n*m*(1-c)-l*s;
           l*n*(1-c)-m*s   m*n*(1-c)+l*s   n*n*(1-c)+c;]
    x = rot * x
    x = x .+ vec(point)
    return x
end

function rotate(obj::T; angle=0.0, point=[0.0, 0.0, 0.0], vector_=[1.0, 0.0, 0.0]) where T <: Union{Shape, PostOpObj}
    # x, v, y, vol, type
    function func(out)
        x = out[:x]
        x = rotate(x; angle=angle, point=point, vector_=vector_)
        return repack!(out, [:x], [x])
    end
    if isa(obj, Shape)
        return PostOpObj(obj, func)
    elseif isa(obj, PostOpObj)
        push!(obj.operations, func)
        return obj
    else
        error("Not allowed.")
    end
end

max_stars_repo_name: sajidmannan/PDMesh

max_issues_repo_name: sajidmannan/PDMesh

max_forks_repo_name: sajidmannan/PDMesh


text:
# --- Utilities ---
# That could be contributed to Julia.

@inline ⩔(f, g) = x -> f(x) | g(x)

"""
    readuntil(predicate, io::IO)

Read characters until matching a predicate.

Based on [this](https://github.com/JuliaLang/julia/issues/21355#issue-221121166).

```jldoctest
julia> readuntil(isspace, IOBuffer("pineapple: 🍍..."))
"pineapple:"
```
"""
@inline function Base.readuntil(predicate, io::IO)
    buf = IOBuffer()
    while !eof(io)
        c = read(io, Char)
        if predicate(c)
            skip(io, -ncodeunits(c))
            break
        end
        write(buf, c)
    end
    return String(take!(buf))
end

"""
    startswith(io::IO, prefix::Union{AbstractString,Base.Chars})

Check if an `IO` object starts with a prefix.

Based on [this](https://github.com/JuliaLang/julia/issues/40616#issue-867861851).

```jldoctest
julia> io = IOBuffer("pineapple: 🍍...");

julia> startswith(io, "pine")
true

julia> startswith(io, "apple")
false
```
"""
@inline function Base.startswith(io::IO, prefix::Union{AbstractString,Base.Chars})
    pos = position(io)
    s = _getminprefix(io, prefix)
    seek(io, pos)
    return startswith(s, prefix)
end
@inline _getminprefix(io::IO, prefix::Union{AbstractString,AbstractChar}) = String(read(io, length(prefix)))
@inline _getminprefix(io::IO, chars::Union{Tuple{Vararg{<:AbstractChar}},AbstractVector{<:AbstractChar},Set{<:AbstractChar}}) = _getminprefix(io, first(chars))

"""
    skipstartswith(io::IO, prefix::Union{AbstractString,Base.Chars})

Check if an `IO` object starts with a prefix and skip it.

```jldoctest
julia> io = IOBuffer("pineapple: 🍍...");

julia> HTM.skipstartswith(io, "pine")
true

julia> read(io, String)
"apple: 🍍..."

julia> io = IOBuffer("pineapple: 🍍...");

julia> HTM.skipstartswith(io, "apple")
false

julia> read(io, String)
"pineapple: 🍍..."
```
"""
@inline skipstartswith(io::IO, prefix::Union{AbstractString,AbstractChar}) = startswith(io, prefix) ? (skip(io, length(prefix)); true) : false
# @inline skipstartswith(io::IO, chars::Union{Tuple{Vararg{<:AbstractChar}},AbstractVector{<:AbstractChar},Set{<:AbstractChar}}) = startswith(io, chars) ? (skip(io, length(first(chars))); true) : false


text:
using Documenter
using WebIO

# We have to ensure that these modules are loaded because some functions are
# defined behind @require guards.
using IJulia, Mux, Blink

makedocs(
    sitename="WebIO",
    format=Documenter.HTML(),
    modules=[WebIO],
    pages=[
        "index.md",
        "gettingstarted.md",
        "API Reference" => [
            "api/about.md",
            "api/node.md",
            "api/scope.md",
            "api/render.md",
            "api/asset.md",
            "api/jsstring.md",
            "api/observable.md",
        ],
        "Providers" => [
            "providers/ijulia.md",
            "providers/blink.md",
            "providers/mux.md",
        ],
        "extending.md",
    ],
)

deploydocs(
    repo="github.com/JuliaGizmos/WebIO.jl.git",
)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule CPUInfo_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("CPUInfo")
JLLWrappers.@generate_main_file("CPUInfo", UUID("48a17e46-e8cd-5edf-9ac6-c264fe76544e"))
end  # module CPUInfo_jll


text:
import AbstractFFTs: plan_fft, plan_fft!, plan_bfft, plan_bfft!,
    plan_rfft, plan_brfft, plan_inv, normalization, fft, bfft, ifft, rfft,
    Plan, ScaledPlan

using LinearAlgebra

const ROCFFT_FORWARD = true
const ROCFFT_INVERSE = false

abstract type ROCFFTPlan{T,K,inplace} <: Plan{T} end

Base.unsafe_convert(::Type{rocfft_plan}, p::ROCFFTPlan) = p.handle

function unsafe_free!(plan::ROCFFTPlan)
    rocfft_plan_destroy(plan.handle)
    rocfft_execution_info_destroy(plan.execution_info)
end

# TODO: Real to Complex full not possible atm
# K is flag for forward/inverse
mutable struct cROCFFTPlan{T,K,inplace,N} <: ROCFFTPlan{T,K,inplace}
    handle::rocfft_plan
    workarea::ROCArray{Int8,1}
    execution_info::rocfft_execution_info
    sz::NTuple{N,Int} # Julia size of input array
    osz::NTuple{N,Int} # Julia size of output array
    xtype::rocfft_transform_type
    region
    pinv::ScaledPlan # required by AbstractFFTs API

    function cROCFFTPlan{T,K,inplace,N}(handle::rocfft_plan, workarea::ROCArray{Int8,1}, X::ROCArray{T,N}, sizey::Tuple,
                                      xtype::rocfft_transform_type, region) where {T,inplace,N,K}
        info_ref = Ref{rocfft_execution_info}()
        rocfft_execution_info_create(info_ref)
        info = info_ref[]
        rocfft_execution_info_set_work_buffer(info, workarea, length(workarea))
        p = new(handle, workarea, info, size(X), sizey, xtype, region)
        finalizer(unsafe_free!, p)
        p
    end
end

mutable struct rROCFFTPlan{T,K,inplace,N} <: ROCFFTPlan{T,K,inplace}
    handle::rocfft_plan
    workarea::ROCArray{Int8,1}
    execution_info::rocfft_execution_info
    sz::NTuple{N,Int} # Julia size of input array
    osz::NTuple{N,Int} # Julia size of output array
    xtype::rocfft_transform_type
    region
    pinv::ScaledPlan # required by AbstractFFTs API

    function rROCFFTPlan{T,K,inplace,N}(handle::rocfft_plan, workarea::ROCArray{Int8,1}, X::ROCArray{T,N}, sizey::Tuple,
                                      xtype::rocfft_transform_type, region) where {T,inplace,N,K}
        info_ref = Ref{rocfft_execution_info}()
        rocfft_execution_info_create(info_ref)
        info = info_ref[]
        rocfft_execution_info_set_work_buffer(info, workarea, length(workarea))
        p = new(handle, workarea, info, size(X), sizey, xtype, region)
        finalizer(unsafe_free!, p)
        p
    end
end

const xtypenames = ("complex forward", "complex inverse", "real forward", "real inverse")

function showfftdims(io, sz, T)
    if isempty(sz)
        print(io,"0-dimensional")
    elseif length(sz) == 1
        print(io, sz[1], "-element")
    else
        print(io, join(sz, "×"))
    end
    print(io, " ROCArray of ", T)
end

function Base.show(io::IO, p::ROCFFTPlan{T,K,inplace}) where {T,K,inplace}
    print(io, inplace ? "rocFFT in-place " : "rocFFT ",
          xtypenames[Int(p.xtype)+1],
          " plan for ")
    showfftdims(io, p.sz, T)
end

Base.size(p::ROCFFTPlan) = p.sz

function create_plan(xtype::rocfft_transform_type, xdims, T, inplace, region)
    precision = (real(T) == Float64) ? rocfft_precision_double : rocfft_precision_single
    sz = [xdims[i] for i in region]
    csz = copy(sz)
    cxdims = [xdims...]
    csz[1] = div(sz[1],2) + 1
    cxdims[region[1]] = csz[1]
    nrank = length(region)
    batch = prod(xdims) ÷ prod(sz)

    handle_ref = Ref{rocfft_plan}()
    worksize_ref = Ref{Csize_t}()
    placement = inplace ? rocfft_placement_inplace : rocfft_placement_notinplace
    rsz = (length(sz) > 1) ? rsz = reverse(sz) : sz
    if batch == 1
        rocfft_plan_create(handle_ref, placement, xtype, precision, nrank, sz, 1, C_NULL)
    else
        cdims = collect(xdims)
        cdims[region[1]] = div(cdims[region[1]],2)+1
        plan_desc_ref = Ref{rocfft_plan_description}()
        rocfft_plan_description_create(plan_desc_ref)
        description = plan_desc_ref[]
        if xtype == rocfft_transform_type_real_forward
            in_array_type = rocfft_array_type_real
            # TODO: output to full array
            out_array_type = rocfft_array_type_hermitian_interleaved
        elseif xtype == rocfft_transform_type_real_inverse
            # TODO: also for complex_interleaved
            in_array_type = rocfft_array_type_hermitian_interleaved
            out_array_type = rocfft_array_type_real
        else
            in_array_type = rocfft_array_type_complex_interleaved
            out_array_type = rocfft_array_type_complex_interleaved
        end
        if ((region...,) == ((1:nrank)...,))
            # handle simple case ... simply! (for robustness)
            rocfft_plan_description_set_data_layout(description, in_array_type, out_array_type,
                                                    C_NULL, C_NULL,
                                                    nrank, C_NULL, 0,
                                                    nrank, C_NULL, 0)
            rocfft_plan_create(handle_ref, placement, xtype, precision, nrank, sz, batch, description)
        else
            strides = [prod(xdims[1:region[k]-1]) for k=1:nrank]
            real_strides = [prod(cdims[1:region[k]-1]) for k=1:nrank]
            if nrank==1 || all(diff(collect(region)) .== 1)
                # _stride: successive elements in dimension of region
                # _dist: distance between first elements of batches
                if region[1] == 1
                    idist = prod(sz)
                    cdist = prod(csz)
                else
                    if region[end] != length(xdims)
                        throw(ArgumentError("batching dims must be sequential"))
                    end
                    idist = 1
                    cdist = 1
                end
                ostrides = copy(strides)
                istrides = copy(strides)
                if xtype == rocfft_transform_type_real_forward
                    odist = cdist
                    ostrides .= real_strides
                else
                    odist = idist
                end
                if xtype == rocfft_transform_type_real_inverse
                    idist = cdist
                    istrides .= real_strides
                end
            else
                if any(diff(collect(region)) .< 1)
                    throw(ArgumentError("region must be an increasing sequence"))
                end

                if region[1] == 1
                    ii=1
                    while (ii < nrank) && (region[ii] == region[ii+1]-1)
                        ii += 1
                    end
                    idist = prod(xdims[1:ii])
                    cdist = prod(cdims[1:ii])
                    ngaps = 0
                else
                    istride = prod(xdims[1:region[1]-1])
                    idist = 1
                    cdist = 1
                    ngaps = 1
                end
                nem = ones(Int,nrank)
                cem = ones(Int,nrank)
                id = 1
                for ii=1:nrank-1
                    if region[ii+1] > region[ii]+1
                        ngaps += 1
                    end
                    while id < region[ii+1]
                        nem[ii] *= xdims[id]
                        cem[ii] *= cdims[id]
                        id += 1
                    end
                    @assert nem[ii] >= sz[ii]
                end
                if region[end] < length(xdims)
                    ngaps += 1
                end
                # CUFFT represents batches by a single stride (_dist)
                # ROCFFT can have multiple strides, but non sequential batch dims are also not working
                # so we must verify that region is consistent with this:
                if ngaps > 1
                    throw(ArgumentError("batch regions must be sequential"))
                end

                ostrides = copy(strides)
                istrides = copy(strides)
                if xtype == rocfft_transform_type_real_forward
                    odist = cdist
                    ostrides .= real_strides
                else
                    odist = idist
                end
                if xtype == rocfft_transform_type_real_inverse
                    idist = cdist
                    istrides .= real_strides
                end
            end
            rocfft_plan_description_set_data_layout(description, in_array_type, out_array_type,
                                                    C_NULL, C_NULL,
                                                    length(istrides), istrides, idist,
                                                    length(ostrides), ostrides, odist)
            rocfft_plan_create(handle_ref, placement, xtype, precision, nrank, sz, batch, description)
            rocfft_plan_description_destroy(description)
        end
    end
    rocfft_plan_get_work_buffer_size(handle_ref[], worksize_ref)
    # TODO allow empty array in ../array.jl
    workarea = worksize_ref[]>0 ? ROCArray{Int8}(undef, (Int(worksize_ref[]),)) : ROCArray{Int8}(undef, (1,))
    return handle_ref[], workarea
end

# promote to a complex floating-point type (out-of-place only),
# so implementations only need Complex{Float} methods
for f in (:fft, :bfft, :ifft)
    pf = Symbol("plan_", f)
    @eval begin
        $f(x::ROCArray{<:Real}, region=1:ndims(x)) = $f(complexfloat(x), region)
        $pf(x::ROCArray{<:Real}, region) = $pf(complexfloat(x), region)
        $f(x::ROCArray{<:Complex{<:Union{Integer,Rational}}}, region=1:ndims(x)) = $f(complexfloat(x), region)
        $pf(x::ROCArray{<:Complex{<:Union{Integer,Rational}}}, region) = $pf(complexfloat(x), region)
    end
end
rfft(x::ROCArray{<:Union{Integer,Rational}}, region=1:ndims(x)) = rfft(realfloat(x), region)
plan_rfft(x::ROCArray{<:Real}, region) = plan_rfft(realfloat(x), region)

# inplace/notinplace complex
for (f,xtype,inplace,forward) in ((:plan_fft!, :rocfft_transform_type_complex_forward, :true, :true),
                          (:plan_bfft!, :rocfft_transform_type_complex_inverse, :true, :false),
                          (:plan_fft, :rocfft_transform_type_complex_forward, :false, :true),
                          (:plan_bfft, :rocfft_transform_type_complex_inverse, :false, :false))
    @eval begin
        function $f(X::ROCArray{T,N}, region) where {T<:rocfftComplexes,N}
            _inplace = $(inplace)
            _xtype = $(xtype)
            K = $(forward)
            pp = create_plan(_xtype, size(X), T, _inplace, region)
            return cROCFFTPlan{T,K,_inplace,N}(pp..., X, size(X), _xtype, region)
        end
    end # eval
end

function plan_rfft(X::ROCArray{T,N}, region) where {T<:rocfftReals,N}
    inplace = false
    xtype = rocfft_transform_type_real_forward
    pp = create_plan(xtype, size(X), T, inplace, region)
    ydims = collect(size(X))
    ydims[region[1]] = div(ydims[region[1]],2) + 1
    return rROCFFTPlan{T,ROCFFT_FORWARD,inplace,N}(pp..., X, (ydims...,), xtype, region)
end

function plan_brfft(X::ROCArray{T,N}, d::Integer, region::Any) where {T<:rocfftComplexes,N}
    inplace = false
    xtype = rocfft_transform_type_real_inverse
    ydims = collect(size(X))
    ydims[region[1]] = d
    pp = create_plan(xtype, (ydims...,), T, inplace, region)
    return rROCFFTPlan{T,ROCFFT_INVERSE,inplace,N}(pp..., X, (ydims...,), xtype, region)
end

# FIXME: plan_inv methods allocate needlessly (to provide type parameters and normalization function)
# Perhaps use FakeArray types to avoid this.
function plan_inv(p::cROCFFTPlan{T,ROCFFT_FORWARD,inplace,N}) where {T<:rocfftComplexes,N,inplace}
    X = ROCArray{T}(undef, p.sz)
    xtype = rocfft_transform_type_complex_inverse
    pp = create_plan(xtype, p.sz, T, inplace, p.region)
    ScaledPlan(cROCFFTPlan{T,ROCFFT_INVERSE,inplace,N}(pp..., X, p.sz, xtype, p.region),
               normalization(X, p.region))
end

function plan_inv(p::cROCFFTPlan{T,ROCFFT_INVERSE,inplace,N}) where {T<:rocfftComplexes,N,inplace}
    X = ROCArray{T}(undef, p.sz)
    xtype = rocfft_transform_type_complex_forward
    pp = create_plan(xtype, p.sz, T, inplace, p.region)
    ScaledPlan(cROCFFTPlan{T,ROCFFT_FORWARD,inplace,N}(pp..., X, p.sz, xtype, p.region),
               normalization(X, p.region))
end

function plan_inv(p::rROCFFTPlan{T,ROCFFT_FORWARD,inplace,N}) where {T<:rocfftReals,N,inplace}
    X = ROCArray{complex(T)}(undef, p.osz)
    Y = ROCArray{T}(undef, p.sz)
    xtype = rocfft_transform_type_real_inverse
    pp = create_plan(xtype, p.sz, T, inplace, p.region)
    scale = normalization(Y, p.region)
    ScaledPlan(rROCFFTPlan{complex(T),ROCFFT_INVERSE,inplace,N}(pp..., X, p.sz, xtype, p.region), scale)
end

function plan_inv(p::rROCFFTPlan{T,ROCFFT_INVERSE,inplace,N}) where {T<:rocfftComplexes,N,inplace}
    X = ROCArray{real(T)}(undef, p.osz)
    xtype = rocfft_transform_type_real_forward
    pp = create_plan(xtype, p.osz, T, inplace, p.region)
    scale = normalization(X, p.region)
    ScaledPlan(rROCFFTPlan{real(T),ROCFFT_FORWARD,inplace,N}(pp..., X, p.sz, xtype, p.region), scale)
end

## plan execution

function assert_applicable(p::ROCFFTPlan{T}, X::ROCArray{T}) where T
    (size(X) == p.sz) ||
        throw(ArgumentError("rocFFT plan applied to wrong-size input"))
end

function assert_applicable(p::ROCFFTPlan{T,K}, X::ROCArray{T}, Y::ROCArray{Ty}) where {T,K,Ty}
    assert_applicable(p, X)
    (size(Y) == p.osz) ||
        throw(ArgumentError("RocFFT plan applied to wrong-size output"))
    # type errors should be impossible by dispatch, but just in case:
    if p.xtype == rocfft_transform_type_real_inverse
        (Ty == real(T)) ||
            throw(ArgumentError("Type mismatch for argument Y"))
    elseif p.xtype == rocfft_transform_type_real_forward
        (Ty == complex(T)) ||
            throw(ArgumentError("Type mismatch for argument Y"))
    else
        (Ty == T) ||
            throw(ArgumentError("Type mismatch for argument Y"))
    end
end

function unsafe_execute!(plan::cROCFFTPlan{T,K,true,N}, X::ROCArray{T,N}) where {T,K,N}
    rocfft_execute(plan, [pointer(X),], C_NULL, plan.execution_info)
end

function unsafe_execute!(plan::cROCFFTPlan{T,K,false,N}, X::ROCArray{T,N}, Y::ROCArray{T}) where {T,N,K}
    Xcopy = copy(X) # since input array can also be modified
    rocfft_execute(plan, [pointer(Xcopy),], [pointer(Y),], plan.execution_info)
end

function unsafe_execute!(plan::rROCFFTPlan{T,ROCFFT_FORWARD,false,N}, X::ROCArray{T,N}, Y::ROCArray{<:rocfftComplexes,N}) where {T<:rocfftReals,N}
    @assert plan.xtype == rocfft_transform_type_real_forward
    Xcopy = copy(X)
    rocfft_execute(plan, [pointer(Xcopy),], [pointer(Y),], plan.execution_info)
end

function unsafe_execute!(plan::rROCFFTPlan{T,ROCFFT_INVERSE,false,N}, X::ROCArray{T,N}, Y::ROCArray{<:rocfftReals,N}) where {T<:rocfftComplexes,N}
    @assert plan.xtype == rocfft_transform_type_real_inverse
    Xcopy = copy(X)
    rocfft_execute(plan, [pointer(Xcopy),], [pointer(Y),], plan.execution_info)
end

function LinearAlgebra.mul!(y::ROCArray{Ty}, p::ROCFFTPlan{T,K,false}, x::ROCArray{T}) where {T,Ty,K}
    assert_applicable(p, x, y)
    unsafe_execute!(p, x, y)
    # 0 is default stream, TODO: make this more scalable
    hipStreamSynchronize(Ptr{Cvoid}(UInt64(0)))
    return y
end

function Base.:(*)(p::cROCFFTPlan{T,K,true,N}, x::ROCArray{T,N}) where {T,K,N}
    assert_applicable(p, x)
    unsafe_execute!(p, x)
    hipStreamSynchronize(Ptr{Cvoid}(UInt64(0)))
    x
end

function Base.:(*)(p::cROCFFTPlan{T,K,false,N}, x::ROCArray{T,N}) where {T,K,N}
    y = ROCArray{T}(undef, p.osz)
    mul!(y, p, x)
end

function Base.:(*)(p::rROCFFTPlan{T,ROCFFT_FORWARD,false,N}, x::ROCArray{T,N}) where {T<:rocfftReals,N}
    @assert p.xtype == rocfft_transform_type_real_forward
    y = ROCArray{complex(T)}(undef, p.osz)
    mul!(y, p, x)
end

function Base.:(*)(p::rROCFFTPlan{T,ROCFFT_INVERSE,false,N}, x::ROCArray{T,N}) where {T<:rocfftComplexes,N}
    @assert p.xtype == rocfft_transform_type_real_inverse
    y = ROCArray{real(T)}(undef, p.osz)
    mul!(y, p, x)
end


text:
using DirectSearch
using Test
using REPL

# include("./test_OrthoMADS.jl")
# include("./test_LTMADS.jl")
# include("./test_Cache.jl")
# include("./test_Search.jl")
# include("./test_Constraints.jl")
# include("./test_Core.jl")
# include("./test_StoppingCondition.jl")
# include("./test_Report.jl")
# include("./solve_tests.jl")
include("./test_BiMADS.jl")


text:
export MatrixDirichlet

import SpecialFunctions: digamma, loggamma
import Base: eltype

struct MatrixDirichlet{T <: Real, A <: AbstractMatrix{T} } <: ContinuousMatrixDistribution
    a :: A
end

Distributions.mean(dist::MatrixDirichlet) = dist.a ./ sum(dist.a, dims = 1)

Base.eltype(::MatrixDirichlet{T}) where T = T

vague(::Type{ <: MatrixDirichlet }, dims::Int)              = MatrixDirichlet(ones(dims, dims))
vague(::Type{ <: MatrixDirichlet }, dims1::Int, dims2::Int) = MatrixDirichlet(ones(dims1, dims2))
vague(::Type{ <: MatrixDirichlet }, dims::Tuple)            = MatrixDirichlet(ones(dims))

function Distributions.entropy(dist::MatrixDirichlet)
    return mapreduce(+, eachcol(dist.a)) do column
        scolumn = sum(column)
        -sum((column .- 1.0) .* (digamma.(column) .- digamma.(scolumn))) - loggamma(scolumn) + sum(loggamma.(column))
    end
end

mean(::typeof(log), dist::MatrixDirichlet) = digamma.(dist.a) .- digamma.(sum(dist.a, dims = 1))

prod_analytical_rule(::Type{ <: MatrixDirichlet }, ::Type{ <: MatrixDirichlet }) = ProdAnalyticalRuleAvailable()

function prod(::ProdAnalytical, left::MatrixDirichlet, right::MatrixDirichlet)
    T = promote_type(eltype(left), eltype(right))
    return MatrixDirichlet(left.a + right.a .- one(T))
end


text:
"""
    CLIMA driver configurations

Use to run CLIMA using the CLIMA driver's `CLIMA.invoke()`. User-customized
configurations can use these as templates.
"""

using MPI

using ..Atmos
using ..HydrostaticBoussinesq
using ..DGmethods
using ..DGmethods.NumericalFluxes
using ..Mesh.Topologies
using ..Mesh.Grids
using ..ODESolvers
using ..PlanetParameters

abstract type AbstractSolverType end
struct ExplicitSolverType <: AbstractSolverType
    solver_method::Function
    ExplicitSolverType(;solver_method=LSRK54CarpenterKennedy) = new(solver_method)
end
struct IMEXSolverType <: AbstractSolverType
    linear_model::Type
    linear_solver::Type
    solver_method::Function
    function IMEXSolverType(;linear_model=AtmosAcousticGravityLinearModel,
                            linear_solver=ManyColumnLU,
                            solver_method=ARK2GiraldoKellyConstantinescu)
        return new(linear_model, linear_solver, solver_method)
    end
end
DefaultSolverType = IMEXSolverType

struct DriverConfiguration{FT}
    name::String
    N::Int
    array_type
    solver_type::AbstractSolverType
    #
    # AtmosModel details
    bl::BalanceLaw
    #
    # execution details
    mpicomm::MPI.Comm
    #
    # mesh details
    grid::DiscontinuousSpectralElementGrid
    #
    # DGModel details
    numfluxnondiff::NumericalFluxNonDiffusive
    numfluxdiff::NumericalFluxDiffusive
    gradnumflux::NumericalFluxGradient

    function DriverConfiguration(name::String, N::Int, FT, array_type,
                                 solver_type::AbstractSolverType,
                                 bl::BalanceLaw,
                                 mpicomm::MPI.Comm,
                                 grid::DiscontinuousSpectralElementGrid,
                                 numfluxnondiff::NumericalFluxNonDiffusive,
                                 numfluxdiff::NumericalFluxDiffusive,
                                 gradnumflux::NumericalFluxGradient)
        return new{FT}(name, N, array_type, solver_type, bl, mpicomm, grid,
                       numfluxnondiff, numfluxdiff, gradnumflux)
    end
end

function Atmos_LES_Configuration(
        name::String,
        N::Int,
        (Δx, Δy, Δz)::NTuple{3,FT},
        xmax::Int, ymax::Int, zmax::Int,
        init_LES!;
        xmin           = 0,
        ymin           = 0,
        zmin           = 0,
        array_type     = CLIMA.array_type(),
        solver_type    = IMEXSolverType(linear_solver=SingleColumnLU),
        model          = AtmosModel{FT}(AtmosLESConfiguration;
                                        init_state=init_LES!),
        mpicomm        = MPI.COMM_WORLD,
        boundary       = ((0,0), (0,0), (1,2)),
        periodicity    = (true, true, false),
        meshwarp       = (x...)->identity(x),
        numfluxnondiff = Rusanov(),
        numfluxdiff    = CentralNumericalFluxDiffusive(),
        gradnumflux    = CentralNumericalFluxGradient()
    ) where {FT<:AbstractFloat}

    @info @sprintf("""Establishing Atmos LES configuration for %s
                   precision        = %s
                   polynomial order = %d
                   grid             = %dx%dx%d
                   resolution       = %dx%dx%d
                   MPI ranks        = %d""",
                   name, FT, N,
                   xmax, ymax, zmax,
                   Δx, Δy, Δz,
                   MPI.Comm_size(mpicomm))

    brickrange = (grid1d(xmin, xmax, elemsize=Δx*N),
                  grid1d(ymin, ymax, elemsize=Δy*N),
                  grid1d(zmin, zmax, elemsize=Δz*N))
    topology = StackedBrickTopology(mpicomm, brickrange,
                                    periodicity=periodicity,
                                    boundary=boundary)

    grid = DiscontinuousSpectralElementGrid(topology,
                                            FloatType=FT,
                                            DeviceArray=array_type,
                                            polynomialorder=N,
                                            meshwarp=meshwarp)

    return DriverConfiguration(name, N, FT, array_type, solver_type, model,
                               mpicomm, grid, numfluxnondiff, numfluxdiff,
                               gradnumflux)
end

function Atmos_GCM_Configuration(
        name::String,
        N::Int,
        (nelem_horz, nelem_vert)::NTuple{2,Int},
        domain_height::FT,
        init_GCM!;
        array_type         = CLIMA.array_type(),
        solver_type        = DefaultSolverType(),
        model              = AtmosModel{FT}(AtmosGCMConfiguration;
                                             init_state=init_GCM!),
        mpicomm            = MPI.COMM_WORLD,
        meshwarp::Function = cubedshellwarp,
        numfluxnondiff     = Rusanov(),
        numfluxdiff        = CentralNumericalFluxDiffusive(),
        gradnumflux        = CentralNumericalFluxGradient()
    ) where {FT<:AbstractFloat}
    @info @sprintf("""Establishing Atmos GCM configuration for %s
                   precision        = %s
                   polynomial order = %d
                   #horiz elems     = %d
                   #vert_elems      = %d
                   domain height    = %.2e
                   MPI ranks        = %d""",
                   name, FT, N,
                   nelem_horz, nelem_vert, domain_height,
                   MPI.Comm_size(mpicomm))

    vert_range = grid1d(FT(planet_radius), FT(planet_radius+domain_height), nelem=nelem_vert)

    topology = StackedCubedSphereTopology(mpicomm, nelem_horz, vert_range)

    grid = DiscontinuousSpectralElementGrid(topology,
                                            FloatType=FT,
                                            DeviceArray=array_type,
                                            polynomialorder=N,
                                            meshwarp=meshwarp)

    return DriverConfiguration(name, N, FT, array_type, solver_type, model,
                               mpicomm, grid, numfluxnondiff, numfluxdiff,
                               gradnumflux)
end

function Ocean_BoxGCM_Configuration(
    name::String,
    N::Int,
    (Nˣ, Nʸ, Nᶻ)::NTuple{3,Int},
    model::HydrostaticBoussinesqModel;
    FT             = Float64,
    array_type     = CLIMA.array_type(),
    solver_type    = ExplicitSolverType(solver_method=LSRK144NiegemannDiehlBusch),
    mpicomm        = MPI.COMM_WORLD,
    numfluxnondiff = Rusanov(),
    numfluxdiff    = CentralNumericalFluxDiffusive(),
    gradnumflux    = CentralNumericalFluxGradient(),
    periodicity    = (false, false, false),
    boundary       = ((1, 1), (1, 1), (2, 3))
    )

    brickrange = (range(FT(0);  length=Nˣ+1, stop=model.problem.Lˣ),
                  range(FT(0);  length=Nʸ+1, stop=model.problem.Lʸ),
                  range(FT(-model.problem.H); length=Nᶻ+1, stop=0))

    topology = StackedBrickTopology(mpicomm, brickrange;
                                    periodicity = periodicity,
                                    boundary = boundary)

    grid = DiscontinuousSpectralElementGrid(topology,
                                            FloatType = FT,
                                            DeviceArray = array_type,
                                            polynomialorder = N)

    return DriverConfiguration(name, N, FT, array_type, solver_type, model,
                               mpicomm, grid, numfluxnondiff, numfluxdiff,
                               gradnumflux)
end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
using CSV
using JuMP
using CPLEX
using DataFrames
using Random
#import Pkg;Pkg.add("CSV"); Pkg.add("JuMP");Pkg.add("CPLEX");Pkg.add("DataFrames");Pkg.add("Random");

include("functions.jl")

dataSet = "flowers"
#dataSet="kidney"
dataFolder = "C:\\Users\\thibaut\\Documents\\Ponts_3A\\CNAM\\RODM\\data\\"
#dataFolder = "../data/"
resultsFolder = "C:\\Users\\thibaut\\Documents\\Ponts_3A\\CNAM\\RODM\\res\\"
#resultsFolder = "../res/"

rawDataPath = dataFolder * dataSet * ".csv"
#println(isfile(rawDataPath))
#print(cd(readdir,dataFolder))

#modelTest= Model(with_optimizer(CPLEX.Optimizer))
println("début prgm")
train, test = createFeatures(dataFolder, dataSet)

rulesMain = createRules(dataSet, resultsFolder, train)
println("=============Règles trouvées==========")
#println("rulesMain",rulesMain)

println("=============Début Tri==========")
timeLimitInSeconds = 6000
orderedRules = sortRules(dataSet, resultsFolder, train, rulesMain, timeLimitInSeconds)
println("=============Fin Tri==========")

println("-- Train results")
goodShowStatistics(orderedRules, train)
#badShowStatistics(orderedRules, train)
println("-- Test results")
goodShowStatistics(orderedRules, test)
#badShowStatistics(orderedRules, test)

println("ordered rules :\n", orderedRules)


text:
struct Identity <: LocalOperator
end

kernelvals(biop::Identity, x) = nothing
integrand(op::Identity, kernel, x, g, f) = dot(f[1], g[1])
scalartype(op::Identity) = Union{}

struct NCross <: LocalOperator
end

kernelvals(op::NCross, mp) = nothing
integrand(op::NCross, kernel, x, g, f) = dot(g[1], normal(x) × f[1])
scalartype(op::NCross) = Union{}

function _alloc_workspace(qd, g, f, tels, bels)
    q = qd[1]
    τ = tels[1]
    w, p = q[1], neighborhood(τ,q[2])
    a = (w, p, g(p), f(p))
    A = Vector{typeof(a)}(undef,length(qd))
end

# defaultquadstrat(op::LocalOperator, tfs, bfs) = defaultquadstrat(op, refspace(tfs), refspace(bfs))
# defaultquadstrat(op::LocalOperator, tfs::Space, bfs::DirectProductSpace) = defaultquadstrat(op, tfs, bfs.factors[1])
# defaultquadstrat(op::LocalOperator, tfs::DirectProductSpace, bfs::Space) = defaultquadstrat(op, tfs.factors[1], bfs)
# defaultquadstrat(op::LocalOperator, tfs::DirectProductSpace, bfs::DirectProductSpace) = defaultquadstrat(op, tfs.factors[1], bfs.factors[1])

const LinearRefSpaceTriangle = Union{RTRefSpace, NDRefSpace, BDMRefSpace}
defaultquadstrat(::LocalOperator, ::LinearRefSpaceTriangle, ::LinearRefSpaceTriangle) = SingleNumQStrat(6)
function quaddata(op::LocalOperator, g::LinearRefSpaceTriangle, f::LinearRefSpaceTriangle, tels, bels,
        qs::SingleNumQStrat)

    u, w = trgauss(qs.quad_rule)
    qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)]
    A = _alloc_workspace(qd, g, f, tels, bels)

    return qd, A
end

defaultquadstrat(::LocalOperator, ::subReferenceSpace, ::subReferenceSpace) = SingleNumQStrat(6)
function quaddata(op::LocalOperator, g::subReferenceSpace, f::subReferenceSpace, tels, bels,
        qs::SingleNumQStrat)

    u, w = trgauss(qs.quad_rule)
    qd = [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)]
    A = _alloc_workspace(qd, g, f, tels, bels)
    return qd, A
end

const LinearRefSpaceTetr = Union{NDLCCRefSpace, NDLCDRefSpace}
defaultquadstrat(::LocalOperator, ::LinearRefSpaceTetr, ::LinearRefSpaceTetr) = SingleNumQStrat(3)
function quaddata(op::LocalOperator, g::LinearRefSpaceTetr, f::LinearRefSpaceTetr, tels, bels, qs::SingleNumQStrat)

    o, x, y, z = CompScienceMeshes.euclidianbasis(3)
    reftet = simplex(x,y,z,o)
    qps = quadpoints(reftet, qs.quad_rule)
    qd = [(w, parametric(p)) for (p,w) in qps]
    A = _alloc_workspace(qd, g, f, tels, bels)
    return qd, A
end

defaultquadstrat(::LocalOperator, ::LagrangeRefSpace{T,D1,2}, ::LagrangeRefSpace{T,D2,2}) where {T,D1,D2} = SingleNumQStrat(6)
function quaddata(op::LocalOperator, g::LagrangeRefSpace{T,Deg,2} where {T,Deg},
    f::LagrangeRefSpace, tels::Vector, bels::Vector, qs::SingleNumQStrat)

    u, w = legendre(qs.quad_rule, 0.0, 1.0)
    qd = [(w[i],u[i]) for i in eachindex(w)]
    A = _alloc_workspace(qd, g, f, tels, bels)
    return qd, A
end

defaultquadstrat(::LocalOperator, ::LagrangeRefSpace{T,D1,3}, ::LagrangeRefSpace{T,D2,3}) where {T,D1,D2} = SingleNumQStrat(6)
function quaddata(op::LocalOperator, g::LagrangeRefSpace{T,Deg,3} where {T,Deg},
    f::LagrangeRefSpace, tels::Vector, bels::Vector, qs::SingleNumQStrat)

    u, w = trgauss(qs.quad_rule)
    qd = [(w[i], SVector(u[1,i], u[2,i])) for i in 1:length(w)]
    A = _alloc_workspace(qd, g, f, tels, bels)
    return qd, A
end

defaultquadstrat(::LocalOperator, ::LagrangeRefSpace{T,D1,4}, ::LagrangeRefSpace{T,D2,4}) where {T,D1,D2} = SingleNumQStrat(6)
function quaddata(op::LocalOperator, g::LagrangeRefSpace{T,Deg,4} where {T,Deg},
    f::LagrangeRefSpace, tels::Vector, bels::Vector, qs::SingleNumQStrat)

     o, x, y, z = CompScienceMeshes.euclidianbasis(3)
     reftet = simplex(x,y,z,o)
     qps = quadpoints(reftet, qs.quad_rule)
     qd = [(w, parametric(p)) for (p,w) in qps]
     A = _alloc_workspace(qd, g, f, tels, bels)
     return qd, A
end

function quadrule(op::LocalOperator, ψ::RefSpace, ϕ::RefSpace, τ, (qd,A), qs::SingleNumQStrat)
    for i in eachindex(qd)
        q = qd[i]
        w, p = q[1], neighborhood(τ,q[2])
        A[i] = (w, p, ψ(p), ϕ(p))
    end
    return A
end


text:
using DiffResults, StaticArrays, Test

using DiffResults: DiffResult, GradientResult, JacobianResult, HessianResult,
                   value, value!, derivative, derivative!, gradient, gradient!,
                   jacobian, jacobian!, hessian, hessian!

@testset "DiffResult" begin

    k = 4
    n0, n1, n2 = rand(), rand(), rand()
    x0, x1, x2 = rand(k), rand(k, k), rand(k, k, k)
    s0, s1, s2 = SVector{k}(rand(k)), SMatrix{k,k}(rand(k, k)), SArray{Tuple{k,k,k}}(rand(k, k, k))
    rn = DiffResult(n0, n1, n2)
    rx = DiffResult(x0, x1, x2)
    rs = DiffResult(s0, s1, s2)
    rsmix = DiffResult(n0, s0, s1)

    issimilar(x, y) = typeof(x) == typeof(y) && size(x) == size(y)
    issimilar(x::DiffResult, y::DiffResult) = issimilar(value(x), value(y)) && all(issimilar, zip(x.derivs, y.derivs))
    issimilar(t::Tuple) = issimilar(t...)

    @test rn === DiffResult(n0, n1, n2)
    @test rx == DiffResult(x0, x1, x2)
    @test rs === DiffResult(s0, s1, s2)
    @test rsmix === DiffResult(n0, s0, s1)

    @test issimilar(GradientResult(x0), DiffResult(first(x0), x0))
    @test issimilar(JacobianResult(x0), DiffResult(x0, similar(x0, k, k)))
    @test issimilar(JacobianResult(similar(x0, k + 1), x0), DiffResult(similar(x0, k + 1), similar(x0, k + 1, k)))
    @test issimilar(HessianResult(x0), DiffResult(first(x0), x0, similar(x0, k, k)))

    @test GradientResult(s0) === DiffResult(first(s0), s0)
    @test JacobianResult(s0) === DiffResult(s0, zeros(SMatrix{k,k,Float64}))
    @test JacobianResult(SVector{k+1}(vcat(s0, 0.0)), s0) === DiffResult(SVector{k+1}(vcat(s0, 0.0)), zeros(SMatrix{k+1,k,Float64}))
    @test HessianResult(s0) === DiffResult(first(s0), s0, zeros(SMatrix{k,k,Float64}))

    @test eltype(rn) === typeof(n0)
    @test eltype(rx) === eltype(x0)
    @test eltype(rs) === eltype(s0)

    rn_copy = copy(rn)
    @test rn == rn_copy
    @test rn === rn_copy

    rx_copy = copy(rx)
    @test rx == rx_copy
    @test rx !== rx_copy

    rs_copy = copy(rs)
    @test rs == rs_copy
    @test rs === rs_copy

    rsmix_copy = copy(rsmix)
    @test rsmix == rsmix_copy
    @test rsmix === rsmix_copy

    @testset "value/value!" begin

        @test value(rn) === n0
        @test value(rx) === x0
        @test value(rs) === s0
        @test value(rsmix) === n0

        rn = value!(rn, n1)
        @test value(rn) === n1
        rn = value!(rn, n0)

        x0_new, x0_copy = rand(k), copy(x0)
        rx = value!(rx, x0_new)
        @test value(rx) === x0 == x0_new
        rx = value!(rx, x0_copy)

        s0_new = rand(k)
        rs = value!(rs, s0_new)
        @test value(rs) == s0_new
        @test typeof(value(rs)) === typeof(s0)
        rs = value!(rs, s0)

        rsmix = value!(rsmix, n1)
        @test value(rsmix) === n1
        rsmix = value!(rsmix, n0)

        rn = value!(exp, rn, n1)
        @test value(rn) === exp(n1)
        rn = value!(rn, n0)

        x0_new, x0_copy = rand(k), copy(x0)
        rx = value!(exp, rx, x0_new)
        @test value(rx) === x0 == exp.(x0_new)
        rx = value!(rx, x0_copy)

        s0_new = rand(k)
        rs = value!(exp, rs, s0_new)
        @test value(rs) == exp.(s0_new)
        @test typeof(value(rs)) === typeof(s0)
        rs = value!(rs, s0)

        rsmix = value!(exp, rsmix, n1)
        @test value(rsmix) === exp(n1)
        rsmix = value!(rsmix, n0)

        ksqrt = Int(sqrt(k))
        T = typeof(SMatrix{ksqrt,ksqrt}(rand(ksqrt,ksqrt)))
        rs_new = value!(rs, convert(T, value(rs)))
        @test rs_new === rs
    end

    @testset "derivative/derivative!" begin

        @test derivative(rn) === n1
        @test derivative(rn, Val{2}) === n2

        @test derivative(rx) === x1
        @test derivative(rx, Val{2}) === x2

        @test derivative(rs) === s1
        @test derivative(rs, Val{2}) === s2

        @test derivative(rsmix) === s0
        @test derivative(rsmix, Val{2}) === s1

        rn = derivative!(rn, n0)
        @test derivative(rn) === n0
        rn = derivative!(rn, n1)

        x1_new, x1_copy = rand(k, k), copy(x1)
        rx = derivative!(rx, x1_new)
        @test derivative(rx) === x1 == x1_new
        rx = derivative!(rx, x1_copy)

        s1_new = rand(k, k)
        rs = derivative!(rs, s1_new)
        @test derivative(rs) == s1_new
        @test typeof(derivative(rs)) === typeof(s1)
        rs = derivative!(rs, s1)

        s0_new = rand(k)
        rsmix = derivative!(rsmix, s0_new)
        @test derivative(rsmix) == s0_new
        @test typeof(derivative(rsmix)) === typeof(s0)
        rsmix = derivative!(rsmix, s0)

        rn = derivative!(rn, n1, Val{2})
        @test derivative(rn, Val{2}) === n1
        rn = derivative!(rn, n2, Val{2})

        x2_new, x2_copy = rand(k, k, k), copy(x2)
        rx = derivative!(rx, x2_new, Val{2})
        @test derivative(rx, Val{2}) === x2 == x2_new
        rx = derivative!(rx, x2_copy, Val{2})

        s2_new = rand(k, k, k)
        rs = derivative!(rs, s2_new, Val{2})
        @test derivative(rs, Val{2}) == s2_new
        @test typeof(derivative(rs, Val{2})) === typeof(s2)
        rs = derivative!(rs, s2, Val{2})

        s1_new = rand(k, k)
        rsmix = derivative!(rsmix, s1_new, Val{2})
        @test derivative(rsmix, Val{2}) == s1_new
        @test typeof(derivative(rsmix, Val{2})) === typeof(s1)
        rsmix = derivative!(rsmix, s1, Val{2})

        rn = derivative!(exp, rn, n0)
        @test derivative(rn) === exp(n0)
        rn = derivative!(rn, n1)

        x1_new, x1_copy = rand(k, k), copy(x1)
        rx = derivative!(exp, rx, x1_new)
        @test derivative(rx) === x1 == exp.(x1_new)
        rx = derivative!(exp, rx, x1_copy)

        s1_new = rand(k, k)
        rs = derivative!(exp, rs, s1_new)
        @test derivative(rs) == exp.(s1_new)
        @test typeof(derivative(rs)) === typeof(s1)
        rs = derivative!(exp, rs, s1)

        s0_new = rand(k)
        rsmix = derivative!(exp, rsmix, s0_new)
        @test derivative(rsmix) == exp.(s0_new)
        @test typeof(derivative(rsmix)) === typeof(s0)
        rsmix = derivative!(exp, rsmix, s0)

        rn = derivative!(exp, rn, n1, Val{2})
        @test derivative(rn, Val{2}) === exp(n1)
        rn = derivative!(rn, n2, Val{2})

        x2_new, x2_copy = rand(k, k, k), copy(x2)
        rx = derivative!(exp, rx, x2_new, Val{2})
        @test derivative(rx, Val{2}) === x2 == exp.(x2_new)
        rx = derivative!(exp, rx, x2_copy, Val{2})

        s2_new = rand(k, k, k)
        rs = derivative!(exp, rs, s2_new, Val{2})
        @test derivative(rs, Val{2}) == exp.(s2_new)
        @test typeof(derivative(rs, Val{2})) === typeof(s2)
        rs = derivative!(exp, rs, s2, Val{2})

        s1_new = rand(k, k)
        rsmix = derivative!(exp, rsmix, s1_new, Val{2})
        @test derivative(rsmix, Val{2}) == exp.(s1_new)
        @test typeof(derivative(rsmix, Val{2})) === typeof(s1)
        rsmix = derivative!(exp, rsmix, s1, Val{2})
    end

    @testset "gradient/gradient!" begin

        x1_new, x1_copy = rand(k, k), copy(x1)
        rx = gradient!(rx, x1_new)
        @test gradient(rx) === x1 == x1_new
        rx = gradient!(rx, x1_copy)

        s1_new = rand(k, k)
        rs = gradient!(rs, s1_new)
        @test gradient(rs) == s1_new
        @test typeof(gradient(rs)) === typeof(s1)
        rs = gradient!(rs, s1)

        x1_new, x1_copy = rand(k, k), copy(x1)
        rx = gradient!(exp, rx, x1_new)
        @test gradient(rx) === x1 == exp.(x1_new)
        rx = gradient!(exp, rx, x1_copy)

        s0_new = rand(k)
        rsmix = gradient!(exp, rsmix, s0_new)
        @test gradient(rsmix) == exp.(s0_new)
        @test typeof(gradient(rsmix)) === typeof(s0)
        rsmix = gradient!(exp, rsmix, s0)

        T = typeof(SVector{k*k}(rand(k*k)))
        rs_new = gradient!(rs, convert(T, gradient(rs)))
        @test rs_new === rs
    end

    @testset "jacobian/jacobian!"  begin

        x1_new, x1_copy = rand(k, k), copy(x1)
        rx = jacobian!(rx, x1_new)
        @test jacobian(rx) === x1 == x1_new
        rx = jacobian!(rx, x1_copy)

        s1_new = rand(k, k)
        rs = jacobian!(rs, s1_new)
        @test jacobian(rs) == s1_new
        @test typeof(jacobian(rs)) === typeof(s1)
        rs = jacobian!(rs, s1)

        x1_new, x1_copy = rand(k, k), copy(x1)
        rx = jacobian!(exp, rx, x1_new)
        @test jacobian(rx) === x1 == exp.(x1_new)
        rx = jacobian!(exp, rx, x1_copy)

        s0_new = rand(k)
        rsmix = jacobian!(exp, rsmix, s0_new)
        @test jacobian(rsmix) == exp.(s0_new)
        @test typeof(jacobian(rsmix)) === typeof(s0)
        rsmix = jacobian!(exp, rsmix, s0)

        T = typeof(SVector{k*k}(rand(k*k)))
        rs_new = jacobian!(rs, convert(T, jacobian(rs)))
        @test rs_new === rs
    end

    @testset "hessian/hessian!" begin

        x2_new, x2_copy = rand(k, k, k), copy(x2)
        rx = hessian!(rx, x2_new)
        @test hessian(rx) === x2 == x2_new
        rx = hessian!(rx, x2_copy)

        s2_new = rand(k, k, k)
        rs = hessian!(rs, s2_new)
        @test hessian(rs) == s2_new
        @test typeof(hessian(rs)) === typeof(s2)
        rs = hessian!(rs, s2)

        x2_new, x2_copy = rand(k, k, k), copy(x2)
        rx = hessian!(exp, rx, x2_new)
        @test hessian(rx) === x2 == exp.(x2_new)
        rx = hessian!(exp, rx, x2_copy)

        s1_new = rand(k, k)
        rsmix = hessian!(exp, rsmix, s1_new)
        @test hessian(rsmix) == exp.(s1_new)
        @test typeof(hessian(rsmix)) === typeof(s1)
        rsmix = hessian!(exp, rsmix, s1)

        T = typeof(SVector{k*k*k}(rand(k*k*k)))
        rs_new = hessian!(rs, convert(T, hessian(rs)))
        @test rs_new === rs

        @test size(gradient(HessianResult(x0))) == size(x0)
        @test size(gradient(HessianResult(x1))) == size(x1)
        @test size(gradient(HessianResult(x2))) == size(x2)

        @test HessianResult(Float32.(x0)).derivs[1] isa Vector{Float32}
    end
end


text:
mutable struct Horizontal <: Geometry
    label::String
    shape::Vector{Box}
    nsample::Int
end

function Horizontal( ; label, shape, nsample)
    return Horizontal(label, shape, nsample)
end


text:
const TZS = Set(timezone_names())
const GMTTZ = TimeZone("GMT")

function toepoch(ts)
    fixedts = astimezone(ZonedDateTime(ts, localzone()), GMTTZ) |> DateTime
    Dates.value(fixedts) - Dates.UNIXEPOCH
end
curts() = toepoch(Dates.now())

struct Opts
    token::String
    host::String
    port::Int
end
const OPTS = Ref(Opts("", "127.0.0.1", 9174))
read_port(port::Integer) = port
read_port(port) = parse(Int, port)

function setupbot!(; token = OPTS[].token,
                     host = OPTS[].host,
                     port = OPTS[].port,
                     email = "",
                     apikey = "",
                     ep = "")
    OPTS[] = Opts(token, host, read_port(port))
    ZulipClient(email = email, apikey = apikey, ep = ep)
end

########################################
# Processing
########################################

function process(obj::ZulipRequest, db, channel, ts, opts = OPTS[])
    status, resp = validate(obj, opts)
    status || return isempty(resp) ? resp : JSON3.write((; content = resp))

    data = obj.data
    resp = if startswith(data, "help")
        process_help(obj, db, opts)
    elseif startswith(data, "list")
        process_list(obj, db, opts)
    elseif startswith(data, "remove")
        process_remove(obj, db, channel, ts, opts)
    elseif startswith(data, "timezone")
        process_timezone(obj, db, opts)
    else
        process_reminder(obj, db, channel, ts, opts)
    end

    return JSON3.write((; content = resp))
end

########################################
# Server
########################################

function heartbeat!(sched, ts, input, output)
    sorted = true
    lock(input)
    while isready(input)
        sorted = false
        datain = take!(input)
        if datain[2] == 1
            push!(sched, datain[1])
        else
            idx = findfirst(x -> x.id == datain[1].id && x.msg.sender_id == datain[1].msg.sender_id, sched)
            idx === nothing && continue
            deleteat!(sched, idx)
        end
    end
    unlock(input)
    if !sorted
        sort!(sched, by = x -> x.exects, rev = true)
    end
    isempty(sched) && return nothing

    while !isempty(sched)
        if ts >= sched[end].exects
            tmsg = pop!(sched)
            put!(output, tmsg)
        else
            break
        end
    end

    nothing
end

function cron_worker(input, output, sched = TimedMessage[], sleepduration = 1)
    while true
        sleep(sleepduration)
        ts = curts()
        heartbeat!(sched, ts, input, output)
    end
end

function populate(db, input)
    msgs = select(db, Vector{TimedMessage})
    @info "populating $(length(msgs)) messages"
    for msg in msgs
        @debug msg
        put!(input, (msg, 1))
    end
end

function msg_worker(db, input)
    while true
        try
            tmsg = take!(input)
            @debug tmsg
            msg = tmsg.msg
            content = base64decode(msg.content) |> String
            if msg.type == "private"
                resp = sendMessage(type = "private", to = JSON3.write([msg.sender_id]), content = content)
            else
                stream = base64decode(msg.stream) |> String
                topic = base64decode(msg.topic) |> String
                resp = sendMessage(type = "stream", to = stream, topic = topic, content = content)
            end
            delete(db, tmsg)
            @debug resp
        catch err
            @error err
        end
    end
end

function run(db, opts = OPTS[])
    inmsg_channel = Channel{Tuple{TimedMessage, Int}}(Inf)
    outmsg_channel = Channel{TimedMessage}(Inf)
    @async cron_worker(inmsg_channel, outmsg_channel)
    @async msg_worker(db, outmsg_channel)

    host = opts.host
    port = opts.port
    @info "Starting Reminder Bot server on $host:$port"

    populate(db, inmsg_channel)

    HTTP.serve(opts.host, opts.port) do http
        ts = Dates.now()
        obj = String(HTTP.payload(http))
        @debug obj
        obj = JSON3.read(obj, ZulipRequest)
        @debug obj
        resp = process(obj, db, inmsg_channel, ts, opts)

        isempty(resp) || return HTTP.Response(resp)
        return nothing
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule xdm_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("xdm")
JLLWrappers.@generate_main_file("xdm", UUID("4805d9eb-6a63-5fb7-bff2-245b636b81dc"))
end  # module xdm_jll


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule CSPICE_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("CSPICE")
JLLWrappers.@generate_main_file("CSPICE", UUID("07f52509-e9d9-513c-a20d-3b911885bf96"))
end  # module CSPICE_jll


text:
include("BODMAS.jl")
include("Strings.jl")
include("LaTeX.jl")