text:
using Dates
using Statistics
using Test

using MarketData

using TimeSeries

@testset "apply" begin

@testset "time series methods" begin
    @testset "lag takes previous day and timestamps it to next day" begin
        ta = lag(cl)
        @test values(ta)[1] ≈ 111.94 atol=.01
        @test timestamp(ta)[1] == Date(2000, 1, 4)
    end

    @testset "lag accepts other offset values" begin
        ta = lag(cl, 9)
        @test timestamp(ta)[1] == Date(2000, 1, 14)
    end

    @testset "lag operates on 2d arrays" begin
        ta = lag(ohlc, 9)
        @test timestamp(ta)[1] == Date(2000, 1, 14)
    end

    @testset "lag returns 1d from 1d time arrays" begin
        @test ndims(values(lag(cl))) == 1
    end

    @testset "lag returns 2d from 2d time arrays" begin
        @test ndims(values(lag(ohlc))) == 2
    end

    @testset "lead takes next day and timestamps it to current day" begin
        ta = lead(cl)
        @test values(ta)[1] ≈ 102.5 atol=.1
        @test timestamp(ta)[1] == Date(2000,1,3)
    end

    @testset "lead accepts other offset values" begin
        ta = lead(cl, 9)
        @test values(ta)[1] ≈ 100.44 atol=.1
        @test timestamp(ta)[1] == Date(2000,1,3)
    end

    @testset "lead operates on 2d arrays" begin
        @test timestamp(lead(ohlc, 9))[1] == Date(2000,1,3)
    end

    @testset "lead returns 1d from 1d time arrays" begin
        @test ndims(values(lead(cl))) == 1
    end

    @testset "lead returns 2d from 2d time arrays" begin
        @test ndims(values(lead(ohlc))) == 2
    end

    @testset "diff calculates 1st-order differences" begin
        @test timestamp(diff(op))                 == timestamp(diff(op, padding=false))
        @test values(diff(op))                    == values(diff(op, padding=false))
        @test values(diff(op, padding=false))[1]  == values(op[2])[1] .- values(op[1])[1]
        @test values(diff(op, padding=true))[1]    ≡ NaN
        @test values(diff(op, padding=true))[2]   == values(diff(op))[1]
        @test values(diff(op, padding=true))[2]   == values(op[2])[1] .- values(op[1])[1]
    end

    @testset "diff calculates 1st-order differences for multi-column ts" begin
        @test timestamp(diff(ohlc))                   == timestamp(diff(ohlc, padding=false))
        @test values(diff(ohlc))                      == values(diff(ohlc, padding=false))
        @test values(diff(ohlc, padding=false))[1,:]  == values(ohlc)[2, :] .- values(ohlc)[1, :]
        @test values(diff(ohlc, padding=true))[2,:]   == values(diff(ohlc))[1, :]
        @test values(diff(ohlc, padding=true))[2,:]   == values(ohlc)[2, :] .- values(ohlc)[1, :]
        @test all(x -> isnan(x), values(diff(ohlc, padding=true))[1, :])
    end

    @testset "diff calculates 2nd-order differences" begin
        @test timestamp(diff(op, differences=2))                  == timestamp(diff(op, padding=false, differences=2))
        @test timestamp(diff(diff(op)))                           == timestamp(diff(op, padding=false, differences=2))

        @test values(diff(op, differences=2))                     == values(diff(op, padding=false, differences=2))
        @test values(diff(diff(op)))                              == values(diff(op, padding=false, differences=2))

        @test values(diff(op, padding=true, differences=2))[3]    == values(diff(op, differences=2))[1]
        @test values(diff(op, padding=true, differences=2)[2])[1]  ≡ NaN
        @test values(diff(op, padding=true, differences=2)[1])[1]  ≡ NaN
    end

    @testset "diff calculates 2nd-order differences for multi-column ts" begin
        @test timestamp(diff(ohlc, differences=2))                  == timestamp(diff(ohlc, padding=false, differences=2))
        @test timestamp(diff(diff(ohlc)))                           == timestamp(diff(ohlc, padding=false, differences=2))

        @test values(diff(ohlc, differences=2))                     == values(diff(ohlc, padding=false, differences=2))
        @test values(diff(diff(ohlc)))                              == values(diff(ohlc, padding=false, differences=2))
        @test values(diff(ohlc, padding=true, differences=2))[3, :] == values(diff(ohlc, differences=2))[1, :]

        @test all(x -> isnan(x), values(diff(ohlc, padding=true, differences=2))[2,:])
        @test all(x -> isnan(x), values(diff(ohlc, padding=true, differences=2))[1,:])
    end

    @testset "diff calculates 3rd-order differences" begin
        @test diff(op, differences=3) |> timestamp               == diff(op, padding=false, differences=3) |> timestamp
        @test diff(op, differences=3) |> values                  == diff(op, padding=false, differences=3) |> values
        @test diff(diff(diff(op))) |> timestamp                  == diff(op, padding=false, differences=3) |> timestamp
        @test diff(diff(diff(op))) |> values                     == diff(op, padding=false, differences=3) |> values
        @test values(diff(op, padding=true, differences=3))[4]   == values(diff(op, differences=3))[1]
        @test values(diff(op, padding=true, differences=3))[3]    ≡ NaN
        @test values(diff(op, padding=true, differences=3))[2]    ≡ NaN
        @test values(diff(op, padding=true, differences=3))[1]    ≡ NaN
    end

    @testset "diff calculates 3rd-order differences for multi-column ts" begin
        @test diff(ohlc, differences=3) |> timestamp                 == diff(ohlc, padding=false, differences=3) |> timestamp
        @test diff(ohlc, differences=3) |> values                    == diff(ohlc, padding=false, differences=3) |> values
        @test diff(diff(diff(ohlc))) |> timestamp                    == diff(ohlc, padding=false, differences=3) |> timestamp
        @test diff(diff(diff(ohlc))) |> values                       == diff(ohlc, padding=false, differences=3) |> values
        @test values(diff(ohlc, padding=true, differences=3))[4, :]  == values(diff(ohlc, differences=3))[1, :]
        @test all(x -> isnan(x), values(diff(ohlc, padding=true, differences=3))[3,:])
        @test all(x -> isnan(x), values(diff(ohlc, padding=true, differences=3))[2,:])
        @test all(x -> isnan(x), values(diff(ohlc, padding=true, differences=3))[1,:])
    end

    @testset "diff n lag" begin
        let ta = diff(cl, 5)
            ans = cl .- lag(cl, 5)

            @test values(ta)    == values(ans)
            @test timestamp(ta) == timestamp(ans)
        end

        let ta = diff(ohlc, 5)
            ans = ohlc .- lag(ohlc, 5)

            @test values(ta)    == values(ans)
            @test timestamp(ta) == timestamp(ans)
        end
    end  # @testset "diff n lag"

    @testset "simple return value" begin
        @test values(percentchange(cl, :simple)) == values(percentchange(cl))
        @test values(percentchange(cl))          == values(percentchange(cl, padding=false))
        @test values(percentchange(cl))[1]        ≈ (102.5 - 111.94) / 111.94 atol=.01
        @test values(percentchange(ohlc))[1, :]   ≈ (values(ohlc)[2, :] - values(ohlc)[1, :]) ./ values(ohlc)[1, :]
        @test isnan.(values(percentchange(cl, padding=true))[1])
        @test values(percentchange(cl, padding=true))[2]      ≈ (102.5 - 111.94) / 111.94 atol=.01
        @test values(percentchange(ohlc, padding=true))[2, :] ≈ (values(ohlc)[2,:] - values(ohlc)[1,:]) ./ values(ohlc)[1,:]
    end

    @testset "log return value" begin
        @test values(percentchange(cl, :log)) == values(percentchange(cl, :log, padding=false))
        @test values(percentchange(cl, :log))[1]       ≈ log(102.5) - log(111.94) atol=.01
        @test values(percentchange(ohlc, :log))[1, :]  ≈ log.(values(ohlc)[2,:]) .- log.(values(ohlc)[1,:]) atol=.01
        @test isnan.(values(percentchange(cl, :log, padding=true))[1])
        @test values(percentchange(cl, :log, padding=true))[2]      ≈ log(102.5) - log(111.94)
        @test values(percentchange(ohlc, :log, padding=true))[2, :] ≈ log.(values(ohlc)[2,:]) .- log.(values(ohlc)[1,:])
    end

    @testset "moving supplies correct window length" begin
        let
            ta = moving(mean, cl, 10)
            @test values(ta)       == values(moving(mean, cl, 10, padding=false))
            @test timestamp(ta)[1] == Date(2000, 1, 14)
            @test values(ta)[1]     ≈ mean(values(cl)[1:10])
        end
        let
            ta = moving(mean, cl, 10, padding = true)
            @test timestamp(ta)   == timestamp(cl)
            @test values(ta)[1]   ≡ NaN
            @test values(ta)[10] == values(moving(mean, cl, 10))[1]
        end

        let
            ta = moving(mean, ohlc, 10)
            @test values(ta)        == values(moving(mean, ohlc, 10, padding=false))
            @test values(ta)[1, :]'  ≈ mean(values(ohlc)[1:10, :], dims = 1)
        end
        let
            ta = moving(mean, ohlc, 10, padding = true)
            @test all(isnan, values(ta)[1, :])
            @test values(ta)[10, :] == values(moving(mean, ohlc, 10))[1, :]
        end

        @testset "moving with do syntax" begin
            moving(cl, 10) do x
                @test isa(x, Array{Float64, 1})
                x[1]
            end

            moving(ohlc, 10) do x
                @test isa(x, Array{Float64, 1})
                x[1]
            end
        end
    end

    @testset "upto method accumulates" begin
        @test values(upto(sum, cl))[10]  ≈ sum(values(cl)[1:10])
        @test values(upto(mean, cl))[10] ≈ mean(values(cl)[1:10])

        @test timestamp(upto(sum, cl))[10]  == Date(2000, 1, 14)
        @test timestamp(upto(mean, cl))[10] == Date(2000, 1, 14)

        # transpose the upto value output from column to row vector but values are identical
        @test values(upto(sum, ohlc))[10, :]'  ≈ sum(values(ohlc)[1:10, :], dims = 1)
        @test values(upto(mean, ohlc))[10, :]' ≈ mean(values(ohlc)[1:10, :], dims = 1)
    end
end

@testset "basecall works with Base methods" begin
    @testset "cumsum works" begin
        @test values(basecall(cl, cumsum))[2] == values(cl)[1] + values(cl)[2]
    end
end

@testset "adding/removing missing rows works" begin
    uohlc = uniformspace(ohlc[1:8])

    @testset "uniform spacing detection works" begin
        @test uniformspaced(datetime1)
        @test uniformspaced(uohlc)
        @test !uniformspaced(cl)
        @test !uniformspaced(ohlc)
    end

    @testset "forcing uniform spacing works" begin
        @test length(uohlc)    == 10
        @test values(uohlc[5]) == values(ohlc[5])
        @test all(isnan, values(uohlc[6]))
        @test all(isnan, values(uohlc[7]))
        @test values(uohlc[8]) == values(ohlc[6])
    end

    @testset "dropnan works" begin
        nohlc = TimeArray(timestamp(ohlc), copy(values(ohlc)), colnames(ohlc), meta(ohlc))
        values(nohlc)[7:12, 2] .= NaN

        @test timestamp(dropnan(uohlc))       == timestamp(dropnan(uohlc, :all))
        @test values(dropnan(uohlc))          == values(dropnan(uohlc, :all))

        @test values(dropnan(ohlc, :all))     == values(ohlc)
        @test timestamp(dropnan(nohlc, :all)) == timestamp(ohlc)
        @test timestamp(dropnan(uohlc, :all)) == timestamp(ohlc[1:8])
        @test values(dropnan(uohlc, :all))    == values(ohlc[1:8])

        @test values(dropnan(ohlc, :any))     == values(ohlc)
        @test timestamp(dropnan(nohlc, :any)) == timestamp(ohlc)[[1:6;13:end]]
        @test values(dropnan(nohlc, :any))    == values(ohlc)[[1:6;13:end], :]
        @test timestamp(dropnan(uohlc, :any)) == timestamp(ohlc[1:8])
        @test values(dropnan(uohlc, :any))    == values(ohlc[1:8])
    end
end

end  # @testset "apply"


text: @info("Add custom startup code HERE in ./startup.jl")

ext: jl

max_stars_repo_name: ma-laforge/SampleConventionalApp

max_issues_repo_name: ma-laforge/SampleConventionalApp

max_forks_repo_name: ma-laforge/SampleConventionalApp


text:
@testset "mexp1" begin
    @markov midplane(lam, mu) begin
        @tr begin
            :up => :down, lam
            :down => :up, mu
        end
        @init begin
            :up, 1
        end
    end

    begin
        lam = 0.4
        mu = 1.5
        delta = 1.5
    end

    Q, x, _, states = generate(midplane(lam, mu))
    println(tprob(1.0, Q, x))
end

@testset "mexp2" begin
    @markov midplane(lam, mu) begin
        @tr begin
            :up => :down, lam
            :down => :up, mu
        end
        @init begin
            :up, 1
        end
    end

    @bind begin
        lam = 0.4
        mu = 1.5
        delta = 1.5
    end

    Q, x, _, states = generate(midplane(lam, mu))
    println(seval(tprob(1.0, Q, x)))
end

@testset "mexp3" begin
    @markov midplane(lam, mu) begin
        @tr begin
            :up => :down, lam
            :down => :up, mu
        end
        @init begin
            :up, 1
        end
    end

    @bind begin
        lam = 0.4
        mu = 1.5
        delta = 1.5
    end

    Q, x, _, states = generate(midplane(lam, mu))
    println(seval(tprob(1.0, Q, x), :lam))
end

@testset "mexp3" begin
    @markov midplane(lam, mu, p) begin
        @tr begin
            :up => :down, lam
            :down => :up, mu
        end
        @init begin
            :up, p
            :down, 1-p
        end
    end

    @bind begin
        lam = 0.4
        mu = 1.5
        delta = 1.5
        p = 1.0
    end

    Q, x, _, states = generate(midplane(lam, mu, p))
    println(seval(tprob(1.0, Q, x), :p))
    println(seval(tprob(1.0, Q, x), :lam))
    println(seval(tprob(1.0, Q, x), (:lam, :p)))
    println(seval(ctprob(100.0, Q, x), :p))
    println(seval(ctprob(100.0, Q, x), :lam))
    println(seval(ctprob(100.0, Q, x), (:lam, :p)))
end


text:
using Sockets: TCPServer
using Flux: destructure, params

"""
  `fedavg!(models_dict::Dict{String, Any}, round_id::String)` -> Flux model

Averages models stored at `models_dict[round_id]`
Deletes those models from `models_dict`
Returns averaged model

Example:
```
# Let m1, m2, m3 be Flux models
models_dict = Dict("some_id" => [m1, m2, m3])
averaged_model = fedavg!(models_dict, "some_id")
```
"""
function fedavg!(models_dict::Dict{String, Any}, round_id::String)
  models = models_dict[round_id]
  local restructure

  # Obtain model weights
  Ws = []
  for model in models
    W, restructure = destructure(model)
    push!(Ws, W)
  end

  # Average weights and form new model
  average(X) = sum(X) / length(X)
  averaged_model = restructure(average(Ws))

  # Remove used models from `model_dict`
  delete!(models_dict, round_id)

  return averaged_model
end

"""
  `update_model(server_model, client_model, learning_rate = 0.1)`

Obtains weighted average of two models and returns resulting model
"""
function update_model(server_model, client_model, learning_rate = 0.1)
  server_W, restructure = destructure(server_model)
  client_W, restructure = destructure(client_model)

  server_degree = Float32(1 - learning_rate)
  client_degree = Float32(learning_rate)

  result_W = (server_W .* server_degree) .+ (client_W .* client_degree)
  return restructure(result_W) # resulting model
end

"""
  `establish_round_id(con::TCPSocket, shared_key::Bytes)`

Establishes Round ID for given client
"""
function establish_round_id(con::TCPSocket, shared_key::Bytes)
  # TODO
  return "not_implemented"
end

"""
  `serve_model(server::TCPServer, model, round_len=3, learning_rate=0.1, verbose=true)`

Sends Flux model to clients, asynchronously

`round_len`: number of clients per round
`learning_rate`: degree that averaged client models replace server model per round
`verbose`: if true, print parameter information
"""
function serve_model(server::TCPServer, model, round_len = 3, learning_rate = 0.1, verbose=true)
  num_connected = 0
  client_models = Dict{String, Any}()

  while true
    con = accept(server)
    if isopen(con)
      num_connected += 1
      println("Clients connected: $num_connected")
    end

    @async while isopen(con)
      try
        # Establish AES-256 Key
        shared_key = handshake_server(con)

        # Securely send model to client
        model_bytes = model_to_bytes(model)
        secure_send_bytes(con, model_bytes, shared_key)

        # Establish Round ID
        round_id = establish_round_id(con, shared_key)

        # Receive Updated Model
        client_update_bytes = secure_receive_bytes(con, shared_key)
        client_update = model_from_bytes(client_update_bytes)

        # Store Updated Model
        if round_id in keys(client_models)
          push!(client_models[round_id], client_update)
        else
          client_models[round_id] = [client_update]
        end

        # Close Connection
        num_connected -= 1
        close(con)
        println("Clients connected: $num_connected")

        # Check number of models stored for each `round_id`
        for round_id in keys(client_models)
          if verbose
            println("Checking Round ID: $round_id (length: $(length(client_models[round_id])))")
          end

          # Average models, if models received >= round_len
          if length(client_models[round_id]) >= round_len
            averaged_model = fedavg!(client_models, round_id)
            model = update_model(model, averaged_model, learning_rate)

            if verbose
              println("Averaging Round ID: $round_id")
              println("\nAveraged Client Params\n", params(averaged_model))
              println("\nUpdated Model Params\n", params(model))
            end
          end
        end

      catch err
        close(con)
        throw(err)
      end
    end
  end
end


text:
"""
    OptSummary

Summary of an `NLopt` optimization

# Fields

* `initial`: a copy of the initial parameter values in the optimization
* `lowerbd`: lower bounds on the parameter values
* `ftol_rel`: as in NLopt
* `ftol_abs`: as in NLopt
* `xtol_rel`: as in NLopt
* `xtol_abs`: as in NLopt
* `initial_step`: as in NLopt
* `maxfeval`: as in NLopt (`maxeval`)
* `maxtime`: as in NLopt
* `final`: a copy of the final parameter values from the optimization
* `fmin`: the final value of the objective
* `feval`: the number of function evaluations
* `optimizer`: the name of the optimizer used, as a `Symbol`
* `returnvalue`: the return value, as a `Symbol`
* `nAGQ`: number of adaptive Gauss-Hermite quadrature points in deviance evaluation for GLMMs
* `REML`: use the REML criterion for LMM fits
* `sigma`: a priori value for the residual standard deviation for LMM

The latter three fields are model characteristics and not related directly to the `NLopt` package or algorithms.
"""
mutable struct OptSummary{T<:AbstractFloat}
    initial::Vector{T}
    lowerbd::Vector{T}
    finitial::T
    ftol_rel::T
    ftol_abs::T
    xtol_rel::T
    xtol_abs::Vector{T}
    initial_step::Vector{T}
    maxfeval::Int
    maxtime::T
    final::Vector{T}
    fmin::T
    feval::Int
    optimizer::Symbol
    returnvalue::Symbol
    nAGQ::Integer           # don't really belong here but I needed a place to store them
    REML::Bool
    sigma::Union{T,Nothing}
end
function OptSummary(
    initial::Vector{T},
    lowerbd::Vector{T},
    optimizer::Symbol;
    ftol_rel::T=zero(T),
    ftol_abs::T=zero(T),
    xtol_rel::T=zero(T),
    initial_step::Vector{T}=T[],
    maxfeval=-1,
    maxtime=T(-1),
) where {T<:AbstractFloat}
    return OptSummary(
        initial,
        lowerbd,
        T(Inf),
        ftol_rel,
        ftol_abs,
        xtol_rel,
        zero(initial),
        initial_step,
        maxfeval,
        maxtime,
        copy(initial),
        T(Inf),
        -1,
        optimizer,
        :FAILURE,
        1,
        false,
        nothing,
    )
end

function Base.show(io::IO, ::MIME"text/plain", s::OptSummary)
    println(io, "Initial parameter vector: ", s.initial)
    println(io, "Initial objective value:  ", s.finitial)
    println(io)
    println(io, "Optimizer (from NLopt):   ", s.optimizer)
    println(io, "Lower bounds:             ", s.lowerbd)
    println(io, "ftol_rel:                 ", s.ftol_rel)
    println(io, "ftol_abs:                 ", s.ftol_abs)
    println(io, "xtol_rel:                 ", s.xtol_rel)
    println(io, "xtol_abs:                 ", s.xtol_abs)
    println(io, "initial_step:             ", s.initial_step)
    println(io, "maxfeval:                 ", s.maxfeval)
    println(io, "maxtime:                  ", s.maxtime)
    println(io)
    println(io, "Function evaluations:     ", s.feval)
    println(io, "Final parameter vector:   ", s.final)
    println(io, "Final objective value:    ", s.fmin)
    return println(io, "Return code:              ", s.returnvalue)
end

Base.show(io::IO, s::OptSummary) = Base.show(io, MIME"text/plain"(), s)

function NLopt.Opt(optsum::OptSummary)
    lb = optsum.lowerbd

    opt = NLopt.Opt(optsum.optimizer, length(lb))
    NLopt.ftol_rel!(opt, optsum.ftol_rel) # relative criterion on objective
    NLopt.ftol_abs!(opt, optsum.ftol_abs) # absolute criterion on objective
    NLopt.xtol_rel!(opt, optsum.xtol_rel) # relative criterion on parameter values
    if length(optsum.xtol_abs) == length(lb)  # not true for fast=false optimization in GLMM
        NLopt.xtol_abs!(opt, optsum.xtol_abs) # absolute criterion on parameter values
    end
    NLopt.lower_bounds!(opt, lb)
    NLopt.maxeval!(opt, optsum.maxfeval)
    NLopt.maxtime!(opt, optsum.maxtime)
    if isempty(optsum.initial_step)
        optsum.initial_step = NLopt.initial_step(opt, optsum.initial, similar(lb))
    else
        NLopt.initial_step!(opt, optsum.initial_step)
    end
    return opt
end

StructTypes.StructType(::Type{<:OptSummary}) = StructTypes.Mutable()
StructTypes.excludes(::Type{<:OptSummary}) = (:lowerbd,)

const _NLOPT_FAILURE_MODES = [
    :FAILURE,
    :INVALID_ARGS,
    :OUT_OF_MEMORY,
    :FORCED_STOP,
    :MAXEVAL_REACHED,
    :MAXTIME_REACHED,
]

function _check_nlopt_return(ret)
    ret == :ROUNDOFF_LIMITED && @warn("NLopt was roundoff limited")
    if ret ∈ _NLOPT_FAILURE_MODES
        @warn("NLopt optimization failure: $ret")
    end
end


text:
using .TurbulenceClosures: ▶z_aaf

using Oceananigans.Operators

const Tx = 16 # CUDA threads per x-block
const Ty = 16 # CUDA threads per y-block

"""
    time_step!(model, Nt, Δt; init_with_euler=true)

Step forward `model` `Nt` time steps with step size `Δt`.

If `init_with_euler` is set to true, then the first step will be taken using a first-order
forward Euler method.
"""
function time_step!(model, Nt, Δt; init_with_euler=true)

    if model.clock.iteration == 0
        [ run_diagnostic(model, diag) for diag in values(model.diagnostics) ]
        [ write_output(model, out)    for out  in values(model.output_writers) ]
    end

    FT = eltype(model.grid)
    U, C, Gⁿ, G⁻, K, P = datatuples(model.velocities, model.tracers, model.timestepper.Gⁿ,
                                    model.timestepper.G⁻, model.diffusivities, model.pressures)

    for n in 1:Nt
        χ = ifelse(init_with_euler && n==1, FT(-0.5), model.timestepper.χ)

        adams_bashforth_time_step!(model, U, C, P, K, Gⁿ, G⁻, Δt, χ)

        [ time_to_run(model.clock, diag) && run_diagnostic(model, diag) for diag in values(model.diagnostics) ]
        [ time_to_run(model.clock, out) && write_output(model, out) for out in values(model.output_writers) ]
    end

    return nothing
end

time_step!(model; Nt, Δt, kwargs...) = time_step!(model, Nt, Δt; kwargs...)

"""
Step forward one time step with a 2nd-order Adams-Bashforth method and pressure-correction
substep.
"""
function adams_bashforth_time_step!(model, U, C, P, K, Gⁿ, G⁻, Δt, χ)

    arch = model.architecture
    grid = model.grid
    bc_args = (model.clock.time, model.clock.iteration, U, C, model.parameters)

    # Pre-computations:
    store_previous_source_terms!(G⁻, arch, grid, Gⁿ)
    fill_halo_regions!(merge(U, C), model.boundary_conditions.solution, arch, grid, bc_args...)

    calculate_diffusivities!(K, arch, grid, model.closure, model.buoyancy, U, C)
    fill_halo_regions!(K, model.boundary_conditions.pressure, arch, grid) # diffusivities share bcs with pressure.

    @launch device(arch) config=launch_config(grid, 2) update_hydrostatic_pressure!(P.pHY′, grid, model.buoyancy, C)
    fill_halo_regions!(P.pHY′, model.boundary_conditions.pressure, arch, grid)

    # Calculate tendency terms (minus non-hydrostatic pressure, which is updated in a pressure correction step):
    calculate_interior_source_terms!(Gⁿ, arch, grid, model.coriolis, model.closure, U, C, P.pHY′, K, model.forcing,
                                     model.parameters, model.clock.time)
    calculate_boundary_source_terms!(Gⁿ, model.boundary_conditions.solution, arch, grid, bc_args...)

    # Complete explicit substep:
    adams_bashforth_update_source_terms!(Gⁿ, arch, grid, χ, G⁻)

    # Start pressure correction substep with a pressure solve:
    fill_halo_regions!(Gⁿ[1:3], model.boundary_conditions.tendency[1:3], arch, grid)
    @launch device(arch) config=launch_config(grid, 3) calculate_poisson_right_hand_side!(model.poisson_solver.storage, arch, grid,
                                                                                model.poisson_solver.bcs, Δt, U, Gⁿ)
    solve_for_pressure!(P.pNHS, arch, grid, model.poisson_solver, model.poisson_solver.storage)
    fill_halo_regions!(P.pNHS, model.boundary_conditions.pressure, arch, grid)

    # Complete pressure correction step:
    update_solution!(U, C, arch, grid, Gⁿ, P.pNHS, Δt)

    # Recompute vertical velocity w from continuity equation to ensure incompressibility
    fill_halo_regions!(U, model.boundary_conditions.solution[1:3], arch, grid, bc_args...)
    @launch device(arch) config=launch_config(grid, 2) compute_w_from_continuity!(U, grid)

    model.clock.time += Δt
    model.clock.iteration += 1

    return nothing
end

#####
##### Navier-Stokes and tracer advection equations
#####

""" Calculate the right-hand-side of the u-momentum equation. """
function calculate_Gu!(Gu, grid, coriolis, closure, U, C, K, F, pHY′, parameters, time)
    @loop for k in (1:grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                @inbounds Gu[i, j, k] = (-u∇u(grid, U.u, U.v, U.w, i, j, k)
                                            - x_f_cross_U(i, j, k, grid, coriolis, U)
                                            - ∂x_p(i, j, k, grid, pHY′)
                                            + ∂ⱼ_2ν_Σ₁ⱼ(i, j, k, grid, closure, U.u, U.v, U.w, K)
                                            + F.u(i, j, k, grid, time, U, C, parameters))
            end
        end
    end
    return nothing
end

""" Calculate the right-hand-side of the v-momentum equation. """
function calculate_Gv!(Gv, grid, coriolis, closure, U, C, K, F, pHY′, parameters, time)
    @loop for k in (1:grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                @inbounds Gv[i, j, k] = (-u∇v(grid, U.u, U.v, U.w, i, j, k)
                                            - y_f_cross_U(i, j, k, grid, coriolis, U)
                                            - ∂y_p(i, j, k, grid, pHY′)
                                            + ∂ⱼ_2ν_Σ₂ⱼ(i, j, k, grid, closure, U.u, U.v, U.w, K)
                                            + F.v(i, j, k, grid, time, U, C, parameters))
            end
        end
    end
    return nothing
end

""" Calculate the right-hand-side of the w-momentum equation. """
function calculate_Gw!(Gw, grid, coriolis, closure, U, C, K, F, parameters, time)
    @loop for k in (1:grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                @inbounds Gw[i, j, k] = (-u∇w(grid, U.u, U.v, U.w, i, j, k)
                                            - z_f_cross_U(i, j, k, grid, coriolis, U)
                                            + ∂ⱼ_2ν_Σ₃ⱼ(i, j, k, grid, closure, U.u, U.v, U.w, K)
                                            + F.w(i, j, k, grid, time, U, C, parameters))
            end
        end
    end
    return nothing
end

""" Calculate the right-hand-side of the tracer advection-diffusion equation. """
function calculate_Gc!(Gc, grid, closure, c, tracer_index, U, C, K, Fc, parameters, time)
    @loop for k in (1:grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                @inbounds Gc[i, j, k] = (-div_flux(grid, U.u, U.v, U.w, c, i, j, k)
                                            + ∇_κ_∇c(i, j, k, grid, c, tracer_index, closure, K)
                                            + Fc(i, j, k, grid, time, U, C, parameters))
            end
        end
    end
    return nothing
end

""" Store previous value of the source term and calculate current source term. """
function calculate_interior_source_terms!(G, arch, grid, coriolis, closure, U, C, pHY′, K, F, parameters, time)

    Bx, By, Bz = floor(Int, grid.Nx/Tx), floor(Int, grid.Ny/Ty), grid.Nz  # Blocks in grid

    @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_Gu!(G.u, grid, coriolis, closure, U, C, K, F,
                                                                            pHY′, parameters, time)

    @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_Gv!(G.v, grid, coriolis, closure, U, C, K, F,
                                                                            pHY′, parameters, time)

    @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_Gw!(G.w, grid, coriolis, closure, U, C, K, F,
                                                                            parameters, time)

    for tracer_index in 1:length(C)
        @inbounds Gc = G[tracer_index+3]
        @inbounds Fc = F[tracer_index+3]
        @inbounds  c = C[tracer_index]
        @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_Gc!(Gc, grid, closure, c, Val(tracer_index),
                                                                                U, C, K, Fc, parameters, time)
    end

    return nothing
end

""" Apply boundary conditions by adding flux divergences to the right-hand-side. """
function calculate_boundary_source_terms!(Gⁿ, bcs, arch, grid, args...)

    # Velocity fields
    for (i, ubcs) in enumerate(bcs[1:3])
        apply_z_bcs!(Gⁿ[i], arch, grid, ubcs.z.bottom, ubcs.z.top, args...)
    end

    # Tracer fields
    for (i, cbcs) in enumerate(bcs[4:end])
        apply_z_bcs!(Gⁿ[i+3], arch, grid, cbcs.z.bottom, cbcs.z.top, args...)
    end

    return nothing
end

#####
##### Pressure-related functions
#####

"Solve the Poisson equation for non-hydrostatic pressure on the CPU."
function solve_for_pressure!(pressure, ::CPU, grid, poisson_solver, ϕ)
    solve_poisson_3d!(poisson_solver, grid)
    view(pressure, 1:grid.Nx, 1:grid.Ny, 1:grid.Nz) .= real.(ϕ)
    return nothing
end

"Solve the Poisson equation for non-hydrostatic pressure on the GPU."
function solve_for_pressure!(pressure, ::GPU, grid, poisson_solver, ϕ)
    solve_poisson_3d!(poisson_solver, grid)
    @launch device(GPU()) config=launch_config(grid, 3) idct_permute!(pressure, grid, poisson_solver.bcs, ϕ)
    return nothing
end

"""
Update the hydrostatic pressure perturbation pHY′. This is done by integrating
the `buoyancy_perturbation` downwards:

    `pHY′ = ∫ buoyancy_perturbation dz` from `z=0` down to `z=-Lz`
"""
function update_hydrostatic_pressure!(pHY′, grid, buoyancy, C)
    @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
        @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
            @inbounds pHY′[i, j, grid.Nz] = - ▶z_aaf(i, j, grid.Nz, grid, buoyancy_perturbation, buoyancy, C) * grid.Δz
            @unroll for k in grid.Nz-1 : -1 : 1
                @inbounds pHY′[i, j, k] =
                    pHY′[i, j, k+1] - ▶z_aaf(i, j, k+1, grid, buoyancy_perturbation, buoyancy, C) * grid.Δz
            end
        end
    end
    return nothing
end

"""
Calculate the right-hand-side of the Poisson equation for the non-hydrostatic
pressure

    `∇²ϕ_{NH}^{n+1} = (∇·u^n)/Δt + ∇·(Gu, Gv, Gw)`
"""
function calculate_poisson_right_hand_side!(RHS, ::CPU, grid, ::PoissonBCs, Δt, U, G)
    @loop for k in (1:grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                # Calculate divergence of the RHS source terms (Gu, Gv, Gw).
                @inbounds RHS[i, j, k] = div_f2c(grid, U.u, U.v, U.w, i, j, k) / Δt +
                                         div_f2c(grid, G.u, G.v, G.w, i, j, k)
            end
        end
    end

    return nothing
end

"""
Calculate the right-hand-side of the Poisson equation for the non-hydrostatic
pressure and in the process apply the permutation

    [a, b, c, d, e, f, g, h] -> [a, c, e, g, h, f, d, b]

in the z-direction which is required by the GPU fast cosine transform algorithm for
horizontally periodic model configurations.
"""
function calculate_poisson_right_hand_side!(RHS, ::GPU, grid, ::PPN, Δt, U, G)
    Nx, Ny, Nz = grid.Nx, grid.Ny, grid.Nz
    @loop for k in (1:Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
        @loop for j in (1:Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                if (k & 1) == 1  # isodd(k)
                    k′ = convert(UInt32, CUDAnative.floor(k/2) + 1)
                else
                    k′ = convert(UInt32, Nz - CUDAnative.floor((k-1)/2))
                end
                @inbounds RHS[i, j, k′] = div_f2c(grid, U.u, U.v, U.w, i, j, k) / Δt +
                                          div_f2c(grid, G.u, G.v, G.w, i, j, k)
            end
        end
    end

    return nothing
end

"""
Calculate the right-hand-side of the Poisson equation for the non-hydrostatic
pressure and in the process apply the permutation

    [a, b, c, d, e, f, g, h] -> [a, c, e, g, h, f, d, b]

in the y- and z-directions which is required by the GPU fast cosine transform algorithm for
reentrant channel model configurations.
"""
function calculate_poisson_right_hand_side!(RHS, ::GPU, grid, ::PNN, Δt, U, G)
    Nx, Ny, Nz = grid.Nx, grid.Ny, grid.Nz
    @loop for k in (1:grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
        @loop for j in (1:Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                if (k & 1) == 1  # isodd(k)
                    k′ = convert(UInt32, CUDAnative.floor(k/2) + 1)
                else
                    k′ = convert(UInt32, Nz - CUDAnative.floor((k-1)/2))
                end

                if (j & 1) == 1  # isodd(j)
                    j′ = convert(UInt32, CUDAnative.floor(j/2) + 1)
                else
                    j′ = convert(UInt32, Ny - CUDAnative.floor((j-1)/2))
                end

                @inbounds RHS[i, j′, k′] = div_f2c(grid, U.u, U.v, U.w, i, j, k) / Δt +
                                           div_f2c(grid, G.u, G.v, G.w, i, j, k)
            end
        end
    end

    return nothing
end

"""
Copy the non-hydrostatic pressure into `pNHS` and undo the permutation

    [a, b, c, d, e, f, g, h] -> [a, c, e, g, h, f, d, b]

along the z-direction.
"""
function idct_permute!(pNHS, grid, ::PPN, ϕ)
    Nx, Ny, Nz = grid.Nx, grid.Ny, grid.Nz
    @loop for k in (1:Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
        @loop for j in (1:Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                if k <= Nz/2
                    @inbounds pNHS[i, j, 2k-1] = real(ϕ[i, j, k])
                else
                    @inbounds pNHS[i, j, 2(Nz-k+1)] = real(ϕ[i, j, k])
                end
            end
        end
    end

    return nothing
end

"""
Copy the non-hydrostatic pressure into `pNHS` and undo the permutation

    [a, b, c, d, e, f, g, h] -> [a, c, e, g, h, f, d, b]

along the y- and z-direction.
"""
function idct_permute!(pNHS, grid, ::PNN, ϕ)
    Nx, Ny, Nz = grid.Nx, grid.Ny, grid.Nz
    @loop for k in (1:grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
        @loop for j in (1:Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                if k <= Nz/2
                    k′ = 2k-1
                else
                    k′ = 2(Nz-k+1)
                end

                if j <= Ny/2
                    j′ = 2j-1
                else
                    j′ = 2(Ny-j+1)
                end

                @inbounds pNHS[i, j′, k′] = real(ϕ[i, j, k])
            end
        end
    end

    return nothing
end

#####
##### Adams-Bashforth stuff
#####

""" Store previous source terms for `u`, `v`, and `w` before updating them. """
function store_previous_velocity_source_terms!(G⁻, grid, Gⁿ)
    @loop for k in (1:grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                @inbounds G⁻.u[i, j, k] = Gⁿ.u[i, j, k]
                @inbounds G⁻.v[i, j, k] = Gⁿ.v[i, j, k]
                @inbounds G⁻.w[i, j, k] = Gⁿ.w[i, j, k]
            end
        end
    end
    return nothing
end

""" Store previous source terms for a tracer before updating them. """
function store_previous_tracer_source_term!(Gc⁻, grid, Gcⁿ)
    @loop for k in (1:grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                @inbounds Gc⁻[i, j, k] = Gcⁿ[i, j, k]
            end
        end
    end
    return nothing
end

""" Store previous source terms before updating them. """
function store_previous_source_terms!(G⁻, arch, grid, Gⁿ)

    # Velocity fields
    @launch device(arch) config=launch_config(grid, 3) store_previous_velocity_source_terms!(G⁻, grid, Gⁿ)

    # Tracer fields
    for i in 4:length(G⁻)
        @inbounds Gc⁻ = G⁻[i]
        @inbounds Gcⁿ = Gⁿ[i]
        @launch device(arch) config=launch_config(grid, 3) store_previous_tracer_source_term!(Gc⁻, grid, Gcⁿ)
    end

    return nothing
end

"""
Evaluate the right-hand-side terms at time step n+½ using a weighted 2nd-order
Adams-Bashforth method
    `G^{n+½} = (3/2 + χ)G^{n} - (1/2 + χ)G^{n-1}`
"""
function adams_bashforth_update_velocity_source_terms!(Gⁿ, grid::AbstractGrid{FT}, χ, G⁻) where FT
    @loop for k in (1:grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                @inbounds Gⁿ.u[i, j, k] = (FT(1.5) + χ) * Gⁿ.u[i, j, k] - (FT(0.5) + χ) * G⁻.u[i, j, k]
                @inbounds Gⁿ.v[i, j, k] = (FT(1.5) + χ) * Gⁿ.v[i, j, k] - (FT(0.5) + χ) * G⁻.v[i, j, k]
                @inbounds Gⁿ.w[i, j, k] = (FT(1.5) + χ) * Gⁿ.w[i, j, k] - (FT(0.5) + χ) * G⁻.w[i, j, k]
            end
        end
    end

    return nothing
end

"""
Evaluate the right-hand-side terms at time step n+½ using a weighted 2nd-order
Adams-Bashforth method
    `G^{n+½} = (3/2 + χ)G^{n} - (1/2 + χ)G^{n-1}`
"""
function adams_bashforth_update_tracer_source_term!(Gcⁿ, grid::AbstractGrid{FT}, χ, Gc⁻) where FT
    @loop for k in (1:grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                @inbounds Gcⁿ[i, j, k] = (FT(1.5) + χ) * Gcⁿ[i, j, k] - (FT(0.5) + χ) * Gc⁻[i, j, k]
            end
        end
    end

    return nothing
end

"""
Evaluate the right-hand-side terms for velocity fields and tracer fields
at time step n+½ using a weighted 2nd-order Adams-Bashforth method.
"""
function adams_bashforth_update_source_terms!(Gⁿ, arch, grid, χ, G⁻)
    # Velocity fields
    @launch device(arch) config=launch_config(grid, 3) adams_bashforth_update_velocity_source_terms!(Gⁿ, grid, χ, G⁻)

    # Tracer fields
    for i in 4:length(Gⁿ)
        @inbounds Gcⁿ = Gⁿ[i]
        @inbounds Gc⁻ = G⁻[i]
        @launch device(arch) config=launch_config(grid, 3) adams_bashforth_update_tracer_source_term!(Gcⁿ, grid, χ, Gc⁻)
    end

    return nothing
end

"""
Update the horizontal velocities u and v via
    `u^{n+1} = u^n + (Gu^{n+½} - δₓp_{NH} / Δx) Δt`
Note that the vertical velocity is not explicitly time stepped.
"""
function update_velocities!(U, grid, Δt, G, pNHS)
    @loop for k in (1:grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                @inbounds U.u[i, j, k] += (G.u[i, j, k] - ∂x_p(i, j, k, grid, pNHS)) * Δt
                @inbounds U.v[i, j, k] += (G.v[i, j, k] - ∂y_p(i, j, k, grid, pNHS)) * Δt
            end
        end
    end

    return nothing
end

"""
Update tracers via
    `c^{n+1} = c^n + Gc^{n+½} Δt`
"""
function update_tracer!(c, grid, Δt, Gc)
    @loop for k in (1:grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                @inbounds c[i, j, k] += Gc[i, j, k] * Δt
            end
        end
    end

    return nothing
end

"Update the solution variables (velocities and tracers)."
function update_solution!(U, C, arch, grid, G, pNHS, Δt)
    @launch device(arch) config=launch_config(grid, 3) update_velocities!(U, grid, Δt, G, pNHS)

    for i in 1:length(C)
        @inbounds c = C[i]
        @inbounds Gc = G[i+3]
        @launch device(arch) config=launch_config(grid, 3) update_tracer!(c, grid, Δt, Gc)
    end

    return nothing
end

"""
Compute the vertical velocity w by integrating the continuity equation from the bottom upwards

    `w^{n+1} = -∫ [∂/∂x (u^{n+1}) + ∂/∂y (v^{n+1})] dz`
"""
function compute_w_from_continuity!(U, grid)
    @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
        @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
            # U.w[i, j, 0] = 0 is enforced via halo regions.
            @unroll for k in 2:grid.Nz
                @inbounds U.w[i, j, k] = U.w[i, j, k-1] - grid.Δz * ∇h_u(i, j, k-1, grid, U.u, U.v)
            end
        end
    end

    return nothing
end


text:
# This file was generated, do not modify it. # hide
using MLJ
MLJ.color_off() # hide
using HTTP
using CSV
import DataFrames: DataFrame, select!, Not
req1 = HTTP.get("http://archive.ics.uci.edu/ml/machine-learning-databases/horse-colic/horse-colic.data")
req2 = HTTP.get("http://archive.ics.uci.edu/ml/machine-learning-databases/horse-colic/horse-colic.test")
header = ["surgery", "age", "hospital_number",
    "rectal_temperature", "pulse",
    "respiratory_rate", "temperature_extremities",
    "peripheral_pulse", "mucous_membranes",
    "capillary_refill_time", "pain",
    "peristalsis", "abdominal_distension",
    "nasogastric_tube", "nasogastric_reflux",
    "nasogastric_reflux_ph", "feces", "abdomen",
    "packed_cell_volume", "total_protein",
    "abdomcentesis_appearance", "abdomcentesis_total_protein",
    "outcome", "surgical_lesion", "lesion_1", "lesion_2", "lesion_3",
    "cp_data"]
csv_opts = (header=header, delim=' ', missingstring="?",
            ignorerepeated=true)
data_train = CSV.read(req1.body; csv_opts...)
data_test  = CSV.read(req2.body; csv_opts...)
@show size(data_train)
@show size(data_test)


text:
module FrostHelperCustomSpring

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/SpringFloor" CustomSpring(x::Integer, y::Integer, playerCanUse::Bool=true, directory::String="objects/spring/", speedMult::Number=1.0, oneUse::Bool=false, renderOutline::Bool=true)
@mapdef Entity "FrostHelper/SpringRight" CustomSpringRight(x::Integer, y::Integer, playerCanUse::Bool=true, directory::String="objects/spring/", speedMult::String="1.0,1.0", oneUse::Bool=false, renderOutline::Bool=true)
@mapdef Entity "FrostHelper/SpringLeft" CustomSpringLeft(x::Integer, y::Integer, playerCanUse::Bool=true, directory::String="objects/spring/", speedMult::String="1.0,1.0", oneUse::Bool=false, renderOutline::Bool=true)
@mapdef Entity "FrostHelper/SpringCeiling" CustomSpringCeiling(x::Integer, y::Integer, playerCanUse::Bool=true, directory::String="objects/spring/", speedMult::Number=1.0, oneUse::Bool=false, renderOutline::Bool=true)

const placements = Ahorn.PlacementDict(
    "Custom Spring (Up, Frost Helper)" => Ahorn.EntityPlacement(
        CustomSpring
    ),
    "Custom Spring (Left, Frost Helper)" => Ahorn.EntityPlacement(
        CustomSpringRight
    ),
    "Custom Spring (Right, Frost Helper)" => Ahorn.EntityPlacement(
        CustomSpringLeft
    ),
    "Custom Spring (Down, Frost Helper)" => Ahorn.EntityPlacement(
        CustomSpringCeiling
    ),
)

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

    return Ahorn.Rectangle(x - 6, y - 3, 12, 5)
end

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

    entity.data["speedMult"] = string(get(entity.data, "speedMult", "1.0"))
    return Ahorn.Rectangle(x - 1, y - 6, 5, 12)
end

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

    return Ahorn.Rectangle(x - 4, y - 6, 5, 12)
end

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

    return Ahorn.Rectangle(x - 6, y - 1, 12, 5)
end

sprite = "objects/spring/00.png"

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CustomSpring, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, -8)
function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CustomSpringLeft, room::Maple.Room)
    if get(entity.data, "speedMult", "1.0") isa Number
        entity.data["speedMult"] = string(get(entity.data, "speedMult", "1.0"))
    end
    Ahorn.drawSprite(ctx, sprite, 9, -11, rot=pi / 2)
end
function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CustomSpringRight, room::Maple.Room)
    if get(entity.data, "speedMult", "1.0") isa Number
        entity.data["speedMult"] = string(get(entity.data, "speedMult", "1.0"))
    end

    Ahorn.drawSprite(ctx, sprite, 3, 1, rot=-pi / 2)
end
Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CustomSpringCeiling, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 12, -2, rot=pi)

function getAllAttributes(entity)
    return entity.x, entity.y, Bool(get(entity.data, "playerCanUse", true)), get(entity.data, "directory", "objects/spring/"), string(get(entity.data, "speedMult", 1.0)), Bool(get(entity.data, "oneUse", false)), Bool(get(entity.data, "renderOutline", true))
end

function turnIntoRight(entity)
    x, y, playerCanUse, directory, speedMult, oneUse, renderOutline = getAllAttributes(entity)
    return CustomSpringRight(x, y, playerCanUse, directory, speedMult, oneUse, renderOutline)
end

function turnIntoLeft(entity)
    x, y, playerCanUse, directory, speedMult, oneUse, renderOutline = getAllAttributes(entity)
    return CustomSpringLeft(x, y, playerCanUse, directory, speedMult, oneUse, renderOutline)
end

function turnIntoUp(entity)
    x, y, playerCanUse, directory, speedMult, oneUse, renderOutline = getAllAttributes(entity)
    return CustomSpring(x, y, playerCanUse, directory, parse(Float32, speedMult), oneUse, renderOutline)
end

function turnIntoDown(entity)
    x, y, playerCanUse, directory, speedMult, oneUse, renderOutline = getAllAttributes(entity)
    return CustomSpringCeiling(x, y, playerCanUse, directory, parse(Float32, speedMult), oneUse, renderOutline)
end

function Ahorn.flipped(entity::CustomSpringLeft, horizontal::Bool)
    if horizontal
        return turnIntoRight(entity)
    end

    return turnIntoUp(entity)
end

function Ahorn.flipped(entity::CustomSpringRight, horizontal::Bool)
    if horizontal
        return turnIntoLeft(entity)
    end

    return turnIntoUp(entity)
end

function Ahorn.flipped(entity::CustomSpring, horizontal::Bool)
    if horizontal
        return turnIntoRight(entity)
    end

    return turnIntoDown(entity)
end

function Ahorn.flipped(entity::CustomSpringCeiling, horizontal::Bool)
    if horizontal
        return turnIntoRight(entity)
    end

    return turnIntoUp(entity)
end

function Ahorn.rotated(entity::CustomSpringRight, steps::Int)
    if steps > 0 # r
        return Ahorn.rotated(turnIntoUp(entity), steps - 1)

    elseif steps < 0 # l
        return Ahorn.rotated(turnIntoDown(entity), steps + 1)
    end

    return entity
end

function Ahorn.rotated(entity::CustomSpringLeft, steps::Int)
    if steps > 0 # r
        return Ahorn.rotated(turnIntoDown(entity), steps - 1)

    elseif steps < 0 # l
        return Ahorn.rotated(turnIntoUp(entity), steps + 1)
    end

    return entity
end

function Ahorn.rotated(entity::CustomSpring, steps::Int)
    if steps > 0 # r
        return Ahorn.rotated(turnIntoLeft(entity), steps - 1)

    elseif steps < 0
        return Ahorn.rotated(turnIntoRight(entity), steps + 1)
    end

    return entity
end

function Ahorn.rotated(entity::CustomSpringCeiling, steps::Int)
    if steps > 0 # r
        return Ahorn.rotated(turnIntoRight(entity), steps - 1)

    elseif steps < 0
        return Ahorn.rotated(turnIntoLeft(entity), steps + 1)
    end

    return entity
end

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
∇ = differentiate(f, [x;y])
H = differentiate(∇, [x;y])

g = ∇ ⋅ ∇
v = [-∇[2]; ∇[1]]
h = v' * H * v
dg = differentiate(g, [x;y])
dh = differentiate(h, [x;y])

∇σ = g .* dh - ((3/2) * h).* dg

F₂ = [v ⋅ ∇σ; f]

curv_result = solve(F₂, start_system = :polyhedral)

max_stars_repo_name: Murgio/high-dimensional

max_issues_repo_name: Murgio/high-dimensional

max_forks_repo_name: Murgio/high-dimensional


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

using Test
using OpticSim

using Unitful
using Base.Filesystem
using DataFrames
using StaticArrays

@testset "GlassCat" begin
    # ensure there aren't any ambiguities
    @test isempty(detect_ambiguities(GlassCat))
    @test isempty(detect_unbound_args(GlassCat))

    @testset "Build Tests" begin
        # check that all automatic downloads are working
        for catname in split("HOYA NIKON OHARA SCHOTT Sumita")
            agffile = joinpath(OpticSim.GlassCat.AGF_DIR, catname * ".agf")
            @test isfile(agffile)
        end

        # check that particularly problematic glasses are parsing correctly
        @test !isnan(OpticSim.GlassCat.NIKON.LLF6.C10)
    end

    @testset "generate.jl" begin
        CATALOG_NAME = "TEST_CAT"
        SOURCE_DIR = joinpath(@__DIR__, "..", "..", "test")
        SOURCE_FILE = joinpath(SOURCE_DIR, "$(CATALOG_NAME).agf")

        TMP_DIR = mktempdir()
        MAIN_FILE = joinpath(TMP_DIR, "AGF_TEST_CAT.jl")

        TEST_CAT_VALUES = DataFrame(
            name = split("MG463_ MG523 _5MG448 MG452_STAR"),
            raw_name = split("MG463& MG523 5MG448 MG452*"),
            dispform = repeat([2], 4),
            Nd = repeat([1.0], 4),
            Vd = repeat([0.0], 4),
            exclude_sub = [5, 0, 0, 0],
            status = [6, 0, 0, 0],
            meltfreq = repeat([-1], 4),
            TCE = [19.8, 18.3, 30.2, 25.7],
            p = [4.63, 5.23, 4.48, 4.52],
            ΔPgF = repeat([0.0], 4),
            ignore_thermal_exp = [3, 0, 0, 0],
            C1 = [6.781409960E+00, 4.401393677E+00, 8.427948454E+00, 4.656568861E+00],
            C2 = [9.868538020E-02, 2.234722960E-02, 6.325268390E-02, 1.170534894E-01],
            C3 = [3.469529960E-03, 4.562935070E+00, -2.441576539E+00, 1.409322578E+00],
            C4 = [2.277823700E+00, 3.239318260E-01, 1.997453950E-03, 2.445302672E-04],
            C5 = [2.126055280E+00, 6.863473749E-01, 4.975938104E-01, 6.352738622E-01],
            C6 = [3.227490100E+03, 1.625888344E+03, 1.443307391E+03, 1.397305161E+03],
            C7 = repeat([0.0], 4),
            C8 = repeat([0.0], 4),
            C9 = repeat([0.0], 4),
            C10 = repeat([0.0], 4),
            D₀ = [2.226000000E-05, 1.122401430E-04, -1.718689957E-05, -2.545621409E-07],
            D₁ = [5.150100000E-08, 2.868995096E-13, -6.499841076E-13, 1.258087492E-10],
            D₂ = [1.309700000E-10, -4.157332734E-16, -4.153311088E-16, 1.119703940E-11],
            E₀ = [4.382400000E-05, 9.532969159E-05, 8.220254442E-06, 2.184273853E-05],
            E₁ = [4.966000000E-07, -1.045247704E-11, 6.086168578E-12, -2.788674395E-09],
            λₜₖ = [3.728000000E-01, 6.338475915E-01, 3.610058073E-01, -1.948154058E+00],
            temp = [2.002000000E+01, 2.000000000E+01, 2.000000000E+01, 2.000000000E+01],
            relcost = [1.1, -1.0, -1.0, -1.0],
            CR = [1.2, -1.0, -1.0, -1.0],
            FR = [1.3, -1.0, -1.0, -1.0],
            SR = [1.4, -1.0, -1.0, -1.0],
            AR = [1.5, -1.0, -1.0, -1.0],
            PR = [1.6, -1.0, -1.0, -1.0],
            λmin = repeat([2.0], 4),
            λmax = repeat([14.0], 4),
            transmission = [
                [[2.0, 1.0, 3.0], [14.0, 4.0, 5.0]],
                [[2.0, 1.0, 1.0], [14.0, 1.0, 1.0]],
                [[2.0, 1.0, 1.0], [14.0, 1.0, 1.0]],
                [[2.0, 1.0, 1.0], [14.0, 1.0, 1.0]],
            ],
        )
        FIELDS = names(TEST_CAT_VALUES)[2:end]

        @testset "Parsing Tests" begin
            cat = OpticSim.GlassCat.sourcefile_to_catalog(SOURCE_FILE)

            for glass in eachrow(TEST_CAT_VALUES)
                name = glass["name"]
                @test name ∈ keys(cat)
                for field in FIELDS
                    @test glass[field] == cat[name][field]
                end
            end
        end

        @testset "Module Gen Tests" begin
            OpticSim.GlassCat.generate_jls([CATALOG_NAME], MAIN_FILE, TMP_DIR, SOURCE_DIR, test=true)
            include(MAIN_FILE)

            for row in eachrow(TEST_CAT_VALUES)
                name = row["name"]
                @test "$CATALOG_NAME.$name" ∈ TEST_GLASS_NAMES

                glass = getfield(getfield(Main, Symbol(CATALOG_NAME)), Symbol(name))
                @test glass ∈ TEST_GLASSES

                for field in FIELDS
                    if field === "raw_name"
                    elseif field === "transmission"
                        @test row[field] == getfield(glass, Symbol(field))[1:length(row[field])]
                        @test zero(SVector{100-length(row[field]),SVector{3,Float64}}) == getfield(glass, Symbol(field))[length(row[field])+1:end]
                    else
                        @test row[field] == getfield(glass, Symbol(field))
                    end
                end
            end

            # these used to be in the "Glass Tests" testset, but they rely on the generated AGF_TEST_CAT.jl file
            g = TEST_CAT.MG523
            @test OpticSim.GlassCat.index(g, ((g.λmin + g.λmax) / 2)u"μm") ≈ 3.1560980389455593 atol=1e-14
            @test_throws ErrorException GlassCat.index(g, (g.λmin - 1)u"μm")
            @test_throws ErrorException GlassCat.index(g, (g.λmax + 1)u"μm")
        end
    end

    @testset "Glass Tests" begin
        @test isapprox(OpticSim.GlassCat.absairindex(500 * u"nm"), 1.0002741948670688, atol = 1e-14)
        @test isapprox(OpticSim.GlassCat.absairindex(600 * u"nm", temperature = 35 * u"°C", pressure = 2.0), 1.0005179096900811, atol = 1e-14)
        @test OpticSim.GlassCat.index(OpticSim.GlassCat.Air, 500 * u"nm") == 1.0
        @test OpticSim.GlassCat.index(OpticSim.GlassCat.Air, 600 * u"nm") == 1.0
        # @test isapprox(index(OpticSim.GlassCat.Vacuum, 500 * u"nm"), 0.9997258841958212, atol = 1e-14)

        # test against true values
        g = OpticSim.GlassCat.SCHOTT.N_BK7
        @test isapprox(OpticSim.GlassCat.index(g, 533 * u"nm"), 1.519417351519283, atol = 1e-14)
        @test isapprox(OpticSim.GlassCat.index(g, 533 * u"nm", temperature = 35 * u"°C"), 1.519462486258311, atol = 1e-14)
        @test isapprox(OpticSim.GlassCat.index(g, 533 * u"nm", pressure = 2.0), 1.518994119690216, atol = 1e-14)
        @test isapprox(OpticSim.GlassCat.index(g, 533 * u"nm", temperature = 35 * u"°C", pressure = 2.0), 1.519059871499476, atol = 1e-14)

        # test transmission
        @test OpticSim.GlassCat.absorption(OpticSim.GlassCat.Air, 500 * u"nm") == 0.0
        @test OpticSim.GlassCat.absorption(OpticSim.GlassCat.Air, 600 * u"nm") == 0.0
        # TODO these are currently taken from this package (i.e. regression tests), ideally we would get true values somehow
        @test OpticSim.GlassCat.absorption(g, 500 * u"nm") == 0.0002407228930225164
        @test OpticSim.GlassCat.absorption(g, 600 * u"nm") == 0.00022060722752440445
        @test OpticSim.GlassCat.absorption(g, 3000 * u"nm") == 0.04086604990127926
        @test OpticSim.GlassCat.absorption(g, 600 * u"nm", temperature = 35 * u"°C", pressure = 2.0) == 0.00022075540719494738

        # test that everything is alloc-less
        @test (@allocated OpticSim.GlassCat.absorption(g, 600 * u"nm")) == 0
        @test (@allocated OpticSim.GlassCat.index(g, 533 * u"nm")) == 0
        @test (@allocated OpticSim.GlassCat.SCHOTT.N_BK7) == 0

        @test OpticSim.GlassCat.glassforid(OpticSim.GlassCat.glassid(OpticSim.GlassCat.SCHOTT.N_BK7)) == OpticSim.GlassCat.SCHOTT.N_BK7

        @test OpticSim.GlassCat.isair(OpticSim.GlassCat.Air) == true
        @test OpticSim.GlassCat.isair(OpticSim.GlassCat.SCHOTT.N_BK7) == false

        # test MIL and model glasses
        fit_acc = 0.001
        bk7 = OpticSim.GlassCat.glassfromMIL(517642)
        @test OpticSim.GlassCat.glassforid(OpticSim.GlassCat.glassid(bk7)) == bk7
        @test OpticSim.GlassCat.isapprox(OpticSim.GlassCat.index(bk7, 0.5875618), 1.517, atol = fit_acc)
        @test OpticSim.GlassCat.isapprox(OpticSim.GlassCat.index(bk7, 0.533), 1.519417351519283, atol = fit_acc)
        @test OpticSim.GlassCat.isapprox(OpticSim.GlassCat.index(bk7, 0.743), 1.511997032563557, atol = fit_acc)
        @test OpticSim.GlassCat.isapprox(OpticSim.GlassCat.index(bk7, 533 * u"nm", temperature = 35 * u"°C", pressure = 2.0), 1.519059871499476, atol = fit_acc)

        t2 = OpticSim.GlassCat.glassfromMIL(1.135635)
        @test OpticSim.GlassCat.isapprox(OpticSim.GlassCat.index(t2, 0.5875618), 2.135, atol = fit_acc)

        fit_acc = 0.0001
        bk7 = OpticSim.GlassCat.modelglass(1.5168, 64.167336, -0.0009)
        @test OpticSim.GlassCat.glassforid(OpticSim.GlassCat.glassid(bk7)) == bk7
        @test isapprox(OpticSim.GlassCat.index(bk7, 0.5875618), 1.5168, atol = fit_acc)
        @test isapprox(OpticSim.GlassCat.index(bk7, 0.533), 1.519417351519283, atol = fit_acc)
        @test isapprox(OpticSim.GlassCat.index(bk7, 0.743), 1.511997032563557, atol = fit_acc)
        @test isapprox(OpticSim.GlassCat.index(bk7, 533 * u"nm", temperature = 35 * u"°C", pressure = 2.0), 1.519059871499476, atol = fit_acc)

        # test other glass
        @test OpticSim.GlassCat.index(OpticSim.GlassCat.CARGILLE.OG0608, 0.578) == 1.4596475735607324

        # make sure that the other functions work
        OpticSim.GlassCat.plot_indices(OpticSim.GlassCat.SCHOTT.N_BK7; polyfit = true, fiterror = true)
    end
end # testset GlassCat


text:
module Constants_module

export min_speed, max_speed, dt
export sensing_range, airspace_dim, dest_range, collision_range

const dt = 1. # sec
const min_speed = 163. # ft/s, 50 m/s
const max_speed = 163. # ft/s, 50 m/s
const sensing_range = 3280.8 # ft, 1 km
const airspace_dim = 32808.4 # ft, 10 km (65616.8 ft = 20 km)
const dest_range = 492. # ft, 150 m
const collision_range = 492. # ft, 150 m

end # module

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
# ---
# title: 826. Most Profit Assigning Work
# id: problem826
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Two Pointers
# link: <https://leetcode.com/problems/most-profit-assigning-work/description/>
# hidden: true
# ---
#
# We have jobs: `difficulty[i]` is the difficulty of the `i`th job, and
# `profit[i]` is the profit of the `i`th job.
#
# Now we have some workers. `worker[i]` is the ability of the `i`th worker,
# which means that this worker can only complete a job with difficulty at most
# `worker[i]`.
#
# Every worker can be assigned at most one job, but one job can be completed
# multiple times.
#
# For example, if 3 people attempt the same job that pays /$/1, then the total
# profit will be /$/3.  If a worker cannot complete any job, his profit is /$/0.
#
# What is the most profit we can make?
#
# **Example 1:**
#
#
#
#     Input: difficulty = [2,4,6,8,10], profit = [10,20,30,40,50], worker = [4,5,6,7]
#     Output: 100
#     **Explanation: W** orkers are assigned jobs of difficulty [4,4,6,6] and they get profit of [20,20,30,30] seperately.
#
# **Notes:**
#
#   * `1 <= difficulty.length = profit.length <= 10000`
#   * `1 <= worker.length <= 10000`
#   * `difficulty[i], profit[i], worker[i]`  are in range `[1, 10^5]`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
myDict=Dict{(@compat AbstractString), Int}()
myDict["Hi"] = 37
myDict["ho"] = 23
print(barplot(myDict))
myPlot=barplot(["Please","don't","crash","my","friend"], [10,24,30,13,7])
print(myPlot)
barplot!(myPlot, ["just", "dont", "!"], [2, 49, 15])
@test_throws DimensionMismatch barplot!(myPlot, ["just", "dont", "!"], [2, 49])
print(myPlot)
print(barplot([5,4,3,2,1], [10,24,30,13,7]))
print(barplot([5,4,3,2,1], [0,0,0,0,0]))
print(barplot([:Please,:dont,:crash,:dude], [10,24,1,13]))
print(barplot([:Please,:dont,:crash,:dude], [10,24,1,13], border=:none))
print(barplot([:Please,:dont,:crash,:dude], [10,24,1,13], border=:solid))
print(barplot([:Please,:dont,:crash,:dude], [10,24,1,13], border=:dashed))
print(barplot([:Please,:dont,:crash,:dude], [1,1,1,1000000], color=:red))
print(barplot([:Please,:dont,:crash,:dude], [10,24,1,13], width=10))
print(barplot([:Please,:dont,:crash,:dude], [10,24,1,13], width=10, margin=0))
print(barplot([:Please,:dont,:crash,:dude], [10,24,1,13], width=10, margin=20))
print(barplot([:Please,:dont,:crash,:dude], [10,24,1,13], width=1))
print(barplot([:Please,:dont,:crash,:dude], [10,24,1,0], width=100))
print(barplot([:Please,:dont,:crash,:dude], [10,24,1,0], title="C'mon man, keep on going!"))
print(barplot([:Please,:dont,:crash,:dude], [10,24,1,0], title="No lab", labels = false))
@test_throws ArgumentError barplot([:Please,:dont,:crash,:dude], [10,-1,1,1])
@test_throws DimensionMismatch barplot([:Please,:dont,:crash,:dude], [10,1,1])
@test_throws ArgumentError barplot([:Please,:dont,:crash,:dude], [10,24,1,0], margin=-1)
print(barplot([:Please,:dont,:crash,:dude], [1.,.7,.1,.6], width=10))

x = [-1.,2, 3, 7]
y = [1.,2, 9, 4]
print(lineplot(x, y))
x = [1.,2, 3, 7]
y = [1.,2, -1, 4]
print(lineplot(x, y))
x = [1,2, 3, 7]
y = [1,2, -1, 4]
@test_throws DimensionMismatch barplot(x, [10,1,1])
@test_throws DimensionMismatch lineplot(x, [10,1,1])
print(lineplot(1:10, 1:10))
print(lineplot(x, y))
@test_throws ArgumentError lineplot(x, y, margin=-1)
print(lineplot(x, y, width=50, height=10))
print(lineplot(x, y, width=100, height=10))
print(lineplot(x, y, width=90, height=40))
print(lineplot(x, y, width=50, height=5, margin = 0))
print(lineplot(x, y, width=10, height=20, margin = 20))
print(lineplot(x, y, width=5, height=5))
print(lineplot(x, y, width=50, height=10, title="Hello world", border=:none))
print(lineplot(x, y, width=50, height=10, title="Hello world", border=:dashed))
print(lineplot(x, y, width=50, height=10, title="Hello world", border=:solid))
print(lineplot(x, y, width=50, height=10, title="Hello world", border=:bold))
print(lineplot(x, y, width=50, height=10, title="Hello world", border=:dotted))
print(lineplot(x, y, width=50, height=10, border=:dotted))
x = [-1.,2, 3, 700000]
y = [1.,2, 9, 4000000]
print(lineplot(x, y))
print(lineplot(x, y, width=5, height=5))
@test_throws ArgumentError scatterplot(x, y, margin=-1)
x = [1.,2, 3, 7]
y = [1.,2, -1, 4]
print(scatterplot(x, y))
x = [1,2, 3, 7]
y = [1,2, -1, 4]
print(scatterplot(x, y))
x = [1.,2, 3, 7]
y = [1,2, -1, 4]
print(scatterplot(x, y))
print(scatterplot(x*.001+1.1, y))
print(lineplot(x, y*20+1000))
print(lineplot(x+1000, -y))

print(lineplot(sin, 1:.5:10))
print(lineplot(sin, 1:.5:10, labels = false))
print(lineplot(sin, [1., 1.5, 2, 2.5, 3, 3.5, 4]))
print(lineplot(sin, [1, 2, 3, 4]))
print(lineplot(sin, 1, 4))
print(lineplot([sin, cos], 1, 4))

x = [1,2, 4, 7, 8]
y = [1,3, 4, 2, 7]
print(stairs(x, y))
y = [1,3, 4, 2, 7000]
print(stairs(x, y))

#x=rand(100000)
#y=rand(100000)
#@time stairs(x, y)
#@time scatterplot(x, y)
#@time lineplot(x, y)

canvas = BrailleCanvas(40, 10, # number of columns and rows
                       plotOriginX = 0., plotOriginY = 0., # position in virtual space
                       plotWidth = 1., plotHeight = 1.) # size of the virtual space
drawLine!(canvas, 0., 0., 1., 1., :blue)
setPoint!(canvas, rand(50), rand(50), :red)
setPoint!(canvas, rand(50), rand(50), color = :red)
drawLine!(canvas, 0., 1., .5, 0., :yellow)

drawLine!(canvas, 0., 0., .9, 9999., :yellow)
drawLine!(canvas, 0., 0., 1., 1., color = :blue)
drawLine!(canvas, .3, .7, 1., 0., :red)
drawLine!(canvas, 0., 2., .5, 0., :yellow)
show(STDOUT, canvas)
myPlot = Plot(canvas, title="testtitle")
for i in 1:10
  annotate!(myPlot, :l, i, "$i", color=:green)
  annotate!(myPlot, :r, i, "$i", :yellow)
end
print(scatterplot([1],[1]))
@test_throws ArgumentError annotate!(myPlot, :bl, 5, ":l  5", :red)
annotate!(myPlot, :l, 5, ":l  5", :red)
annotate!(myPlot, :r, 5, "5  :r", color=:red)
annotate!(myPlot, :bl, ":bl", :blue)
annotate!(myPlot, :b, ":b", :green)
annotate!(myPlot, :br, ":br", :magenta)
annotate!(myPlot, :tl, ":tl", :cyan)
annotate!(myPlot, :tr, ":tr")
annotate!(myPlot, :t, ":t", :yellow)
drawLine!(myPlot, 0., 1., 1., 0., :blue)
setPoint!(myPlot, rand(10), rand(10))
setPixel!(myPlot, 1, 1)
setPixel!(myPlot, 1, 1, :green)
setPixel!(myPlot, 1, 1, color=:green)

ttl = "title!(plot, text)"
title!(myPlot, ttl)
@test title(myPlot) == ttl

xlab = "xlabel!(plot, text)"
xlabel!(myPlot, xlab)
@test xlabel(myPlot) == xlab

ylab = "ylabel!(plot, text)"
ylabel!(myPlot, ylab)
@test ylabel(myPlot) == ylab

show(STDOUT, myPlot)

x = [1,2, 4, 7, 8]
y = [1,3, 4, 2, 7]
myPlot = stairs(x, y, style=:pre)
print(myPlot)
myPlot = stairs(x, y)
print(myPlot)
stairs!(myPlot, x- .2, y + 1.5)
print(myPlot)

@test_throws ArgumentError lineplot(sin, 1:.5:12, color=:blue, ylim=[-1.,1., 2.])
myPlot = lineplot(sin, 1:.5:12, color=:blue, ylim=[-1.,1.])
print(myPlot)
lineplot!(myPlot, cos, 1:.5:10, color=:red)
print(myPlot)
myPlot = lineplot(sin, 1:.5:12, color=:blue, ylim=[-1,1], xlim=[-1.,5.])
print(myPlot)

x = [1,2, 4, 7, 8]
y = [1,3, 4, 2, 7]
myPlot = stairs(x, y, width = 10, padding = 3)
annotate!(myPlot, :tl, "Hello")
annotate!(myPlot, :t, "how are")
annotate!(myPlot, :tr, "you?")
annotate!(myPlot, :bl, "Hello")
annotate!(myPlot, :b, "how are")
annotate!(myPlot, :br, "you?")
print(myPlot)
lineplot!(myPlot, 1, .5)
print(myPlot)

print(histogram(rand(1000), bins=10, title="Histogram"))

print(spy(sprand(100,100,.15)))
print(spy(full(sprand(100,100,.15))))

x1, y1 = rand(500)*10, rand(500)*10
x2, y2 = rand(1000)*5+1, rand(1000)*5+1
canvas = BrailleCanvas(40, 20,
                       plotOriginX = 0., plotOriginY = 0.,
                       plotWidth = 10., plotHeight = 10.)
setPoint!(canvas, x1, y1, :red)
setPoint!(canvas, x2, y2, :blue)
show(canvas)
myPlot = densityplot(x1,y1, color = :red)
setPoint!(myPlot, x2, y2, :blue)
show(myPlot)
myPlot = densityplot(randn(1000), randn(1000), color = :blue)
densityplot!(myPlot, randn(1000) + 2, randn(1000) + 2, color = :red)
print(myPlot)

x = rand(100); y = rand(100)
print(scatterplot(x, y, xlim=[minimum(x), maximum(x)], ylim=[minimum(y), maximum(y)]))

x = rand(100) * 1000; y = rand(100) * 1000
print(scatterplot(x, y, xlim=[minimum(x), maximum(x)], ylim=[minimum(y), maximum(y)]))


text:
using LibTOAST
using LinearAlgebra
using Test

# Get some test data from the Toast++ source distribution
meshdir_2D = joinpath(LibTOAST._jl_toast_dir, "test", "2D", "meshes")
meshdir_3D = joinpath(LibTOAST._jl_toast_dir, "test", "3D", "meshes")

@testset "LibTOAST.jl Core" begin

  include("mesh.jl")

  global mesh2D = Mesh(joinpath(meshdir_2D, "circle25_32.msh"));
  global mesh3D = Mesh(joinpath(meshdir_2D, "circle25_32.msh"));

  include("sysmat.jl")
  include("assembly.jl")
  include("raster.jl")
  #include("coefftypes.jl")
  include("regul.jl")

end


text:
"""
    uri_add_path(::URI, ::String...)
Return a new `URI` with the original path extended by the
specified additional components.
"""
function uri_add_path(uri::URI, subs::String...)::URI
    path = ["", URIs.splitpath(uri)..., subs...]
    URI(uri; path=join(path, "/"))
end

"""
    fetch_toml(::URI)
Fetch the contents of the specified `URI` and parse it as TOML.
"""
function fetch_toml(uri::URI)
    response = HTTP.request("GET", uri)
    @assert response.status == 200
    TOML.parse(String(response.body))
end

max_stars_repo_name: MarkNahabedian/WebBasedWorkspace

max_issues_repo_name: MarkNahabedian/WebBasedWorkspace

max_forks_repo_name: MarkNahabedian/WebBasedWorkspace


text:
struct VisualizationData
  grid::Grid
  filebase::AbstractString
  celldata
  nodaldata
  function VisualizationData(grid::Grid,filebase::AbstractString;celldata=Dict(),nodaldata=Dict())
    new(grid,filebase,celldata,nodaldata)
  end
end

"""
This function returns an iterable collection (e.g. a Vector) of VisualizationData objects
"""
function visualization_data(args...;kwargs...)
  @abstractmethod
end

# Visualizing Triangulation, cellfields, etc.

function visualization_data(grid::Grid, filebase::AbstractString;celldata=Dict(),nodaldata=Dict())
  (VisualizationData(grid,filebase;celldata=celldata,nodaldata=nodaldata),)
end

function visualization_data(
  trian::Triangulation, filebase::AbstractString;
  order=-1, nsubcells=-1, celldata=Dict(), cellfields=Dict())

  if order == -1 && nsubcells == -1
    # Use the given cells as visualization cells
    f = (reffe) -> UnstructuredGrid(reffe)
  elseif order != -1 && nsubcells == -1
    # Use cells of given order as visualization cells
    f = (reffe) -> UnstructuredGrid(LagrangianRefFE(Float64,get_polytope(reffe),order))
  elseif order == -1 && nsubcells != -1
    # Use linear sub-cells with nsubcells per direction
    f = (reffe) -> UnstructuredGrid(compute_reference_grid(reffe,nsubcells))
  else
    @unreachable "order and nsubcells kw-arguments can not be given at the same time"
  end

  ref_grids = map(f, get_reffes(trian))
  visgrid = VisualizationGrid(trian,ref_grids)

  cdata = _prepare_cdata(celldata,visgrid.sub_cell_to_cell)
  pdata = _prepare_pdata(trian,cellfields,visgrid.cell_to_refpoints)

  (VisualizationData(visgrid,filebase;celldata=cdata,nodaldata=pdata),)
end

function _prepare_node_to_coords(cell_to_points)
  cell_to_offset = zeros(Int,length(cell_to_points))
  P = eltype(eltype(cell_to_points))
  node_to_coords = P[]
  cache = array_cache(cell_to_points)
  _prepare_node_to_coords_fill!(node_to_coords,cell_to_offset,cache,cell_to_points)
  (node_to_coords, cell_to_offset)
end

function _prepare_node_to_coords_fill!(node_to_coords,cell_to_offset,cache,cell_to_points)
  offset = 0
  for cell in 1:length(cell_to_points)
    cell_to_offset[cell] = offset
    points = getindex!(cache,cell_to_points,cell)
    for point in points
      push!(node_to_coords, point)
      offset += 1
    end
  end
end

function _prepare_sub_cell_to_nodes(
  cell_to_ctype,ctype_to_scell_to_snodes,cell_to_offset)

  sub_cell_to_nodes = Table(Int[],Int32[1])
  sub_cell_to_cell = Int[]

  for (cell, ctype) in enumerate(cell_to_ctype)
    scell_to_snodes = ctype_to_scell_to_snodes[ctype]
    offset = cell_to_offset[cell]
    for (scell, snodes) in enumerate(scell_to_snodes)
      push!(sub_cell_to_cell,cell)
      push!(sub_cell_to_nodes.ptrs,length(snodes))
      for snode in snodes
        node = snode+offset
        push!(sub_cell_to_nodes.data,node)
      end
    end
  end

  length_to_ptrs!(sub_cell_to_nodes.ptrs)

  (sub_cell_to_nodes, sub_cell_to_cell)
end

function _prepare_sctype_to_reffe(ctype_to_r_to_reffe,ctype_to_scell_to_r,cell_to_ctype)

  i_to_reffe = vcat(ctype_to_r_to_reffe...)
  u_to_reffe, i_to_u = _find_unique_with_indices(i_to_reffe)
  ctype_to_nreffes = map(length,ctype_to_r_to_reffe)
  sub_cell_to_u = _prepare_sub_cell_to_u(
    ctype_to_r_to_reffe,ctype_to_nreffes,ctype_to_scell_to_r,cell_to_ctype,i_to_u)

  (u_to_reffe, sub_cell_to_u)
end

function _prepare_sub_cell_to_u(
  ctype_to_r_to_reffe,ctype_to_nreffes,ctype_to_scell_to_r,cell_to_ctype,i_to_u)

  i = 1
  ctype_to_r_to_u = Vector{Int}[]
  for (ctype, r_to_reffes) in enumerate(ctype_to_r_to_reffe)
    nreffes = ctype_to_nreffes[ctype]
    r_to_u = zeros(Int,nreffes)
    for r in 1:nreffes
      r_to_u[r] = i_to_u[i]
      i +=1
    end
    push!(ctype_to_r_to_u,r_to_u)
  end

  sub_cell_to_u = Int8[]
  for ctype in cell_to_ctype
    scell_to_r = ctype_to_scell_to_r[ctype]
    r_to_u = ctype_to_r_to_u[ctype]
    for r in scell_to_r
      push!(sub_cell_to_u,r_to_u[r])
    end
  end

  sub_cell_to_u

end

function _prepare_pdata(trian,cellfields,samplingpoints)
  x = CellPoint(samplingpoints,trian,ReferenceDomain())
  pdata = Dict()
  for (k,v) in cellfields
    _v = CellField(v,trian)
    pdata[k], = _prepare_node_to_coords(evaluate(_v,x))
  end
  pdata
end

function _prepare_cdata(celldata,sub_cell_to_cell)
  cdata = Dict()
  for (k,v) in celldata
    cell_to_a = collect(get_array(v))
    T = eltype(cell_to_a)
    sub_cell_to_a = zeros(T,length(sub_cell_to_cell))
    _fill_sub_cell_to_a!(sub_cell_to_a,cell_to_a,sub_cell_to_cell)
    cdata[k] = sub_cell_to_a
  end
  k2 = "cell"
  @assert ! haskey(cdata,k2) "cell is a reserved key"
  cdata[k2] = sub_cell_to_cell
  cdata
end

function _fill_sub_cell_to_a!(sub_cell_to_a,cell_to_a,sub_cell_to_cell)
  for (sub_cell, cell) in enumerate(sub_cell_to_cell)
    sub_cell_to_a[sub_cell] = cell_to_a[cell]
  end
end

struct VisualizationGrid{Dc,Dp} <: Grid{Dc,Dp}
  sub_grid::UnstructuredGrid{Dc,Dp}
  sub_cell_to_cell::AbstractVector{<:Integer}
  cell_to_refpoints::AbstractVector{<:AbstractVector{<:Point}}
end

get_reffes(g::VisualizationGrid) = get_reffes(g.sub_grid)

get_cell_type(g::VisualizationGrid) = get_cell_type(g.sub_grid)

get_node_coordinates(g::VisualizationGrid) = get_node_coordinates(g.sub_grid)

get_cell_node_ids(g::VisualizationGrid) = get_cell_node_ids(g.sub_grid)

function VisualizationGrid(trian::Triangulation, ref_grids::AbstractArray{<:UnstructuredGrid})

  cell_to_ctype = collect1d(get_cell_type(trian))
  ctype_to_refpoints = map(get_node_coordinates, ref_grids)
  cell_to_refpoints = CompressedArray(ctype_to_refpoints,cell_to_ctype)
  cell_map = get_cell_map(trian)
  cell_to_points = lazy_map(evaluate,cell_map, cell_to_refpoints)

  node_to_coords, cell_to_offset = _prepare_node_to_coords(cell_to_points)

  ctype_to_scell_to_snodes = map(get_cell_node_ids,ref_grids)

  sub_cell_to_nodes, sub_cell_to_cell = _prepare_sub_cell_to_nodes(
    cell_to_ctype,ctype_to_scell_to_snodes,cell_to_offset)

  ctype_to_reffes = map(get_reffes,ref_grids)
  ctype_to_scell_type = map(get_cell_type,ref_grids)
  sctype_to_reffe, sub_cell_to_sctype = _prepare_sctype_to_reffe(
    ctype_to_reffes,ctype_to_scell_type,cell_to_ctype)

  sub_grid = UnstructuredGrid(
    node_to_coords,
    sub_cell_to_nodes,
    sctype_to_reffe,
    sub_cell_to_sctype,
    NonOriented())

  VisualizationGrid(sub_grid,sub_cell_to_cell,cell_to_refpoints)

end

function visualization_data(p::Polytope,filebase)
  map(0:(num_dims(p)-1)) do d
    filebase_d = "$(filebase)_$d"
    grid = Grid(ReferenceFE{d},p)
    VisualizationData(grid,filebase_d)
  end
end

function visualization_data(x::AbstractVector{<:Point}, filebase; kwargs...)
  grid = UnstructuredGrid(x)
  (VisualizationData(grid,filebase; kwargs...),)
end

function visualization_data(reffe::LagrangianRefFE,filebase)
  p = get_polytope(reffe)
  visdata_p = visualization_data(p,filebase)
  node_coords = get_node_coordinates(reffe)
  node_comp_to_dof = get_node_and_comp_to_dof(reffe)
  nodaldata = [
    "dof" => node_comp_to_dof,
    "node" => collect(1:num_nodes(reffe))]
  visdata_n = visualization_data(node_coords, "$(filebase)_nodes"; nodaldata = nodaldata)
  [visdata_p..., visdata_n...]
end

function visualization_data(
  model::DiscreteModel, filebase::AbstractString;labels::FaceLabeling=get_face_labeling(model))
  map(0:num_cell_dims(model)) do d
    grid = Grid(ReferenceFE{d},model)
    cdat = _prepare_cdata_model(labels,d)
    VisualizationData(grid,"$(filebase)_$d";celldata=cdat)
  end
end

function _prepare_cdata_model(labels,d)
  dface_to_entity = get_face_entity(labels,d)
  cdat = []
  for tag in 1:num_tags(labels)
    dface_to_isontag = zeros(Int,num_faces(labels,d))
    for entity in get_tag_entities(labels,tag)
      _set_entity!(dface_to_isontag,dface_to_entity,entity,tag)
    end
    name = get_tag_name(labels,tag)
    push!(cdat, name => dface_to_isontag )
  end
  push!(cdat,"entity" => dface_to_entity)
  cdat
end

function _set_entity!(dface_to_isontag,dface_to_entity,entity,tag)
  for i in 1:length(dface_to_entity)
    if dface_to_entity[i] == entity
      dface_to_isontag[i] = tag
    end
  end
end

function visualization_data(
  cell_to_points::AbstractArray{<:AbstractArray{<:Point}},
  filename; celldata=Dict(), nodaldata=Dict())

  node_to_point, cell_to_offset = _prepare_node_to_coords(cell_to_points)
  nnodes = length(node_to_point)
  cell_to_nodes = identity_table(Int,Int32,nnodes)
  cell_to_ctype = fill(Int8(1),nnodes)
  ctype_to_reffe = [VERTEX1,]
  node_to_cell = _prepare_node_to_cell(cell_to_offset,nnodes)
  grid = UnstructuredGrid(
    node_to_point,
    cell_to_nodes,
    ctype_to_reffe,
    cell_to_ctype,
    Oriented())
  cdata = _prepare_cdata(celldata,node_to_cell)
  pdata = _prepare_pdata_for_cell_points(nodaldata)

  (VisualizationData(grid,filename,celldata=cdata,nodaldata=pdata),)
end

function _prepare_node_to_cell(cell_to_offset,nnodes)
  node_to_cell = zeros(Int,nnodes)
  ncells = length(cell_to_offset)
  for cell in 1:ncells
    nini = cell_to_offset[cell]+1
    if cell < ncells
      nend = cell_to_offset[cell+1]
    else
      nend = nnodes
    end
    for node in nini:nend
      node_to_cell[node] = cell
    end
  end
  node_to_cell
end

function _prepare_pdata_for_cell_points(nodaldata)
  pdata = Dict()
  for (k,v) in nodaldata
    pdata[k], = _prepare_node_to_coords(v)
  end
  pdata
end

function visualization_data(x::CellPoint,filename; cellfields=Dict())
  nodaldata = Dict((
    k=>evaluate(CellField(v,get_triangulation(x),DomainStyle(x)),x) for (k,v) in cellfields ))
  cell_to_points = get_array(x)
  visualization_data(cell_to_points, filename, nodaldata=nodaldata)
end


text:
const DesignResolution = (1080, 1920)

frame_time(scene::AbstractPokemonScene) = scene.frame_time

function get_current_context!(ctx::PokemonContext)::BattleContext
    data = ctx.data
    parsing_context = data.context
    if isnothing(parsing_context.current)
        current = parsing_context.current = BattleContext()
        push!(parsing_context.battles, current.battle)
        if ctx.config.source.parse_type > PARSE_NONE
            current.parsed_battle = ParsedBattle()
            push!(parsing_context.parsed_battles, current.parsed_battle)
        end
    end
    parsing_context.current
end

_parse_scene(::Type{<:AbstractPokemonScene}, ::VsFrame, ::PokemonContext) = nothing
VsRecorderBase.vs_tryparse_scene(
    T::Type{<:AbstractPokemonScene},
    frame::VsFrame,
    ctx::Union{PokemonContext, Nothing} = nothing;
    force::Bool = false
) = try
    if isnothing(ctx)
        ctx = default_context()
    end
    if ctx.current_frame ≢ frame
        ctx.current_frame = frame
    end
    ctx.data.force = force
    _parse_scene(T, frame, ctx)
catch e
    ctx.data.last_error = e
    rethrow(e)
    nothing
end

VsRecorderBase.vs_update!(ctx::PokemonContext, scene::AbstractPokemonScene) = _vs_update!(ctx, scene)

should_include(scene_type::Type{<:AbstractPokemonScene}, ctx::PokemonContext) = true


text:
# ---
# title: 1002. Find Common Characters
# id: problem1002
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array, Hash Table
# link: <https://leetcode.com/problems/find-common-characters/description/>
# hidden: true
# ---
#
# Given an array `A` of strings made only from lowercase letters, return a list
# of all characters that show up in all strings within the list **(including
# duplicates)**.  For example, if a character occurs 3 times in all strings but
# not 4 times, you need to include that character three times in the final
# answer.
#
# You may return the answer in any order.
#
#
#
# **Example 1:**
#
#
#
#     Input: ["bella","label","roller"]
#     Output: ["e","l","l"]
#
#
# **Example 2:**
#
#
#
#     Input: ["cool","lock","cook"]
#     Output: ["c","o"]
#
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 100`
#   2. `1 <= A[i].length <= 100`
#   3. `A[i][j]` is a lowercase letter
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Gadfly

DF68 = 8 # свободное расстояние кода, вариант 68

FCoefs = [3, 15, 58, 201, 655, 2052, 6255, 18687, 54974, 159765]
F68Powers = [8:2:26]

function Q(x)
    0.5 * erfc(x / sqrt(2))
end

function Peb(df, B, D0, F, p0)
    B(df, p0) * F(D0(p0))
end

function dskB(w, p0)
    (1 - p0) / (1 - 2p0) * sqrt(2 / (π * w))
end

function dskD0(p0)
    2 * sqrt(p0 * (1 - p0))
end

function abgshB(w, sg)
    Q(sqrt(2 * w * sg)) * exp(w * sg)
end

function abgshD0(sg)
    exp(-sg)
end

function F68(D)
    (D^8 * (D^4 - 3D^2 + 3)) / (D^4 - 3D^2 + 1)^2
end

function F68_5(D)
    dot(FCoefs[1:5], map(i -> D^i, F68Powers[1:5]))
end

function F68_10(D)
    dot(FCoefs, map(i -> D^i, F68Powers))
end

sdb = [-2:0.01:10] # SNR, dB per bit
sb = 10 .^ (sdb / 10) # SNR per bit
sig = sqrt(0.5 ./ sb)
sig1 = 1 ./ sig
p0 = Q(sig1)

dskPEB68      = map(p -> Peb(DF68, dskB, dskD0, F68, p), p0)
dskPEB68_5    = map(p -> Peb(DF68, dskB, dskD0, F68_5, p), p0)
dskPEB68_10   = map(p -> Peb(DF68, dskB, dskD0, F68_10, p), p0)

abgshPEB68      = map(p -> Peb(DF68, abgshB, abgshD0, F68, p), sb)
abgshPEB68_5    = map(p -> Peb(DF68, abgshB, abgshD0, F68_5, p), sb)
abgshPEB68_10   = map(p -> Peb(DF68, abgshB, abgshD0, F68_10, p), sb)

gplot = plot(layer(x=sdb, y=dskPEB68,      Theme(default_color=colorant"blue"), Geom.line),
             layer(x=sdb, y=dskPEB68_5,    Theme(default_color=colorant"red"), Geom.line),
             layer(x=sdb, y=dskPEB68_10,   Theme(default_color=colorant"green"), Geom.line),
             layer(x=sdb, y=abgshPEB68,    Theme(default_color=colorant"orange"), Geom.line),
             layer(x=sdb, y=abgshPEB68_5,  Theme(default_color=colorant"black"), Geom.line),
             layer(x=sdb, y=abgshPEB68_10, Theme(default_color=colorant"purple"), Geom.line),
             Guide.manual_color_key("Графики",
                                   ["ДСК, F", "ДСК, 5", "ДСК, 10", "АБГШ, F", "АБГШ, 5", "АБГШ, 10"],
                                   ["blue", "red", "green", "orange", "black", "purple"]),
             Theme(default_color=colorant"black"),
             Guide.xticks(ticks=[-2:1:10]), Scale.y_log10,
             Guide.xlabel("E/N0 (dB)"), Guide.ylabel("Вероятность ошибки на бит"))

draw(PDF("errors_all.pdf", 6inch, 4inch), gplot)

max_stars_repo_name: ChShersh/university-courses

max_issues_repo_name: ChShersh/university-courses

max_forks_repo_name: ChShersh/university-courses


text:
function genes_helper!(ex)
    if ex isa Expr && !isempty(ex.args)
        if ex.args[1] ∉ [:Ref, :upstream, :downstream, :neighbours]
            skipgene = ex.args[1] ∉ [:(||), :(&&)]
            skipfirst = ex.head == :call
            genes_special_cases!(ex.args; skipgene = skipgene, skipfirst = skipfirst)
        end
        if ex.head == :call
            if ex.args[1] == :Ref
                # Remove Ref(), but ignore the content
                ex2 = copy(ex)
                ex.head = :ref
                empty!(ex.args)
                push!(ex.args, ex2)
            elseif ex.args[1] == :get
                if length(ex.args) == 3
                    # Expand `get(s::Symbol, default)` to `get(gene, s::Symbol, default)`
                    genes_helper!(ex.args[3])
                    ex.args = vcat(ex.args[1], :gene, ex.args[2:3])
                else
                    genes_helper!(ex.args[2])
                    genes_helper!(ex.args[4])
                end
            elseif ex.args[1] in [:upstream, :downstream, :neighbours]
                # Expand `upstream(i, f)` to `upstream(gene, f, i)`, etc.
                if length(ex.args) == 3
                    feature_function!(ex.args)
                    ex.args = vcat(ex.args[1], :gene, ex.args[2:3])
                end
            else
                for i in eachindex(ex.args)[2:end]
                    if ex.args[i] isa QuoteNode
                        ex.args[i] = :(Base.getproperty(gene, $(ex.args[i])))
                    else
                        genes_helper!(ex.args[i])
                    end
                end
            end
        elseif ex.head == :(.) && ex.args[1] == :get
            if length(ex.args[2].args) == 2
                genes_helper!(ex.args[2].args[2])
                ex.args[2].args = vcat(Expr(:call, :Ref, :gene), ex.args[2].args)
            else
                genes_helper!(ex.args[2].args[1])
                genes_helper!(ex.args[2].args[3])
            end
        else
            for i in eachindex(ex.args)
                if ex.args[i] isa QuoteNode
                    if ex.head != :(.)
                        ex.args[i] = :(Base.getproperty(gene, $(ex.args[i])))
                    end
                else
                    genes_helper!(ex.args[i])
                end
            end
        end
    end
end

"""
    @genes(chr, exs...)

Iterate over and evaluate expressions in `exs` for all genes in `chr.genes`,
returning genes where all expressions evaluate to `true`. Any given symbol `s`
in the expression will be substituted for `gene.s`. The gene itself can be
accessed in the expression as `gene`. Accessing properties of the returned list
of genes returns a view, which can be altered.

Symbols and expressions escaped with `Ref()` will be ignored.

Some short-hand forms are available to make life easier:
    `CDS`, `rRNA`, and `tRNA` expand to `feature(gene) == "..."`,
    `get(s::Symbol, default)` expands to `get(gene, s, default)`

# Examples
```julia
julia> chromosome = readgbk("example.gbk")
Chromosome 'example' (5028 bp) with 6 annotations

julia> @genes(chromosome, iscds) |> length
3

julia> @genes(chromosome, length(gene) < 500)
     CDS             3..206
                     /db_xref="GI:1293614"
                     /locus_tag="tag01"
                     /codon_start="3"
                     /product="TCP1-beta"
                     /protein_id="AAA98665.1"

julia> @genes(chromosome, ismissing(:gene)) |> length
2

julia> @genes(chromosome, ismissing(:gene)).gene .= "Unknown";

julia> @genes(chromosome, ismissing(:gene)) |> length
0
```

All arguments have to evaluate to `true` for a gene to be included, so the following expressions are equivalent:
```julia
@genes(chr, feature(gene) == Ref(:CDS), length(gene) > 300)
@genes(chr, (feature(gene) == Ref(:CDS)) && (length(gene) > 300))
```

`@genes` returns a `Vector{Gene}`. Attributes can be accessed with dot-syntax, and can be assigned to
```julia
@genes(chr, :locus_tag == "tag03")[1].pseudo = true
@genes(chr, CDS, ismissing(:gene)).gene .= "unknown"
```
"""
macro genes(chr, exs...)
    exs = Any[ex for ex in exs]
    genes_special_cases!(exs; skipgene = false, skipfirst = false)
    for (i, ex) in enumerate(exs)
        # Expressions
        if ex isa Expr
            exs[i] = :((x -> ismissing(x) ? false : x)($(exs[i])))
            genes_helper!(exs[i])
        elseif ex isa QuoteNode
            exs[i] = :((x -> ismissing(x) ? false : x)(Base.getproperty(gene, $ex)))
        end
    end
    ex = Expr(:&&, exs...)
    hits = gensym()
    i = gensym()
    return esc(quote
        if $chr isa AbstractVector
            vcat([begin
                    $hits = falses(size(c.genes, 1))
                    for ($i, gene) in enumerate(c.genes)
                        local h = $ex
                        $hits[$i] = ismissing(h) ? false : h
                    end
                    c.genes[$hits]
                end for c in $chr]...)
        else
            $hits = falses(size($chr.genes, 1))
            for ($i, gene) in enumerate($chr.genes)
                local h = $ex
                $hits[$i] = ismissing(h) ? false : h
            end
            $chr.genes[$hits]
        end
    end)
end

upstream(f::Function, gene, cond::Function, i) = f(upstream(gene, cond, i))
upstream(gene, cond::Function, i::Int) = upstream(gene, cond, i:i)[1]
function upstream(gene, cond::Function, i)
    ngenes = length(parent(gene).genes)
    j = 0
    genes = Gene[]
    grange = iscomplement(gene) ?
        parent(gene).genes[mod1.(index(gene) + 1 : index(gene) + ngenes - 2, ngenes)] :
        parent(gene).genes[reverse(mod1.((index(gene) + 1) : (index(gene) + ngenes - 2), ngenes))]
    for g in grange
        if cond(g)
            j += 1
            j in i && push!(genes, g)
        end
    end
    return genes
end

downstream(f::Function, gene, cond::Function, i) = f(downstream(gene, cond, i))
downstream(gene, cond::Function, i::Int) = downstream(gene, cond, i:i)[1]
function downstream(gene, cond::Function, i)
    ngenes = length(parent(gene).genes)
    j = 0
    genes = Gene[]
    grange = !iscomplement(gene) ?
        parent(gene).genes[mod1.(index(gene) + 1 : index(gene) + ngenes - 2, ngenes)] :
        parent(gene).genes[reverse(mod1.((index(gene) + 1) : (index(gene) + ngenes - 2), ngenes))]
    for g in grange
        if cond(g)
            j += 1
            j in i && push!(genes, g)
        end
    end
    return genes
end

neighbours(f::Function, gene, cond::Function, i) = f(neighbours(gene, cond, i))
function neighbours(gene, cond::Function, i)
    vcat(upstream(gene, cond, i), downstream(gene, cond, i))
end

function feature_function!(exs; skipgene = false, skipfirst = true)
    i = 2
    ex = exs[i]
    if ex isa Symbol
        if ex == :CDS
            exs[i] = :(GenomicAnnotations.feature(g) == :CDS)
        elseif ex == :(!CDS)
            exs[i] = :(GenomicAnnotations.feature(g) != :CDS)
        elseif ex == :rRNA
            exs[i] = :(GenomicAnnotations.feature(g) == :rRNA)
        elseif ex == :(!rRNA)
            exs[i] = :(GenomicAnnotations.feature(g) != :rRNA)
        elseif ex == :tRNA
            exs[i] = :(GenomicAnnotations.feature(g) == :tRNA)
        elseif ex == :(!tRNA)
            exs[i] = :(GenomicAnnotations.feature(g) != :tRNA)
        elseif ex == :gene && !skipgene
            exs[i] = :(GenomicAnnotations.feature(g) == :gene)
        elseif ex == :(!gene)
            exs[i] = :(GenomicAnnotations.feature(g) != :gene)
        elseif ex == :regulatory
            exs[i] = :(GenomicAnnotations.feature(g) == :regulatory)
        elseif ex == :(!regulatory)
            exs[i] = :(GenomicAnnotations.feature(g) != :regulatory)
        elseif ex == :source
            exs[i] = :(GenomicAnnotations.feature(g) == :source)
        elseif ex == :(!source)
            exs[i] = :(GenomicAnnotations.feature(g) != :source)
        end
        exs[i] = Expr(:(->), :g, Expr(:block, exs[i]))
    end
end

function genes_special_cases!(exs; skipgene = false, skipfirst = true)
    features = [:CDS, :rRNA, :tRNA, :mRNA, :exon, :intron, :regulatory, :source, :repeat_region]
    for (i, ex) in enumerate(exs)
        skipfirst && i == 1 && continue
        if ex isa Symbol
            if ex == :gene && !skipgene
                exs[i] = :(GenomicAnnotations.feature(gene) == Ref(:gene))
            elseif ex in features
                exs[i] = :(GenomicAnnotations.feature(gene) == Ref($(QuoteNode(ex))))
            end
        elseif ex isa Expr && ex.head == :call && ex.args[1] == :(!)
            feature = ex.args[2]
            if feature == :gene
                exs[i] = :(GenomicAnnotations.feature(gene) != Ref(:gene))
            elseif feature in features
                exs[i] = :(GenomicAnnotations.feature(gene) != Ref($(QuoteNode(feature))))
            end
        end
    end
end

# Features from the GenBank specifications without shorthand forms:
#   assembly_gap
#   C_region
#   centromere
#   D-loop
#   D_segment
#   gap
#   iDNA
#   J_segment
#   mat_peptide
#   misc_binding
#   misc_difference
#   misc_feature
#   misc_recomb
#   misc_RNA
#   misc_structure
#   mobile_element
#   modified_base
#   ncRNA
#   N_region
#   old_sequence
#   operon
#   oriT
#   polyA_site
#   precursor_RNA
#   prim_transcript
#   primer_bind
#   propeptide
#   protein_bind
#   rep_origin
#   S_region
#   sig_peptide
#   stem_loop
#   STS
#   telomere
#   tmRNA
#   transit_peptide
#   unsure
#   V_region
#   V_segment
#   variation
#   3'UTR
#   5'UTR


text:
module NeuroStudies

using DataFrames
using JSON3
using Metadata
using OrderedCollections
using Tables
using Base: @propagate_inbounds

export study, subject, subpaths, session, modality, derived, file

# TODO:
# - https://bids-specification.readthedocs.io/en/stable/99-appendices/10-file-collections.html
# - creating new file with similar entity names
# - Replace DataFrames with `Tables.dictrowtable`

include("entities.jl")
include("types.jl")
include("vectors.jl")
include("paths.jl")
include("utils.jl")

# include("data_type_entities.jl")
# ∘

end