text:
# # Density Current
#
# In this example, we demonstrate the usage of the `ClimateMachine`
#  to solve the density current test by Straka 1993.
# We solve a flow in a box configuration, which is
# representative of a large-eddy simulation. Several versions of the problem
# setup may be found in literature, but the general idea is to examine the
# vertical ascent of a thermal _bubble_ (we can interpret these as simple
# representation of convective updrafts).
#
# ## Description of experiment
# The setup described below is such that the simulation reaches completion
# (timeend = 900 s) in approximately 4 minutes of wall-clock time on 1 GPU
#
# 1) Dry Density Current (circular potential temperature perturbation)
# 2) Boundaries
#    - `Impenetrable(FreeSlip())` - no momentum flux, no mass flux through
#      walls.
#    - `Impermeable()` - non-porous walls, i.e. no diffusive fluxes through
#       walls.
# 3) Domain - 25600m (horizontal) x 10000m (horizontal) x 6400m (vertical)
# 4) Resolution - 100m effective resolution
# 5) Total simulation time - 900s
# 6) Mesh Aspect Ratio (Effective resolution) 1:1
# 7) Overrides defaults for
#    - CPU Initialisation
#    - Time integrator
#    - Sources
#    - Smagorinsky Coefficient _C_smag
# 8) Default settings can be found in `src/Driver/<files>.jl`

#md # !!! note
#md #     This experiment setup assumes that you have installed the
#md #     `ClimateMachine` according to the instructions on the landing page.
#md #     We assume the users' familiarity with the conservative form of the
#md #     equations of motion for a compressible fluid
#md #
#md #     The following topics are covered in this example
#md #     - Package requirements
#md #     - Defining a `model` subtype for the set of conservation equations
#md #     - Defining the initial conditions
#md #     - Applying boundary conditions
#md #     - Applying source terms
#md #     - Choosing a turbulence model
#md #     - Adding tracers to the model
#md #     - Choosing a time-integrator
#md #
#md #     The following topics are not covered in this example
#md #     - Defining new boundary conditions
#md #     - Defining new turbulence models
#md #     - Building new time-integrators
#
# ## Boilerplate (Using Modules)
#
# #### [Skip Section](@ref init-dc)
#
# Before setting up our experiment, we recognize that we need to import some
# pre-defined functions from other packages. Julia allows us to use existing
# modules (variable workspaces), or write our own to do so.  Complete
# documentation for the Julia module system can be found
# [here](https://docs.julialang.org/en/v1/manual/modules/#).

# We need to use the `ClimateMachine` module! This imports all functions
# specific to atmospheric and ocean flow modeling.  While we do not cover the
# ins-and-outs of the contents of each of these we provide brief descriptions
# of the utility of each of the loaded packages.
using ClimateMachine
ClimateMachine.init(parse_clargs = true)

using ClimateMachine.Atmos
using ClimateMachine.Orientations
# - Required so that we inherit the appropriate model types for the large-eddy
#   simulation (LES) and global-circulation-model (GCM) configurations.
using ClimateMachine.ConfigTypes
# - Required so that we may define diagnostics configurations, e.g. choice of
#   file-writer, choice of output variable sets, output-frequency and directory,
using ClimateMachine.Diagnostics
# - Required so that we may define (or utilise existing functions) functions
#   that are `called-back` or executed at frequencies of either timesteps,
#   simulation-time, or wall-clock time.
using ClimateMachine.GenericCallbacks
# - Required so we load the appropriate functions for the time-integration
#   component. Contains ODESolver methods.
using ClimateMachine.ODESolvers
# - Required for utility of spatial filtering functions (e.g. positivity
#   preservation)
using ClimateMachine.Mesh.Filters
# - Required so functions for computation of temperature profiles.
using ClimateMachine.TemperatureProfiles
# - Required so functions for computation of moist thermodynamic quantities and turbulence closures
# are available.
using ClimateMachine.Thermodynamics
using ClimateMachine.TurbulenceClosures
# - Required so we may access our variable arrays by a sensible naming
#   convention rather than by numerical array indices.
using ClimateMachine.VariableTemplates
# - Required so we may access planet parameters
#   ([CLIMAParameters](https://github.com/CliMA/CLIMAParameters.jl)
#   specific to this problem include the gas constant, specific heats,
#   mean-sea-level pressure, gravity and the Smagorinsky coefficient)

# In ClimateMachine we use `StaticArrays` for our variable arrays.
using StaticArrays
# We also use the `Test` package to help with unit tests and continuous
# integration systems to design sensible tests for our experiment to ensure new
# / modified blocks of code don't damage the fidelity of the physics. The test
# defined within this experiment is not a unit test for a specific
# subcomponent, but ensures time-integration of the defined problem conditions
# within a reasonable tolerance. Immediately useful macros and functions from
# this include `@test` and `@testset` which will allow us to define the testing
# parameter sets.
using Test

using CLIMAParameters
using CLIMAParameters.Atmos.SubgridScale: C_smag
using CLIMAParameters.Planet: R_d, cp_d, cv_d, MSLP, grav
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

# ## [Initial Conditions](@id init-dc)
#md # !!! note
#md #     The following variables are assigned in the initial condition
#md #     - `state.ρ` = Scalar quantity for initial density profile
#md #     - `state.ρu`= 3-component vector for initial momentum profile
#md #     - `state.ρe`= Scalar quantity for initial total-energy profile
#md #       humidity
#md #     - `state.tracers.ρχ` = Vector of four tracers (here, for demonstration
#md #       only; we can interpret these as dye injections for visualisation
#md #       purposes)
function init_densitycurrent!(problem, bl, state, aux, (x, y, z), t)
    ## Problem float-type
    FT = eltype(state)

    ## Unpack constant parameters
    R_gas::FT = R_d(bl.param_set)
    c_p::FT = cp_d(bl.param_set)
    c_v::FT = cv_d(bl.param_set)
    p0::FT = MSLP(bl.param_set)
    _grav::FT = grav(bl.param_set)
    γ::FT = c_p / c_v

    ## Define bubble center and background potential temperature
    xc::FT = 0
    yc::FT = 0
    zc::FT = 3000
    rx::FT = 4000
    rz::FT = 2000
    r = sqrt(((x - xc)^2) / rx^2 + ((z - zc)^2) / rz^2)

    ## TODO: clean this up, or add convenience function:
    ## This is configured in the reference hydrostatic state
    θ_ref::FT = bl.ref_state.virtual_temperature_profile.T_surface
    Δθ::FT = 0
    θamplitude::FT = -15.0

    ## Compute temperature difference over bubble region
    if r <= 1
        Δθ = 0.5 * θamplitude * (1 + cospi(r))
    end

    ## Compute perturbed thermodynamic state:
    θ = θ_ref + Δθ                                      ## potential temperature
    π_exner = FT(1) - _grav / (c_p * θ) * z             ## exner pressure
    ρ = p0 / (R_gas * θ) * (π_exner)^(c_v / R_gas)      ## density
    T = θ * π_exner
    e_int = internal_energy(bl.param_set, T)
    ts = PhaseDry(bl.param_set, e_int, ρ)
    ρu = SVector(FT(0), FT(0), FT(0))                   ## momentum
    ## State (prognostic) variable assignment
    e_kin = FT(0)                                       ## kinetic energy
    e_pot = gravitational_potential(bl.orientation, aux)## potential energy
    ρe_tot = ρ * total_energy(e_kin, e_pot, ts)         ## total energy

    ## Assign State Variables
    state.ρ = ρ
    state.ρu = ρu
    state.ρe = ρe_tot
end

# ## [Model Configuration](@id config-helper)
# We define a configuration function to assist in prescribing the physical
# model.
function config_densitycurrent(FT, N, resolution, xmax, ymax, zmax)

    ## Choose an Explicit Single-rate Solver LSRK144 from the existing [ODESolvers](@ref
    ## ODESolvers-docs) options Apply the outer constructor to define the
    ode_solver = ClimateMachine.ExplicitSolverType(
        solver_method = LSRK144NiegemannDiehlBusch,
    )

    ## The model coefficient for the turbulence closure is defined via the
    ## [CLIMAParameters
    ## package](https://CliMA.github.io/CLIMAParameters.jl/dev/) A reference
    ## state for the linearisation step is also defined.
    T_surface = FT(300)
    T_min_ref = FT(0)
    T_profile = DryAdiabaticProfile{FT}(param_set, T_surface, T_min_ref)
    ref_state = HydrostaticState(T_profile)

    ## The fun part! Here we assemble the `AtmosModel`.
    ##md # !!! note
    ##md #     Docs on model subcomponent options can be found here:
    ##md #     - [`param_set`](https://CliMA.github.io/CLIMAParameters.jl/dev/)
    ##md #     - [`turbulence`](@ref Turbulence-Closures-docs)
    ##md #     - [`source`](@ref atmos-sources)
    ##md #     - [`init_state`](@ref init-dc)

    _C_smag = FT(0.21)
    model = AtmosModel{FT}(
        AtmosLESConfigType,                             # Flow in a box, requires the AtmosLESConfigType
        param_set;                                      # Parameter set corresponding to earth parameters
        init_state_prognostic = init_densitycurrent!,   # Apply the initial condition
        ref_state = ref_state,                          # Reference state
        turbulence = Vreman(_C_smag),                   # Turbulence closure model
        moisture = DryModel(),                          # Exclude moisture variables
        source = (Gravity(),),                          # Gravity is the only source term here
        tracers = NoTracers(),                          # Tracer model with diffusivity coefficients
    )

    ## Finally, we pass a `Problem Name` string, the mesh information, and the
    ## model type to  the [`AtmosLESConfiguration`](@ref ClimateMachine.AtmosLESConfiguration) object.
    config = ClimateMachine.AtmosLESConfiguration(
        "DryDensitycurrent",      # Problem title [String]
        N,                        # Polynomial order [Int]
        resolution,               # (Δx, Δy, Δz) effective resolution [m]
        xmax,                     # Domain maximum size [m]
        ymax,                     # Domain maximum size [m]
        zmax,                     # Domain maximum size [m]
        param_set,                # Parameter set.
        init_densitycurrent!,     # Function specifying initial condition
        solver_type = ode_solver, # Time-integrator type
        model = model,            # Model type
        periodicity = (false, false, false),
        boundary = ((1, 1), (1, 1), (1, 1)),   # Set all boundaries to solid walls
    )
    return config
end

#md # !!! note
#md #     `Keywords` are used to specify some arguments (see appropriate source
#md #     files).

function main()
    ## These are essentially arguments passed to the
    ## [`config_densitycurrent`](@ref config-helper) function.  For type
    ## consistency we explicitly define the problem floating-precision.
    FT = Float64
    ## We need to specify the polynomial order for the DG discretization,
    ## effective resolution, simulation end-time, the domain bounds, and the
    ## courant-number for the time-integrator. Note how the time-integration
    ## components `solver_config` are distinct from the spatial / model
    ## components in `driver_config`. `init_on_cpu` is a helper keyword argument
    ## that forces problem initialisation on CPU (thereby allowing the use of
    ## random seeds, spline interpolants and other special functions at the
    ## initialisation step.)
    N = 4
    Δx = FT(100)
    Δy = FT(250)
    Δv = FT(100)
    resolution = (Δx, Δy, Δv)
    xmax = FT(25600)
    ymax = FT(1000)
    zmax = FT(6400)
    t0 = FT(0)
    timeend = FT(100)
    CFL = FT(1.5)

    ## Assign configurations so they can be passed to the `invoke!` function
    driver_config = config_densitycurrent(FT, N, resolution, xmax, ymax, zmax)
    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        init_on_cpu = true,
        Courant_number = CFL,
    )

    ## Invoke solver (calls `solve!` function for time-integrator), pass the driver, solver and diagnostic config
    ## information.
    result =
        ClimateMachine.invoke!(solver_config; check_euclidean_distance = true)

    ## Check that the solution norm is reasonable.
    @test isapprox(result, FT(1); atol = 1.5e-2)
end

# The experiment definition is now complete. Time to run it.
# `julia --project=$CLIMA_HOME tutorials/Atmos/densitycurrent.jl --vtk 1smins`
# to run with VTK output enabled at intervals of 1 simulation minute.
#
# ## References
#
# [1] J. Straka, R. Wilhelmson, L. Wicker, J. Anderson, K. Droegemeier,
# Numerical solution of a nonlinear density current: a benchmark solution and comparisons
# Int. J. Numer. Methods Fluids 17 (1993) 1–22,  https://doi.org/10.1002/fld.1650170103
#
# [2] R. Carpenter, K. Droegemeier, P. Woodward, C. Hane,
# Application of the piecewise parabolic method (PPM) to meteorological modeling,
# Mon. Weather Rev. 118 (1990) 586–612, https://doi.org/10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2

main()

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
pushfirst!(LOAD_PATH, joinpath("..", ".."))

using CUDA
using Oceananigans

using ConvergenceTests.DoublyPeriodicTaylorGreen: setup_and_run

arch = CUDA.has_cuda() ? GPU() : CPU()
Nx = [8, 16, 32, 64, 128]
stop_time = 0.25

# Calculate time step
max_h = 2π / maximum(Nx)
Δt = 0.01 * max_h^2 # satisfy diffusive constraint for finest resolution
Nt = round(Int, stop_time / Δt)
Δt = stop_time / Nt

for N in Nx
    setup_and_run(architecture=arch, Nx=N, Δt=Δt, stop_iteration=Nt)
end

ext: jl


text:
using PyCall
pyimport("d4rl")

function h5toEB(filename)
    # Load the data from the file
    actions = h5read(filename, "actions")
    obs = h5read(filename, "observations")
    rewards = h5read(filename, "rewards")
    dones = Bool.(h5read(filename, "terminals"))
    timeouts = Bool.(h5read(filename, "timeouts"))

    # we are going to remove all of the observations that included a terminal or done signal
    keep = .!(dones .| timeouts)
    keep[end] = false # remove the last one regardless

    next_obs = obs[:, 2:end][:, keep[1:end-1]] # shift the observations and then keep the ones that we want
    # Keep the rest of the elements consistent
    obs = obs[:,keep]
    actions = actions[:,keep]
    rewards = reshape(rewards[keep], 1, :)
    dones = reshape(dones[keep], 1, :)

    ExperienceBuffer(Dict(:s => obs, :sp=>next_obs, :a=>actions, :r=>rewards, :done=>dones))
end

max_stars_repo_path: examples/offline rl/d4rl_utils.jl

max_issues_repo_path: examples/offline rl/d4rl_utils.jl

max_forks_repo_path: examples/offline rl/d4rl_utils.jl


text:
module TestTrap
import Debug
using Compat, Debug, Debug.AST, Debug.Meta
import Debug.Node, Debug.Scope

is_trap(::Node) = false
is_trap(node::ExNode) = isblocknode(parentof(node))
is_trap(::@compat(Union{Event})) = true

ip = 1

function trap(e::Event, s::Scope)
    global ip
#    println(typeof(e), headof(e.node))
    @assert isa(e, answers[ip][1])
    @assert answers[ip][2] === headof(e.node)
    ip += 1
end
function trap(node::ExNode, s::Scope)
    global ip
#    println(headof(node))
    @assert answers[ip] == headof(node)
    ip += 1
end
function output(x)
    global ip
#    println(x)
    @assert answers[ip] == x
    ip += 1
end

macro test_enterleave(ex)
    Debug.code_debug(Debug.Graft.instrument(is_trap, trap, ex))
end

answers = Any[
    :while, (Enter,:while), (Leave,:while),
    :(=),
    :while,(Enter,:while),:+=,:call,1,:+=,:call,2,:+=,:call,3,(Leave,:while),
    :try, (Enter,:try),
        :while, (Enter,:while), (Leave,:while), :call, "try",
    (Leave,:try),
    :for, (Enter,:for), :call, 2, :call, 3, :call, 4, (Leave,:for),
    :call, (Enter,:let), :call, 5, (Leave,:let), 5,
    :call, (Enter,:comprehension),
        :call, 12, :call, 13,
    (Leave,:comprehension), [12,13],
    :function,
    :call, (Enter,:function), :call, 9, :call, (Leave,:function), 81
]

@test_enterleave begin
    while false
    end
    i = 0
    while i < 3
        i += 1
        output(i)
    end

    try
        while false
        end
        output("try")
    end

    for x=2:4
        output(x)
    end

    output(let x=5
        output(x)
        x
    end)

    output([(output(x); x) for x=12:13])

    function f(x)
        output(x)
        x^2
    end

    output(f(9))
end

@assert ip == length(answers) + 1

end # module


text:
module PFB

using REPL.TerminalMenus

function targets(dir::String, dirs::Bool=true)
    options = readdir(dir, join=true)
    if dirs
        options = filter(isdir, options)
    else
        options = filter(options) do opt
            _, ext = splitext(opt)
            ext == ".csv"
        end
    end
    menu = MultiSelectMenu(options)
    return options, request("Select targets:", menu)
end

include("CovPlots.jl")
include("Execs.jl")

end # module

max_stars_repo_name: cyruscyliu/nyx-net-profuzzbench

max_issues_repo_name: cyruscyliu/nyx-net-profuzzbench

max_forks_repo_name: cyruscyliu/nyx-net-profuzzbench


text:
using POMDPIDMModel
using BeliefUpdaters
using Random
using Plots
using Measures
"""
Rules
1) Pay attention that the length of all routes should cover sRange, but not the reverse,
   or will result in error like 'attempt to access 1-element Array{Geometry,1} at index [2]'
2) Curvature should be negative if route extends clockwise
"""

function DPObj(style::Symbol, Rother::UInt16, Rlength::Float64, gauge::Float64, θ::Float64, Stopline::Float64, Smin::Float64)
    # Init  Sts
    EgoCar = CarSt(0.0, 4.0, UInt16(1))
    OtherCar = CarSt(0.0, 3.0, Rother) # OtherCar.r is not used later.
    SsInit = Sts(EgoCar, OtherCar)

    # space range
    sRange = [0.0:0.1:Rlength;]
    vRange = [0.0:1.0:6.0;] # max velocity = 72 km/h
    Aset = (min=-2.0, max=2.0, comfort=1.0)
    # other parameters
    discount_factor = 0.9
    Δt = 1.0

    if style == :Crossroad
        Vgeos = [4.0 4.0 4.0; 4.0 2.0 4.0; 4.0 3.0 4.0; 4.0 4.0 4.0]
        Routes = RoutesGenerator(Vgeos, Rlength, gauge, θ, 0.15, 0.7; style=style)
        rRange = Vector{UInt16}([2, 3, 4]) # potential route index for other car
        return DrivePOMDP(discount_factor, Δt, sRange, vRange, rRange, Aset, Routes, SsInit, Stopline, Smin)
    elseif style == :TJunction
        Vgeos = [4.0 4.0 4.0; 4.0 2.0 4.0; 4.0 3.0 4.0]
        Routes = RoutesGenerator(Vgeos, Rlength, gauge, θ, 0.25, 0.35; style=style)
        rRange = Vector{UInt16}([2, 3])
        return DrivePOMDP(discount_factor, Δt, sRange, vRange, rRange, Aset, Routes, SsInit, Stopline, Smin)
    else
        error("Unsupported Style. Please choose either :Crossroad or :TJunction.")
    end
end

Rlength = 24.0
gauge = 4.0
θ = 110.0/(180/π)
halflen = Rlength/2.0+7.0
style = :TJunction; Stopline = 4.0
#style = :Crossroad; Stopline = 10.0
Smin = 2.0

# init
DP = DPObj(style, UInt16(2), Rlength, gauge, θ, Stopline, Smin)
pltJunction = Junction(halflen, gauge, θ; style=style)
pltMap = PlotRoutes(DP, pltJunction)
DrawStopline!(pltMap, DP.Routes[DP.SsInit.Ego.r], DP.Stopline, gauge)
plot(pltMap, size=(600, 400))
#plot(pltMap, size=(600, 500))
savefig("output/Figure4.6a.pdf")

#=
Svec = [i for i in 0.0:0.1:Rlength]
idx_r12 = findfirst(x->x==3.0, DP.Routes[3].Vref)+1
Svec_r1 = Svec[1:idx_r12]
Vref_r1 = DP.Routes[3].Vref[1:idx_r12]
Vref_r23 = DP.Routes[3].Vref[idx_r12:end]
idx_r23 = findlast(x->x==3.0, Vref_r23)+idx_r12-1
Svec_r2 = Svec[idx_r12:idx_r23]
Vref_r2 = DP.Routes[3].Vref[idx_r12:idx_r23]
Svec_r3 = Svec[idx_r23:end]
Vref_r3 = DP.Routes[3].Vref[idx_r23:end]

Aref_r1 = DP.Routes[3].Aref[1:idx_r12]
Aref_r2 = DP.Routes[3].Aref[idx_r12:idx_r23]
Aref_r3 = DP.Routes[3].Aref[idx_r23:end]

#=
plt1 = plot(Svec_r1, Vref_r1, xlim=(0.0, Rlength+1), ylim=(0.0, 5.0), c=:green, size=(460, 310), grid=true, legend=false, ylabel="V(m/s)", left_margin=-2mm, titleposition=:left, showaxis=true, xlabel="S(m)", xticks =0:3.0:Rlength+1, yticks =0:1.0:5.0)
plot!(plt1, Svec_r2, Vref_r2, c=:red)
plot!(plt1, Svec_r3, Vref_r3, c=:blue)
savefig("output/Vref.pdf")
=#
plt2 = plot(Svec_r1, Aref_r1, xlim=(0.0, Rlength+1), ylim=(-2.0, 4.0), c=:green, size=(460, 310), grid=true, legend=false, ylabel="A(m/s²)", left_margin=-2mm, titleposition=:left, showaxis=true, xlabel="S(m)", xticks =0:3.0:Rlength+1, yticks =-2.0:1.0:4.0)
plot!(plt2, Svec_r2, Aref_r2, c=:red)
plot!(plt2, Svec_r3, Aref_r3, c=:blue)
savefig("output/Aref.pdf")
=#


text:
#-------------------------------------------------------------------------------
#conversions between spectral units

export ν2f, f2ν, ν2λ, λ2ν, λ2f, f2λ

"""
Convert wavenumber [cm``^{-1}``] to frequency [1/s]
"""
ν2f(ν) = 100.0*𝐜*ν

"""
Convert frequency [1/s] to wavenumber [cm``^{-1}``]
"""
f2ν(f) = f/(100.0*𝐜)

"""
Convert wavenumber [cm``^{-1}``] to wavelength [m]
"""
ν2λ(ν) = 0.01/ν

"""
Convert wavelength [m] to wavenumber [cm``^{-1}``]
"""
λ2ν(λ) = 0.01/λ

"""
Convert wavelength [m] to frequency [1/s]
"""
λ2f(λ) = 𝐜/λ

"""
Convert frequency [1/s] to wavelength [m]
"""
f2λ(f) = f/𝐜

#-------------------------------------------------------------------------------
export planck, normplanck, dplanck, stefanboltzmann, equilibriumtemperature

"""
    planck(ν, T)

Compute black body intensity [W/m``^2``/cm``^{-1}``/sr] using [Planck's law](https://en.wikipedia.org/wiki/Planck%27s_law)

# Arguments
* `ν`: wavenumger [cm``^{-1}``]
* `T`: temperature [Kelvin]
"""
function planck(ν, T)
    νₘ = 100.0*ν #convert from cm^-1 to m^-1
    x = 𝐡*𝐜*νₘ/(𝐤*T) #exponent
    p = 2*𝐡*𝐜^2*νₘ^3 #prefactor
    #result, converting back to cm^-1 with the factor of 100
    100.0*p/(exp(x) - 1.0)
end

"""
    normplanck(ν, T)

Compute black body intensity [W/m``^2``/cm``^{-1}``/sr] using [Planck's law](https://en.wikipedia.org/wiki/Planck%27s_law), normalized by the power emitted per unit area at the given temperature ([`stefanboltzmann`](@ref)),

```
B(ν,T)/σT^4
```

yielding units of 1/cm``^{-1}``/sr.

# Arguments
* `ν`: wavenumger [cm``^{-1}``]
* `T`: temperature [Kelvin]
"""
normplanck(ν, T) = planck(ν, T)/stefanboltzmann(T)

"""
    dplanck(ν, T)

Evaluates ``\\frac{\\partial B}{\\partial T}``
"""
function dplanck(ν, T)
    νₘ = 100.0*ν #convert from cm^-1 to m^-1
    x = 𝐡*𝐜*νₘ/(𝐤*T) #exponent
    y = exp(x)
    if isinf(y)
        return zero(y) #fail accurately, Inf/(Inf - 1)^2 → 0
    end
    p = 2*𝐡^2*𝐜^3*νₘ^4/(𝐤*T^2) #prefactor
    #result, converting back to cm^-1 with the factor of 100
    100.0*p*y/(y - 1)^2
end

"""
    stefanboltzmann(T)

Compute black body radiation power using the [Stefan-Boltzmann](https://en.wikipedia.org/wiki/Stefan%E2%80%93Boltzmann_law) law, ``σT^4`` [W/m``^2``].
"""
stefanboltzmann(T) = 𝛔*(T^4)

"""
    equilibriumtemperature(F, A)

Compute the [planetary equilibrium temperature](https://en.wikipedia.org/wiki/Planetary_equilibrium_temperature), or equivalent blackbody temperature of a planet.

``(\\frac{(1 - A)F}{4\\sigma})^{1/4}``

# Arguments
* `F`: stellar flux [W/m``^2``]
* `A`: albedo
"""
equilibriumtemperature(F, A) = ((1 - A)*F/(4*𝛔))^(1/4)

"""
    equilibriumtemperature(L, A, R)

Compute the [planetary equilibrium temperature](https://en.wikipedia.org/wiki/Planetary_equilibrium_temperature), or equivalent blackbody temperature of a planet.

``(\\frac{(1 - A)L}{16 \\sigma \\pi R^2})^{1/4}``

# Arguments
* `L`: stellar luminosity [W]
* `A`: albedo
* `R`: orbital distance [m]
"""
equilibriumtemperature(L, A, R) = (L*(1 - A)/(16*𝛔*π*R^2))^(1/4)

#-------------------------------------------------------------------------------
export dτdP, transmittance, schwarzschild, absorption, emission

"""
    dτdP(σ, g, μ)

Evaluate the differential increase in optical depth in pressure coordinates, equivalent to the [`schwarzschild`](@ref) equation without Planck emission.

``\\frac{dτ}{dP} = σ\\frac{\\textrm{N}_A}{g μ}``

where ``N_A`` is Avogadro's number.

# Arguments
* `σ`: absorption cross-section [cm``^2``/molecule]
* `g`: gravitational acceleration [m/s``^2``]
* `μ`: mean molar mass [kg/mole]
"""
dτdP(σ, g, μ) = 1e-4*σ*𝐍𝐚/(μ*g)

"""
    transmittance(τ)

Evaluate transmittance from optical depth, ``t = e^{-τ}``
"""
transmittance(τ) = exp(-τ)

"""
    schwarzschild(I, ν, σ, T, P)

Evaluate the [Schwarzschild differential equation](https://en.wikipedia.org/wiki/Schwarzschild%27s_equation_for_radiative_transfer) for radiative transfer with units of length/height [m] and assuming the ideal gas law.

``\\frac{dI}{dz} = σ\\frac{P}{k_B T}[B_ν(T) - I]``

where ``B_ν`` is [`planck`](@ref)'s law.

# Arguments
* `I`: radiative intensity [W/m``^2``/cm``^{-1}``/sr]
* `ν`: radiation wavenumber [cm``^{-1}``]
* `σ`: absorption cross-section [cm``^2``/molecule]
* `T`: temperature [K]
* `P`: pressure [Pa]
"""
schwarzschild(I, ν, σ, T, P) = 1e-4*σ*(P/(𝐤*T))*(planck(ν,T) - I)

"""
    schwarzschild(I, ν, σ, g, μ, T)

Evaluate the [Schwarzschild differential equation](https://en.wikipedia.org/wiki/Schwarzschild%27s_equation_for_radiative_transfer) for radiative transfer with pressure units [Pa] and assuming the ideal gas law.

``\\frac{dI}{dP} = σ\\frac{\\textrm{N}_A}{g μ}[B_ν(T) - I]``

where ``B_ν`` is [`planck`](@ref)'s law and ``N_A`` is Avogadro's number.

# Arguments
* `I`: radiative intensity [W/m``^2``/cm``^{-1}``/sr]
* `ν`: radiation wavenumber [cm``^{-1}``]
* `σ`: absorption cross-section [cm``^2``/molecule]
* `g`: gravitational acceleration [m/s``^2``]
* `μ`: mean molar mass [kg/mole]
* `T`: temperature [K]
"""
schwarzschild(I, ν, σ, g, μ, T) = 1e-4*σ*(𝐍𝐚/(μ*g))*(planck(ν,T) - I)

"""
    absorption(I, σ, g, μ)

Evaluate the [Schwarzschild differential equation](https://en.wikipedia.org/wiki/Schwarzschild%27s_equation_for_radiative_transfer) for radiative transfer with pressure units [Pa] and assuming the ideal gas law, but **without Planck emission**. This can be useful if the wavenumber is far away from the emission region of the temperature and the Planck emission is known to be negligible.

``\\frac{dI}{dP} = -σ\\frac{\\textrm{N}_A}{g μ}I``

where ``N_A`` is Avogadro's number.

# Arguments
* `I`: radiative intensity [W/m``^2``/cm``^{-1}``/sr]
* `σ`: absorption cross-section [cm``^2``/molecule]
* `g`: gravitational acceleration [m/s``^2``]
* `μ`: mean molar mass [kg/mole]
"""
absorption(I, σ, g, μ) = -1e-4*σ*(𝐍𝐚/(μ*g))*I

"""
    emission(I, ν, σ, g, μ, T)

Evaluate the [Schwarzschild differential equation](https://en.wikipedia.org/wiki/Schwarzschild%27s_equation_for_radiative_transfer) for radiative transfer with pressure units [Pa] and assuming the ideal gas law, but **without Planck emission**. This can be useful if the wavenumber is far away from the emission region of the temperature and the Planck emission is known to be negligible.

``\\frac{dI}{dP} = -σ\\frac{\\textrm{N}_A}{g μ}B_ν(T)``

where ``B_ν`` is [`planck`](@ref)'s law

# Arguments
* `I`: radiative intensity [W/m``^2``/cm``^{-1}``/sr]
* `σ`: absorption cross-section [cm``^2``/molecule]
* `g`: gravitational acceleration [m/s``^2``]
* `μ`: mean molar mass [kg/mole]
"""
emission(ν, σ, g, μ, T) = 1e-4*σ*(𝐍𝐚/(μ*g))*planck(ν,T)


text:
using Knet
using Images
using ImageDraw

add_dim(x::Array) = reshape(x, (size(x)..., 1));

reduce_dim(x::Array) = reshape(x, (size(x)[1:end-1]));

# Not working as expected
function concatenate_iterables(iterables...)
    return zip(iterables...) |> Iterators.flatten
end

# mutable struct Data{T}; x; y; batchsize; length; partial; imax; indices; shuffle; xsize; ysize; xtype; ytype; end
function append_to_data!(data::Knet.Train20.Data, another_data::Knet.Train20.Data)
    println(size(data.x))
    data.x = cat(data.x, another_data.x, dims = 2)
    data.y = cat(data.y, another_data.y, dims = 2)
    data.length = data.length + another_data.length
    data.imax = data.length - data.batchsize + 1
    data.indices = 1:data.length
end

function model_end_results(
    model,
    train_data,
    validation_data;
    train_data_items,
    validation_data_items,
)
    train_loss = model(train_data)
    validation_loss = model(validation_data)
    train_acc = modelized_naive_pck_sigm(model, train_data)
    validation_acc = modelized_naive_pck_sigm(model, validation_data)

    train_acc_pckh = nothing
    validation_acc_pckh = nothing

    if train_data_items != nothing && validation_data_items != nothing
        train_acc_pckh = modelized_PCKh_sigm(model, train_data, train_data_items)[1]
        validation_acc_pckh =
            modelized_PCKh_sigm(model, validation_data, validation_data_items)[1]
    end

    end_res_dict = Dict(
        "train_loss" => train_loss,
        "validation_loss" => validation_loss,
        "train_acc" => train_acc,
        "validation_acc" => validation_acc,
        "train_acc_PCKh" => train_acc_pckh,
        "validation_PCKh" => validation_acc_pckh,
    )
    for (k, v) in end_res_dict
        println("$(k): $(v)")
    end
end

function display_deepercut_input_image(input_image)
    colorviews = []

    function single_display(single_img)
        perm = permutedims(single_img, [3, 1, 2])
        push!(colorviews, colorview(RGB, perm))
    end

    input_image_sizes = size(input_image)
    if input_image_sizes |> length == 4
        for i = 1:input_image_sizes[end]
            img = input_image[:, :, :, i]
            single_display(img)
        end
    else
        single_display(input_image)
    end

    return colorviews
end

#=
# TODO: This can be used for visualization purposes
arrange_img_damacfix = convert(Array{FixedPointNumbers.Normed{UInt8,8},3}, channelview(loaded_img));
converted_img = convert(Array{Float32}, arrange_img_damacfix);
img_v = @view converted_img[:, 38:75 , 112:177] # blue;
fill!(img_v, 0) # b;
=#

function fetch_loc_reffed_max_point(scmap, joint_id, add_loc_ref)
  offmat = scmap[:, :, global_num_joints+1:global_num_joints*3]
            max_coord = argmax(scmap[:,:, joint_id])
            pred = (max_coord[1], max_coord[2])
            pred_f8 = pred .* preprocess_stride .+ (0.5 * preprocess_stride)
            # TODO: ORIGINAL Be careful about here
            dx = offmat[:, :, 2*joint_id-1][max_coord] * global_locref_stdev
            dy = offmat[:, :, 2*joint_id][max_coord] * global_locref_stdev

            # println("***")
            #println(pred_f8)
            if add_loc_ref
                 pred_f8 = (pred_f8[1] + dx, pred_f8[2] + dy)
            end

          #  println(" after $(pred_f8)")
           # println("dx: $(dx) , dy: $(dy)")
            pred_f8 = CartesianIndex(pred_f8[1] |> round |> Int , pred_f8[2] |> round |> Int)
    return pred_f8
end

function show_scmap_on_image(
    image,
    scmap;
    should_display = false,
    should_use_scmap_size = true,
    display_name = nothing,
        confidence_threshold = 0,
        focus_on_argmax = false,
        add_loc_ref_offset = false,
        return_single_image = false,
        should_color_scmap = true,
        use_sigm = true
)
    if use_sigm
        scmap = sigm.(scmap)
    end

    colored_images = []
    #size_h = should_use_scmap_size ? size(scmap)[1] : size(image)[1]
   # size_w = should_use_scmap_size ? size(scmap)[2] : size(image)[2]
  #  output_sized_image = imresize(image, size_h, size_w)
      size_h = size(image)[1]
    size_w = size(image)[2]
   output_sized_image = image
    perm = permutedims(output_sized_image, [3, 1, 2])

    if should_color_scmap

      all_joints = [[1, 6], [2, 5], [3, 4], [7, 12], [8, 11], [9, 10], [13], [14]]
        for (joint_group) in all_joints
            for joint_id in joint_group

        scmap_slice = Array(scmap[:,:, joint_id])
        # scmap_slice = imresize(scmap_slice, size_h, size_w)

        if focus_on_argmax != true
            scmap_slice = imresize(scmap_slice, size(image)[1],  size(image)[2])
            color_idx = findall(scmap_slice .> confidence_threshold)
                else
         #=
        elseif add_loc_ref_offset # && !should_use_scmap_size
           pred_f8 =  fetch_loc_reffed_max_point(scmap, joint_id)
            color_idx = create_indexes_around_center(pred_f8, size_w)
        else
            max_point = argmax(scmap_slice)
            color_idx = create_indexes_around_center(max_point, max(size_h, size_w))
        end
                =#
        pred_f8 =  fetch_loc_reffed_max_point(scmap, joint_id, add_loc_ref_offset)
        color_idx = create_indexes_around_center(pred_f8,  max(size_h, size_w))
                end
        img_r = @view perm[1, color_idx]
        img_g = @view perm[2, color_idx]
        img_b = @view perm[3, color_idx]

        fill!(img_r, (joint_id + 4) / (global_num_joints + 4) * 1.2)
        fill!(img_b, 1.4 - ( (joint_id + 4) / (global_num_joints + 4) * 1.4))
        fill!(img_g, (joint_id + 4) / (global_num_joints + 4) * 1)

        #fill!(img_r, 1)

            if !return_single_image
                colored = colorview(RGB, perm)
                push!(colored_images, colored)
                perm = permutedims(output_sized_image, [3, 1, 2])
            end

    end
    end
    end

   return return_single_image ? colorview(RGB, perm) : colored_images
end

function create_indexes_around_center(center_point, array_max_dim_size; distance = 5)
    initial_x = center_point[1]
    initial_y = center_point[2]

    top_x = initial_x + distance
    bottom_x = initial_x - distance
    top_y = initial_y + distance
    bottom_y = initial_y - distance

    p1 = (top_x, top_y)
    p2 = (top_x, bottom_y)
    p3 = (bottom_x, top_y)
    p4 = (bottom_x, bottom_y)

    candidates = [p1, p2, p3, p4]

    filtered_candidates = filter(e -> !(true in (e .< 1)) && !(true in (e .> array_max_dim_size)) , candidates)

    if filtered_candidates |> length == 4
        return CartesianIndices((bottom_x:top_x, bottom_y:top_y))
    else
        return CartesianIndices((initial_x:initial_x, initial_y:initial_y))
    end
end

function draw_score_maps(y)
    images = []

    function draw_score_maps_for_single(y_single)
        #= all_joints = Dict(
            "ankle" => [1, 6],
            "knee" => [2, 5],
            "hip" => [3, 4],
           "wrist" => [7, 12],
            "elbow" => [8, 11],
            "shoulder" => [9, 10],
            "chin" => [13],
         "forehead" => [14]
        ) =#
        images_for_single = []
        all_joints = [[1, 6], [2, 5], [3, 4], [7, 12], [8, 11], [9, 10], [13], [14]]
        for (joint_group) in all_joints
            for j in joint_group
                permuted = permutedims(y_single[:, :, j], (1, 2))
                push!(images_for_single, colorview(Gray, permuted))
            end
        end
        mosaic = mosaicview(
            images_for_single...,
            ncol = 2,
            rowmajor = true,
            npad = 1,
            fillvalue = 0.3,
        )
        display(mosaic)
        push!(images, images_for_single...)
    end

    y_sizes = size(y)
    if y_sizes |> length == 4
        for i = 1:y_sizes[end]
            single_y = y[:, :, :, i]
            draw_score_maps_for_single(single_y)
        end
    else
        draw_score_maps_for_single(y)
    end

    return images
end

function mosaicify(images; ncol = 2, rowmajor = true, npad = 1, fillvalue = 0.3)
    mosaic = mosaicview(
        images...,
        ncol = ncol,
        rowmajor = rowmajor,
        npad = npad,
        fillvalue = fillvalue,
    )
    display(mosaic)
end

function get_object_tag(obj)
    obj_field_names = fieldnames(typeof(obj))
    has_tag = :tag in obj_field_names
    if has_tag
        return obj.tag
    else
        return nothing
    end
end

function read_accuracy_results(acc_dist_dict)
    person_types = ["single", "multi", "mixed"]
    h_ranges = [hr for hr = 0:0.1:1]
    joint_ids = [j for j = 1:global_num_joints]
    acc_res = zeros(3 * 11, global_num_joints)
    for (h_i, h_range) in enumerate(h_ranges)
        for joint_id in joint_ids
            for (p_i, person_type) in enumerate(person_types)
                accurate = acc_dist_dict[person_type][h_range][joint_id][1]
                total = acc_dist_dict[person_type][h_range][joint_id][2]
                perc = 100.0 * accurate / total
                acc_res[(h_i-1)*3+p_i, joint_id] = round(perc, digits = 2)
            end
        end
    end
    return acc_res
end

function write_results_to_file(filename, results_kv_tuples...)
    open(filename, "a") do io
        for (k, v) in results_kv_tuples
            write(io, "$(k): $(v)")
            write(io, "\n")
        end
    end
end


text:
struct MixedUpdater <: Updater
    problem::SingleOCPOMDP
    dt::Float64
    states::Vector{SingleOCState}
end

function MixedUpdater(problem::SingleOCPOMDP, dt::Float64 = 0.1)
    return MixedUpdater(problem, dt, states(problem))
end

# efficient version of the discrete updater
function POMDPs.update(bu::MixedUpdater, bold::SingleOCBelief, a::SingleOCAction, o::SingleOCObs)
    bnew = SingleOCBelief()
    pomdp = bu.problem

    ego = o.ego
    Y = AutomotivePOMDPs.get_Y_grid(pomdp)
    V_ped = AutomotivePOMDPs.get_V_ped_grid(pomdp) # try to avoid collect
    n_ped_states = length(Y)*length(V_ped) + 1
    saved_prob = 0.
    for i=1:n_ped_states
        if i == n_ped_states
            ped = get_off_the_grid(pomdp)
        else
            k, l = Tuple(CartesianIndices((length(Y), length(V_ped)))[i])
            y = Y[k]
            v = V_ped[l]
            ped = AutomotivePOMDPs.yv_to_state(pomdp, y, v)
        end
        collision = collision_checker(ego, ped, pomdp.ego_type, pomdp.ped_type)
        sp = SingleOCState(collision, ego, ped)
        probo = obs_weight(o, sp, pomdp)
        # if observation prob is 0.0, then skip rest of update b/c bnew[i] is zero
        if probo ≈ 0.
            continue
        elseif saved_prob < probo
            saved_prob = probo
        end

        b_sum = 0.0 # belief for state sp
        for (j, s) in enumerate(bold.it)
            td = transition(pomdp, s, a, bu.dt)
            # use pdf on the pedestrian only
            pp = partial_pdf(td, sp)
            b_sum += pp * bold.p[j]
        end
        if b_sum != 0.
            push!(bnew.it, sp)
            push!(bnew.p, probo*b_sum)
        end
    end
    # if norm is zero, the update was invalid
    if sum(bnew.p) == 0.0
        # println("Invalid update for: ", bold, " ", a, " ", o)
        # val, ind = fargmax(bold.p)
        # ml_state = bold.it[ind]
        println("Invalid update for: observation ", o,
                "\n action ", a,
                "\n prob", saved_prob,
                "\n old belief ", bold)
        # bnew = SingleOCDistribution([1.0], [o])
        throw("UpdateError")
        # println("Invalid update, probo = $probo")
        # bnew = observation(pomdp, a, o)
    else
        normalize!(bnew.p, 1)
    end
    bnew
end

function partial_pdf(d::SingleOCDistribution, sp::SingleOCState)
    for (i, s) in enumerate(d.it)
        if sp.ped == s.ped
            return d.p[i]
        end
    end
    return 0.
end

"""
    initialize_off_belief(pomdp::OCPOMDP)
Initialize a belief on all the non observable position
"""
function initialize_off_belief(policy::QMDPEval, ego::VehicleState)
    b = Dict{Int64, SingleOCDistribution}()
    b[OFF_KEY] = initialstate(policy.pomdp, ego)
    return b
end

function POMDPs.initialstate(pomdp::OCPOMDP, n::Int64)
    init_dis = DecBelief()
    for i=1:n
        init_dis[i+1] = initialstate(pomdp)
    end
    return init_dis
end

function POMDPs.update(up::MixedUpdater, bold::Dict{Int64, SingleOCDistribution}, a::Float64, o::Dict{Int64, SingleOCObs})
    bnew = Dict{Int64, SingleOCDistribution}()
    for oid in keys(o)
        if haskey(bold, oid) && oid != OFF_KEY
            bnew[oid] = update(up, bold[oid], SingleOCAction(a), o[oid])
        elseif oid == OFF_KEY
            bnew[OFF_KEY] = update(up, bold[OFF_KEY], SingleOCAction(a), o[oid])
        else # cars appeared
            bnew[oid] = update(up, bold[OFF_KEY], SingleOCAction(a), o[oid])
        end
    end
    # print_b_stats(up.problem, bnew[OFF_KEY])
    return bnew
end

function print_b_stats(pomdp::SingleOCPOMDP, b::SingleOCDistribution)
    off_prob = 0.
    for (i,s) in enumerate(b.it)
        if off_the_grid(pomdp, s.ped)
            off_prob = b.p[i]
        end
    end
    max_s = most_likely_state(b)
    println("States to consider: $(length(b.it)) ; off prob $off_prob, Most likely state $max_s ")
end


text:
using Statistics

# Input data
data = map(line -> map(p -> p[1], split(line, "")), readlines("Day10.input"));

# Part 1
scoreparenthesis =
    (c) -> Int.(c .== ('(', '[', '{', '<')) .- Int.(c .== (')', ']', '}', '>'));

tocallstack =
    (f) -> reduce(
        (acc, v) -> (
            any(scoreparenthesis(v) .< 0) ?
            length(acc) == 0 ? push!(acc, (v, false)) :
            (
                last(acc)[2] &&
                all((scoreparenthesis(v) .+ scoreparenthesis(pop!(acc)[1])) .== 0) ?
                acc : push!(acc, (v, false))
            ) : push!(acc, (v, true))
        ),
        f,
        init = [],
    );

ans = sum(
    reduce(
        (acc, v) -> acc .+ (v .* (3, 57, 1197, 25137)),
        map(
            v -> scoreparenthesis(first(v)[1]) .* -1,
            filter(!isempty, map(d -> filter(t -> !t[2], tocallstack(d)), data)),
        ),
        init = (0, 0, 0, 0),
    ),
);
println("Part 1: $ans")

# Part 2
ans = Int(Statistics.median(
    map(
        v -> sum(
            reduce(
                (acc, p) -> (acc .* 5) .+ (scoreparenthesis(p[1]) .* (1, 2, 3, 4)),
                reverse(v),
                init = (0, 0, 0, 0),
            ),
        ),
        filter(s -> all(map(p -> p[2], s)), map(v -> tocallstack(v), data)),
    ),
));
println("Part 2: $ans")


text:
precedence(op::Int) = op < Tokens.end_assignments ?  AssignmentOp :
                       op < Tokens.end_pairarrow ? 2 :
                       op < Tokens.end_conditional ? ConditionalOp :
                       op < Tokens.end_arrow ?       ArrowOp :
                       op < Tokens.end_lazyor ?      LazyOrOp :
                       op < Tokens.end_lazyand ?     LazyAndOp :
                       op < Tokens.end_comparison ?  ComparisonOp :
                       op < Tokens.end_pipe ?        PipeOp :
                       op < Tokens.end_colon ?       ColonOp :
                       op < Tokens.end_plus ?        PlusOp :
                       op < Tokens.end_bitshifts ?   BitShiftOp :
                       op < Tokens.end_times ?       TimesOp :
                       op < Tokens.end_rational ?    RationalOp :
                       op < Tokens.end_power ?       PowerOp :
                       op < Tokens.end_decl ?        DeclarationOp :
                       op < Tokens.end_where ?       WhereOp : DotOp

precedence(kind::Tokens.Kind) = kind == Tokens.DDDOT ? DddotOp :
                        kind < Tokens.begin_assignments ? 0 :
                        kind < Tokens.end_assignments ?   AssignmentOp :
                        kind < Tokens.end_pairarrow ?   2 :
                       kind < Tokens.end_conditional ?    ConditionalOp :
                       kind < Tokens.end_arrow ?          ArrowOp :
                       kind < Tokens.end_lazyor ?         LazyOrOp :
                       kind < Tokens.end_lazyand ?        LazyAndOp :
                       kind < Tokens.end_comparison ?     ComparisonOp :
                       kind < Tokens.end_pipe ?           PipeOp :
                       kind < Tokens.end_colon ?          ColonOp :
                       kind < Tokens.end_plus ?           PlusOp :
                       kind < Tokens.end_bitshifts ?      BitShiftOp :
                       kind < Tokens.end_times ?          TimesOp :
                       kind < Tokens.end_rational ?       RationalOp :
                       kind < Tokens.end_power ?          PowerOp :
                       kind < Tokens.end_decl ?           DeclarationOp :
                       kind < Tokens.end_where ?          WhereOp :
                       kind < Tokens.end_dot ?            DotOp :
                       kind == Tokens.ANON_FUNC ? AnonFuncOp :
                       kind == Tokens.PRIME ?             PrimeOp : 20

precedence(x) = 0
precedence(x::AbstractToken) = precedence(x.kind)
precedence(x::OPERATOR) = precedence(x.kind)

isunaryop(op) = false
isunaryop(op::OPERATOR) = isunaryop(op.kind)
isunaryop(t::AbstractToken) = isunaryop(t.kind)
isunaryop(kind::Tokens.Kind) = kind == Tokens.ISSUBTYPE ||
                  kind == Tokens.ISSUPERTYPE ||
                  kind == Tokens.PLUS ||
                  kind == Tokens.MINUS ||
                  kind == Tokens.NOT ||
                  kind == Tokens.APPROX ||
                  kind == Tokens.NOT_SIGN ||
                  kind == Tokens.AND ||
                  kind == Tokens.SQUARE_ROOT ||
                  kind == Tokens.CUBE_ROOT ||
                  kind == Tokens.QUAD_ROOT ||
                  kind == Tokens.DECLARATION ||
                  kind == Tokens.EX_OR ||
                  kind == Tokens.COLON

isunaryandbinaryop(t) = false
isunaryandbinaryop(t::AbstractToken) = isunaryandbinaryop(t.kind)
isunaryandbinaryop(kind::Tokens.Kind) = kind == Tokens.PLUS ||
                           kind == Tokens.MINUS ||
                           kind == Tokens.EX_OR ||
                           kind == Tokens.ISSUBTYPE ||
                           kind == Tokens.ISSUPERTYPE ||
                           kind == Tokens.AND ||
                           kind == Tokens.APPROX ||
                           kind == Tokens.DECLARATION ||
                           kind == Tokens.COLON

isbinaryop(op) = false
isbinaryop(op::OPERATOR) = isbinaryop(op.kind)
isbinaryop(t::AbstractToken) = isbinaryop(t.kind)
isbinaryop(kind::Tokens.Kind) = isoperator(kind) &&
                    !(kind == Tokens.SQUARE_ROOT ||
                    kind == Tokens.CUBE_ROOT ||
                    kind == Tokens.QUAD_ROOT ||
                    kind == Tokens.NOT ||
                    kind == Tokens.NOT_SIGN)

isassignment(t::AbstractToken) = Tokens.begin_assignments < t.kind < Tokens.end_assignments

function non_dotted_op(t::AbstractToken)
    k = t.kind
    return (k == Tokens.COLON_EQ ||
            k == Tokens.PAIR_ARROW ||
            k == Tokens.EX_OR_EQ ||
            k == Tokens.CONDITIONAL ||
            k == Tokens.LAZY_OR ||
            k == Tokens.LAZY_AND ||
            k == Tokens.ISSUBTYPE ||
            k == Tokens.ISSUPERTYPE ||
            k == Tokens.LPIPE ||
            k == Tokens.RPIPE ||
            k == Tokens.EX_OR ||
            k == Tokens.COLON ||
            k == Tokens.DECLARATION ||
            k == Tokens.IN ||
            k == Tokens.ISA ||
            k == Tokens.WHERE ||
            (isunaryop(k) && !isbinaryop(k) && !(k == Tokens.NOT)))
end

issyntaxcall(op) = false
function issyntaxcall(op::OPERATOR)
    K = op.kind
    P = precedence(K)
    P == AssignmentOp && !(K == Tokens.APPROX || K == Tokens.PAIR_ARROW) ||
    K == Tokens.RIGHT_ARROW ||
    P == LazyOrOp ||
    P == LazyAndOp ||
    K == Tokens.ISSUBTYPE ||
    K == Tokens.ISSUPERTYPE ||
    K == Tokens.COLON ||
    K == Tokens.DECLARATION ||
    K == Tokens.DOT ||
    K == Tokens.DDDOT ||
    K == Tokens.PRIME ||
    K == Tokens.WHERE
end

issyntaxunarycall(op) = false
function issyntaxunarycall(op::OPERATOR)
    K = op.kind
    !op.dot && (K == Tokens.EX_OR ||
    K == Tokens.AND ||
    K == Tokens.DECLARATION ||
    K == Tokens.ISSUBTYPE ||
    K == Tokens.ISSUPERTYPE)
end

LtoR(prec::Int) = AssignmentOp ≤ prec ≤ LazyAndOp || prec == PowerOp

"""
    parse_unary(ps)

Having hit a unary operator at the start of an expression return a call.
"""
function parse_unary(ps::ParseState, op::OPERATOR)
    K,dot = op.kind, op.dot
    if op isa OPERATOR && op.kind == Tokens.COLON
        ret = parse_unary_colon(ps, op)
    elseif (is_plus(op) || is_minus(op)) && (ps.nt.kind == Tokens.INTEGER || ps.nt.kind == Tokens.FLOAT) && isemptyws(ps.ws) && ps.nnt.kind!=Tokens.CIRCUMFLEX_ACCENT
        arg = LITERAL(next(ps))
        ret = LITERAL(op.fullspan + arg.fullspan, (op.fullspan + arg.span), string(is_plus(op) ? "+" : "-" , val(ps.t, ps)), ps.t.kind)
    else
        P = precedence(K)
        prec = P == DeclarationOp ? DeclarationOp :
                    K == Tokens.AND ? DeclarationOp :
                    K == Tokens.EX_OR ? 20 : PowerOp
        arg = @closer ps unary @precedence ps prec parse_expression(ps)
        if issyntaxunarycall(op)
            ret = UnarySyntaxOpCall(op, arg)
        else
            ret = UnaryOpCall(op, arg)
        end
    end

    return ret
end

function parse_unary_colon(ps::ParseState, op::OPERATOR)
    if Tokens.begin_keywords < ps.nt.kind < Tokens.end_keywords
        ret = EXPR{Quotenode}(Any[op, IDENTIFIER(next(ps))])
    elseif Tokens.begin_literal < ps.nt.kind < Tokens.end_literal ||
        isoperator(ps.nt.kind) || ps.nt.kind == Tokens.IDENTIFIER
        ret = EXPR{Quotenode}(Any[op, INSTANCE(next(ps))])
    elseif closer(ps)
        ret = op
    else
        arg = @precedence ps 20 parse_expression(ps)
        ret = EXPR{Quote}(Any[op, arg])
    end
    return ret
end

function parse_operator_eq(ps::ParseState, @nospecialize(ret), op)
    nextarg = @precedence ps AssignmentOp - LtoR(AssignmentOp) parse_expression(ps)

    if is_func_call(ret) && !(nextarg isa EXPR{Begin} || (nextarg isa EXPR{InvisBrackets} && nextarg.args[2] isa EXPR{Block}))
        nextarg = EXPR{Block}(Any[nextarg])
    end
    return BinarySyntaxOpCall(ret, op, nextarg)
end

# Parse conditionals
function parse_operator_cond(ps::ParseState, @nospecialize(ret), op)
    nextarg = @closer ps ifop parse_expression(ps)
    op2 = OPERATOR(next(ps))
    nextarg2 = @closer ps comma @precedence ps 0 parse_expression(ps)

    return ConditionalOpCall(ret, op, nextarg, op2, nextarg2)
end

# Parse comparisons
function parse_comp_operator(ps::ParseState, @nospecialize(ret), op)
    nextarg = @precedence ps ComparisonOp - LtoR(ComparisonOp) parse_expression(ps)

    if ret isa EXPR{Comparison}
        push!(ret, op)
        push!(ret, nextarg)
    elseif ret isa BinaryOpCall && precedence(ret.op) == ComparisonOp
        ret = EXPR{Comparison}(Any[ret.arg1, ret.op, ret.arg2, op, nextarg])
    elseif ret isa BinarySyntaxOpCall && (is_issubt(ret.op) || is_issupt(ret.op))
        ret = EXPR{Comparison}(Any[ret.arg1, ret.op, ret.arg2, op, nextarg])
    elseif (is_issubt(op) || is_issupt(op))
        ret = BinarySyntaxOpCall(ret, op, nextarg)
    else
        ret = BinaryOpCall(ret, op, nextarg)
    end
    return ret
end

# Parse ranges
function parse_operator_colon(ps::ParseState, @nospecialize(ret), op)
    nextarg = @precedence ps ColonOp - LtoR(ColonOp) parse_expression(ps)

    if ret isa BinaryOpCall && is_colon(ret.op)
        ret = EXPR{ColonOpCall}(Any[ret.arg1, ret.op, ret.arg2, op, nextarg])
    else
        ret = BinaryOpCall(ret, op, nextarg)
    end
    return ret
end

# Parse power (special case for preceding unary ops)
function parse_operator_power(ps::ParseState, @nospecialize(ret), op)
    nextarg = @precedence ps PowerOp - LtoR(PowerOp) @closer ps inwhere parse_expression(ps)

    if ret isa UnaryOpCall
        nextarg = BinaryOpCall(ret.arg, op, nextarg)
        ret = UnaryOpCall(ret.op, nextarg)
    else
        ret = BinaryOpCall(ret, op, nextarg)
    end
    return ret
end

# parse where
function parse_operator_where(ps::ParseState, @nospecialize(ret), op)
    nextarg = @precedence ps LazyAndOp @closer ps inwhere parse_expression(ps)

    if nextarg isa EXPR{Braces}
        args = nextarg.args
    else
        args = Any[nextarg]
    end
    return WhereOpCall(ret, op, args)
end

function parse_operator_dot(ps::ParseState, @nospecialize(ret), op)
    if ps.nt.kind == Tokens.LPAREN
        @static if VERSION > v"1.1-"
            iserred = ps.ws.kind != Tokens.EMPTY_WS
            sig = @default ps parse_call(ps, ret)
            nextarg = EXPR{TupleH}(sig.args[2:end])
            if iserred
                nextarg = ErrorToken(nextarg)
            end
        else
            sig = @default ps parse_call(ps, ret)
            nextarg = EXPR{TupleH}(sig.args[2:end])
        end
    elseif iskw(ps.nt) || ps.nt.kind == Tokens.IN || ps.nt.kind == Tokens.ISA || ps.nt.kind == Tokens.WHERE
        nextarg = IDENTIFIER(next(ps))
    elseif ps.nt.kind == Tokens.COLON
        op2 = OPERATOR(next(ps))
        if ps.nt.kind == Tokens.LPAREN
            nextarg = @closeparen ps @precedence ps DotOp - LtoR(DotOp) parse_expression(ps)
            nextarg = EXPR{Quote}(Any[op2, nextarg])
        else
            nextarg = @precedence ps DotOp - LtoR(DotOp) parse_unary(ps, op2)
        end
    elseif ps.nt.kind == Tokens.EX_OR && ps.nnt.kind == Tokens.LPAREN
        op2 = OPERATOR(next(ps))
        nextarg = parse_call(ps, op2)
    else
        nextarg = @precedence ps DotOp - LtoR(DotOp) parse_expression(ps)
    end

    if nextarg isa IDENTIFIER || nextarg isa EXPR{Vect} || (nextarg isa UnarySyntaxOpCall && is_exor(nextarg.arg1))
        ret = BinarySyntaxOpCall(ret, op, Quotenode(nextarg))
    elseif nextarg isa EXPR{MacroCall}
        mname = BinarySyntaxOpCall(ret, op, Quotenode(nextarg.args[1]))
        ret = EXPR{MacroCall}(Any[mname])
        for i = 2:length(nextarg.args)
            push!(ret, nextarg.args[i])
        end
    else
        ret = BinarySyntaxOpCall(ret, op, nextarg)
    end
    return ret
end

function parse_operator_anon_func(ps::ParseState, @nospecialize(ret), op)
    arg = @closer ps comma @precedence ps 0 parse_expression(ps)

    if !(arg isa EXPR{Begin} || (arg isa EXPR{InvisBrackets} && arg.args[2] isa EXPR{Block}))
        arg = EXPR{Block}(Any[arg])
    end
    return BinarySyntaxOpCall(ret, op, arg)
end

function parse_operator(ps::ParseState, @nospecialize(ret), op)
    K,dot = op.kind, op.dot
    P = precedence(K)

    if ret isa EXPR{ChainOpCall} && (is_star(op) || is_plus(op)) && op.kind == ret.args[2].kind
        nextarg = @precedence ps P - LtoR(P) parse_expression(ps)
        push!(ret, op)
        push!(ret, nextarg)
        ret = ret
    elseif ret isa BinaryOpCall && (is_star(op) || is_plus(op)) && op.kind == ret.op.kind && !ret.op.dot && ret.op.span > 0
        nextarg = @precedence ps P - LtoR(P) parse_expression(ps)
        ret = EXPR{ChainOpCall}(Any[ret.arg1, ret.op, ret.arg2, op, nextarg])
    elseif is_eq(op)
        ret = parse_operator_eq(ps, ret, op)
    elseif is_cond(op)
        ret = parse_operator_cond(ps, ret, op)
    elseif is_colon(op)
        ret = parse_operator_colon(ps, ret, op)
    elseif is_where(op)
        ret = parse_operator_where(ps, ret, op)
    elseif is_anon_func(op)
        ret = parse_operator_anon_func(ps, ret, op)
    elseif is_dot(op)
        ret = parse_operator_dot(ps, ret, op)
    elseif is_dddot(op) || is_prime(op)
        ret = UnarySyntaxOpCall(ret, op)
    elseif P == ComparisonOp
        ret = parse_comp_operator(ps, ret, op)
    elseif P == PowerOp
        ret = parse_operator_power(ps, ret, op)
    else
        ltor = K == Tokens.LPIPE ? true : LtoR(P)
        nextarg = @precedence ps P - ltor parse_expression(ps)

        if issyntaxcall(op)
            ret = BinarySyntaxOpCall(ret, op, nextarg)
        else
            ret = BinaryOpCall(ret, op, nextarg)
        end
    end
    return ret
end


text:
module Data

import SHA

using Random: shuffle!
using Base: @propagate_inbounds

export CMUDict, cmudict

deps(path...) = joinpath(@__DIR__, "..", "..", "deps", path...)

function download_and_verify(url, path, hash)
    tmppath = tempname()
    download(url, tmppath)
    hash_download = open(tmppath) do f
        bytes2hex(SHA.sha256(f))
    end
    if hash_download !== hash
        msg  = "Hash Mismatch!\n"
        msg *= "  Expected sha256:   $hash\n"
        msg *= "  Calculated sha256: $hash_download"
        error(msg)
    end
    mv(tmppath, path; force=true)
end

function __init__()
  mkpath(deps())
end

include("dataloader.jl")
export DataLoader

include("mnist.jl")
export MNIST

include("fashion-mnist.jl")
export FashionMNIST

include("cmudict.jl")
using .CMUDict

include("tree.jl")
include("sentiment.jl")
using .Sentiment

include("iris.jl")
export Iris

include("housing.jl")
export Housing

@deprecate DataLoader(x...; kws...) DataLoader(x; kws...)

end


text:
mutable struct CoupledArray{T,T2} <: AbstractArray{Float64,1}
  u::T
  u_control::T2
  order::Bool
end

Base.length(A::CoupledArray) = length(A.u) + length(A.u_control)
Base.size(A::CoupledArray) = (length(A),)
@inline function Base.getindex(A::CoupledArray,i::Int)
  if A.order == true
    i <= length(A.u) ? A.u[i] : A.u_control[i-length(A.u)]
  else
    i <= length(A.u) ? A.u_control[i] : A.u[i-length(A.u)]
  end
end

@inline function Base.getindex(A::CoupledArray,I...)
  A[CartesianIndices(A.u,I...)]
end

@inline function Base.getindex(A::CoupledArray,I::CartesianIndex{1})
  A[I[1]]
end

@inline Base.setindex!(A::CoupledArray,v,I...) = (A[CartesianIndices(A.u,I...)] = v)
@inline Base.setindex!(A::CoupledArray,v,I::CartesianIndex{1}) = (A[I[1]] = v)
@inline function Base.setindex!(A::CoupledArray,v,i::Int)
  if A.order == true
    i <= length(A.u) ? (A.u[i] = v) : (A.u_control[i-length(A.u)] = v)
  else
    i <= length(A.u) ? (A.u_control[i] = v) : (A.u[i-length(A.u)] = v)
  end
end

Base.IndexStyle(::Type{<:CoupledArray}) = IndexLinear()
Base.similar(A::CoupledArray) = CoupledArray(similar(A.u),similar(A.u_control),A.order)
Base.similar(A::CoupledArray,::Type{S}) where {S} = CoupledArray(similar(A.u,S),similar(A.u_control,S),A.order)

function recursivecopy!(dest::T, src::T) where T<:CoupledArray
  recursivecopy!(dest.u,src.u)
  recursivecopy!(dest.u_control,src.u_control)
  dest.order = src.order
end

add_idxs1(::Type{T},expr) where {T<:CoupledArray} = :($(expr).u)
add_idxs2(::Type{T},expr) where {T<:CoupledArray} = :($(expr).u_control)
@generated function Base.broadcast!(f,A::CoupledArray,B::Union{Number,CoupledArray}...)
  exs1 = ((add_idxs1(B[i],:(B[$i])) for i in eachindex(B))...,)
  exs2 = ((add_idxs2(B[i],:(B[$i])) for i in eachindex(B))...,)
  res = quote
      broadcast!(f,A.u,$(exs1...));broadcast!(f,A.u_control,$(exs2...))
    end
  res
end

Base.show(io::IO,A::CoupledArray) = show(io,A.u)
TreeViews.hastreeview(x::CoupledArray) = true
plot_indices(A::CoupledArray) = eachindex(A)
flip_u!(A::CoupledArray) = (A.order = !A.order)


text:
# ---
# title: 1030. Matrix Cells in Distance Order
# id: problem1030
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Sort
# link: <https://leetcode.com/problems/matrix-cells-in-distance-order/description/>
# hidden: true
# ---
#
# We are given a matrix with `R` rows and `C` columns has cells with integer
# coordinates `(r, c)`, where `0 <= r < R` and `0 <= c < C`.
#
# Additionally, we are given a cell in that matrix with coordinates `(r0, c0)`.
#
# Return the coordinates of all cells in the matrix, sorted by their distance
# from `(r0, c0)` from smallest distance to largest distance.  Here, the
# distance between two cells `(r1, c1)` and `(r2, c2)` is the Manhattan
# distance, `|r1 - r2| + |c1 - c2|`.  (You may return the answer in any order
# that satisfies this condition.)
#
#
#
# **Example 1:**
#
#
#
#     Input: R = 1, C = 2, r0 = 0, c0 = 0
#     Output: [[0,0],[0,1]]
#     Explanation: The distances from (r0, c0) to other cells are: [0,1]
#
#
# **Example 2:**
#
#
#
#     Input: R = 2, C = 2, r0 = 0, c0 = 1
#     Output: [[0,1],[0,0],[1,1],[1,0]]
#     Explanation: The distances from (r0, c0) to other cells are: [0,1,1,2]
#     The answer [[0,1],[1,1],[0,0],[1,0]] would also be accepted as correct.
#
#
# **Example 3:**
#
#
#
#     Input: R = 2, C = 3, r0 = 1, c0 = 2
#     Output: [[1,2],[0,2],[1,1],[0,1],[1,0],[0,0]]
#     Explanation: The distances from (r0, c0) to other cells are: [0,1,1,2,2,3]
#     There are other answers that would also be accepted as correct, such as [[1,2],[1,1],[0,2],[1,0],[0,1],[0,0]].
#
#
#
#
# **Note:**
#
#   1. `1 <= R <= 100`
#   2. `1 <= C <= 100`
#   3. `0 <= r0 < R`
#   4. `0 <= c0 < C`
#
#
## @lc code=start
using LeetCode

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


text:
module SpringCollab2020RemoveLightSources

using ..Ahorn, Maple

@mapdef Trigger "SpringCollab2020/RemoveLightSourcesTrigger" RemoveLightSourcesTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight,
    enableLightSources::Bool=false, fade::Bool=false)

const placements = Ahorn.PlacementDict(
    "Remove Light Sources (Spring Collab 2020)" => Ahorn.EntityPlacement(
        RemoveLightSourcesTrigger,
        "rectangle",
    ),
)

end


text:
cd(@__DIR__); allNGWPs = false; include("setups/toronto.jl")
gr(dpi = 200)

## pedestrian signal 16 most important LP-NGWP vectors
LP_NGWP = lp_ngwp(𝚽, Gstar_Sig.W, GP_dual; ϵ = 0.3)
fp = load("../datasets/new_toronto.jld", "fp")
G_Sig.f = reshape(fp, (N, 1))
dmatrix_LP = ngwp_analysis(G_Sig, LP_NGWP)
dvec_lp_ngwp, BS_lp_ngwp = ngwp_bestbasis(dmatrix_LP, GP_dual)
important_idx = sortperm(dvec_lp_ngwp[:].^2; rev = true)
println("================fp-LP-NGWP-top-basis-vectors=================")
for i in 1:16
    dr, dc = BS_lp_ngwp.levlist[important_idx[i]]
    w = LP_NGWP[dr, dc, :]
    j, k, l = NGWP_jkl(GP_dual, dr, dc)
    println("(j, k, l) = ($(j), $(k), $(l))")
    sgn = (maximum(w) > -minimum(w)) * 2 - 1
    gplot(A, X; width=1)
    scatter_gplot!(X; marker = sgn .* w, plotOrder = :s2l, ms = 3)
    plt = plot!(cbar = false, clims = (-0.075,0.075))
    savefig(plt,
        "../figs/Toronto_fp_DAG_LP_NGWP_ibv$(lpad(i,2,"0"))_j$(j)_k$(k)_l$(l).png")
end


text:
# Attempt to compute a reasonable default mu: at the starting
# position, the gradient of the input function should dominate the
# gradient of the barrier.
function initial_mu(gfunc::AbstractArray{T}, gbarrier::AbstractArray{T}, mu0factor::T = 0.001, mu0::T = convert(T, NaN)) where T
    if isnan(mu0)
        gbarriernorm = sum(abs, gbarrier)
        if gbarriernorm > 0
            mu = mu0factor*sum(abs, gfunc)/gbarriernorm
        else
            # Presumably, there is no barrier function
            mu = zero(T)
        end
    else
        mu = mu0
    end
    return mu
end

function barrier_box(g, x::AbstractArray{T}, l::AbstractArray{T}, u::AbstractArray{T}) where T
    calc_g = g !== nothing

    v = zero(T)
    @inbounds for i in eachindex(x)
        thisl = l[i]
        if isfinite(thisl)
            dx = x[i] - thisl
            if dx <= zero(T)
                return convert(T, Inf)
            end
            v -= log(dx)
            if calc_g
                g[i] = -one(T)/dx
            end
        else
            if calc_g
                g[i] = zero(T)
            end
        end
        thisu = u[i]
        if isfinite(thisu)
            dx = thisu - x[i]
            if dx <= zero(T)
                return convert(T, Inf)
            end
            v -= log(dx)
            if calc_g
                g[i] += one(T)/dx
            end
        end
    end
    return v
end

function function_barrier(gfunc, gbarrier, x::AbstractArray, f, fbarrier)
    vbarrier = fbarrier(gbarrier, x)
    return (isfinite(vbarrier) ? f(gfunc, x) : vbarrier), vbarrier
end

function barrier_combined(gfunc, gbarrier, g, x::AbstractArray,
                          fb, mu::Ref{<:Real})
    valfunc, valbarrier = fb(gbarrier, x, gfunc)
    if g !== nothing
        g .= gfunc .+ mu[].*gbarrier
    end
    return convert(eltype(x), valfunc + mu[]*valbarrier) # FIXME make this unnecessary
end

function limits_box(x::AbstractArray{T}, d::AbstractArray{T},
                    l::AbstractArray{T}, u::AbstractArray{T}) where T
    alphamax = convert(T, Inf)
    @inbounds for i in eachindex(x)
        if d[i] < 0
            alphamax = min(alphamax, ((l[i]-x[i])+eps(l[i]))/d[i])
        elseif d[i] > 0
            alphamax = min(alphamax, ((u[i]-x[i])-eps(u[i]))/d[i])
        end
    end
    epsilon = eps(max(alphamax, one(T)))
    if !isinf(alphamax) && alphamax > epsilon
        alphamax -= epsilon
    end
    return alphamax
end

# Default preconditioner for box-constrained optimization
# This creates the inverse Hessian of the barrier penalty
function precondprepbox!(P, x, l, u, mu)
    @. P.diag = 1/(mu[]*(1/(x-l)^2 + 1/(u-x)^2) + 1)
end

struct Fminbox{O<:AbstractOptimizer, T, P} <: AbstractConstrainedOptimizer
    method::O
    mu0::T
    mufactor::T
    precondprep::P
end

"""
# Fminbox
## Constructor
```julia
Fminbox(method;
        mu0=NaN,
        mufactor=0.0001,
        precondprep(P, x, l, u, mu) -> precondprepbox!(P, x, l, u, mu))
```
"""
function Fminbox(method::AbstractOptimizer = LBFGS();
                 mu0::Real = NaN, mufactor::Real = 0.001,
                 precondprep = (P, x, l, u, mu) -> precondprepbox!(P, x, l, u, mu))
    if method isa Newton || method isa NewtonTrustRegion
        throw(ArgumentError("Newton is not supported as the Fminbox optimizer."))
    end
    Fminbox(method, promote(mu0, mufactor)..., precondprep) # default optimizer
end

Base.summary(F::Fminbox) = "Fminbox with $(summary(F.method))"

# barrier_method() constructs an optimizer to solve the barrier problem using m = Fminbox.method as the reference.
# Essentially it only updates the P and precondprep fields of `m`.

# fallback
barrier_method(m::AbstractOptimizer, P, precondprep) =
    error("You need to specify a valid inner optimizer for Fminbox, $m is not supported. Please consult the documentation.")

barrier_method(m::ConjugateGradient, P, precondprep) =
    ConjugateGradient(eta = m.eta, alphaguess = m.alphaguess!,
                      linesearch = m.linesearch!, P = P,
                      precondprep = precondprep)

barrier_method(m::LBFGS, P, precondprep) =
    LBFGS(alphaguess = m.alphaguess!, linesearch = m.linesearch!, P = P,
          precondprep = precondprep)

barrier_method(m::GradientDescent, P, precondprep) =
    GradientDescent(alphaguess = m.alphaguess!, linesearch = m.linesearch!, P = P,
                    precondprep = precondprep)

barrier_method(m::Union{NelderMead, SimulatedAnnealing, ParticleSwarm, BFGS, AbstractNGMRES},
               P, precondprep) = m # use `m` as is

function optimize(f,
                  g,
                  l::AbstractArray{T},
                  u::AbstractArray{T},
                  initial_x::AbstractArray{T},
                  F::Fminbox = Fminbox(),
                  options = Options(); inplace = true, autodiff = :finite) where T<:AbstractFloat

    g! = inplace ? g : (G, x) -> copyto!(G, g(x))
    od = OnceDifferentiable(f, g!, initial_x, zero(T))

    optimize(od, l, u, initial_x, F, options)
end

function optimize(f,
                  l::AbstractArray{T},
                  u::AbstractArray{T},
                  initial_x::AbstractArray{T},
                  F::Fminbox = Fminbox(),
                  options = Options(); inplace = true, autodiff = :finite) where T<:AbstractFloat

    od = OnceDifferentiable(f, initial_x, zero(T); autodiff = autodiff)
    optimize(od, l, u, initial_x, F, options)
end

function optimize(
        df::OnceDifferentiable,
        l::AbstractArray{T},
        u::AbstractArray{T},
        initial_x::AbstractArray{T},
        F::Fminbox = Fminbox(),
        options = Options()) where T<:AbstractFloat

    outer_iterations = options.outer_iterations
    allow_outer_f_increases = options.allow_outer_f_increases
    show_trace, store_trace, extended_trace = options.show_trace, options.store_trace, options.extended_trace

    x = copy(initial_x)
    fbarrier = (gbarrier, x) -> barrier_box(gbarrier, x, l, u)
    fb = (gbarrier, x, gfunc) -> function_barrier(gfunc, gbarrier, x, df.fdf, fbarrier)
    gfunc = similar(x)
    gbarrier = similar(x)
    P = InverseDiagonal(similar(initial_x))
    # to be careful about one special case that might occur commonly
    # in practice: the initial guess x is exactly in the center of the
    # box. In that case, gbarrier is zero. But since the
    # initialization only makes use of the magnitude, we can fix this
    # by using the sum of the absolute values of the contributions
    # from each edge.
    boundaryidx = Vector{Int}()
    for i in eachindex(gbarrier)
        thisx = x[i]
        thisl = l[i]
        thisu = u[i]

        if thisx == thisl
            thisx = 0.99*thisl+0.01*thisu
            x[i] = thisx
            push!(boundaryidx,i)
        elseif thisx == thisu
            thisx = 0.01*thisl+0.99*thisu
            x[i] = thisx
            push!(boundaryidx,i)
        elseif thisx < thisl || thisx > thisu
            throw(ArgumentError("Initial x[$(Tuple(CartesianIndices(x)[i]))]=$thisx is outside of [$thisl, $thisu]"))
        end

        gbarrier[i] = (isfinite(thisl) ? one(T)/(thisx-thisl) : zero(T)) + (isfinite(thisu) ? one(T)/(thisu-thisx) : zero(T))
    end
    if length(boundaryidx) > 0
        @warn("Initial position cannot be on the boundary of the box. Moving elements to the interior.\nElement indices affected: $boundaryidx")
    end

    gradient!(df, x)
    gfunc .= gradient(df)

    mu = Ref(initial_mu(gfunc, gbarrier, T(F.mufactor), T(F.mu0)))

    # Use the barrier-aware preconditioner to define
    # barrier-aware optimization method instance (precondition relevance)
    _optimizer = barrier_method(F.method, P, (P, x) -> F.precondprep(P, x, l, u, mu))

    if show_trace > 0
        println("Fminbox")
        println("-------")
        print("Initial mu = ")
        show(IOContext(stdout, :compact=>true), "text/plain", mu[])
        println("\n")
    end

    g = similar(x)
    fval_all = Vector{Vector{T}}()

    # Count the total number of outer iterations
    iteration = 0

    # define the function (dfbox) to optimize by the inner optimizer
    funcc = (g, x) -> barrier_combined(gfunc, gbarrier, g, x, fb, mu)
    dfbox = OnceDifferentiable(x -> funcc(nothing, x),
                               (g, x) -> (funcc(g, x); g),
                               funcc, initial_x, zero(T))

    xold = similar(x)
    converged = false
    local results
    first = true
    fval0 = zero(T)

    while !converged && iteration < outer_iterations
        # Increment the number of steps we've had to perform
        iteration += 1

        copyto!(xold, x)
        # Optimize with current setting of mu
        fval0 = funcc(nothing, x)
        if show_trace > 0
            header_string = "Fminbox iteration $iteration"
            println(header_string)
            println("-"^length(header_string))
            print("Calling inner optimizer with mu = ")
            show(IOContext(stdout, :compact=>true), "text/plain", mu[])
            println("\n")
            println("(numbers below include barrier contribution)")
        end
        resultsnew = optimize(dfbox, x, _optimizer, options)
        if first
            results = resultsnew
            first = false
        else
            append!(results, resultsnew)
        end
        copyto!(x, minimizer(results))
        if show_trace > 0
            println()
            println("Exiting inner optimizer with x = ", x)
            print("Current distance to box: ")
            show(IOContext(stdout, :compact=>true), "text/plain", min(minimum(x-l), minimum(u-x)))
            println()
            println("Decreasing barrier term μ.\n")
        end

        # Decrease mu
        mu[] *= F.mufactor

        # Test for convergence
        g .= gfunc .+ mu[].*gbarrier

        results.x_converged, results.f_converged,
        results.g_converged, converged, f_increased = assess_convergence(x, xold, minimum(results), fval0, g,
                                                                         options.outer_x_tol, options.outer_f_tol, options.outer_g_tol)
        if f_increased && !allow_outer_f_increases
            @warn("f(x) increased: stopping optimization")
            break
        end
    end

    return MultivariateOptimizationResults(F, initial_x, minimizer(results), df.f(minimizer(results)),
            iteration, results.iteration_converged,
            results.x_converged, results.x_tol, norm(x - xold),
            results.f_converged, results.f_tol, f_abschange(minimum(results), fval0),
            results.g_converged, results.g_tol, norm(g, Inf),
            results.f_increased, results.trace, results.f_calls,
            results.g_calls, results.h_calls)
end


text:
# Cards
using PlayingCards
card = A♡
string(card)
suit(A♡)
rank(A♠)
high_value(A♢)
high_value(J♣)
low_value(A♡)

# Deck
using PlayingCards
deck = ordered_deck()
shuffle!(deck)
hand = pop!(deck, 2)

for card in deck
  @show card
end


text:
"""
    SphereBoxCollision

    collision between sphere and box

    origin_sphere:    position of sphere contact relative to body center of mass
    origin_box_a:     position of box corner contact a relative to body center of mass
    origin_box_b:     position of box corner contact b relative to body center of mass
    radius_sphere:    radius of sphere contact
"""
mutable struct SphereBoxCollision{T,O,I,OI} <: Collision{T,O,I,OI}
    origin_sphere::SVector{I,T}
    width::T
    depth::T
    height::T
    radius_sphere::T
end

# distance
function distance(collision::SphereBoxCollision, xp, qp, xc, qc)
    # contact origin points
    cop = contact_point_origin(xp, qp, collision.origin_sphere)
    coc = contact_point_box(cop, xc, qc, collision.width, collision.depth, collision.height)

    # distance between contact origins
    d = norm(cop - coc, 2)

    # minimum distance between spheres
    return d - collision.radius_sphere
end

function ∂distance∂x(gradient::Symbol, collision::SphereBoxCollision, xp, qp, xc, qc)
    # contact origin points
    cop = contact_point_origin(xp, qp, collision.origin_sphere)
    coc = contact_point_box(cop, xc, qc, collision.width, collision.depth, collision.height)

    # distance between contact origins
    d = norm(cop - coc, 2)
    ∂norm∂d = ∂norm∂x(cop - coc)

    if gradient == :parent
        D = ∂norm∂d * (1.0 * ∂contact_point_origin∂x(xp, qp, collision.origin_sphere) - ∂contact_point_box∂p(cop, xc, qc, collision.width, collision.depth, collision.height) * ∂contact_point_origin∂x(xp, qp, collision.origin_sphere))
    elseif gradient == :child
        D = ∂norm∂d * -1.0 * ∂contact_point_box∂x(cop, xc, qc, collision.width, collision.depth, collision.height)
    end

    # if gradient == :parent
    #     FD = FiniteDiff.finite_difference_jacobian(x -> distance(collision, x, qp, xc, qc), xp)
    # elseif gradient == :child
    #     FD = FiniteDiff.finite_difference_jacobian(x -> distance(collision, xp, qp, x, qc), xc)
    # end

    # return FD

    # @assert norm(D - FD, Inf) < 1.0e-5

    return D
end

function ∂distance∂q(gradient::Symbol, collision::SphereBoxCollision, xp, qp, xc, qc)
    # contact origin points
    cop = contact_point_origin(xp, qp, collision.origin_sphere)
    coc = contact_point_box(cop, xc, qc, collision.width, collision.depth, collision.height)

    # distance between contact origins
    d = norm(cop - coc, 2)
    ∂norm∂d = ∂norm∂x(cop - coc)

    if gradient == :parent
        D = ∂norm∂d * (1.0 * ∂contact_point_origin∂q(xp, qp, collision.origin_sphere) - ∂contact_point_box∂p(cop, xc, qc, collision.width, collision.depth, collision.height) * ∂contact_point_origin∂q(xp, qp, collision.origin_sphere))
    elseif gradient == :child
        D = ∂norm∂d * -1.0 * ∂contact_point_box∂q(cop, xc, qc, collision.width, collision.depth, collision.height)
    end

    # if gradient == :parent
    #     FD = FiniteDiff.finite_difference_jacobian(q -> distance(collision, xp, Quaternion(q..., false), xc, qc), vector(qp))
    # elseif gradient == :child
    #     FD = FiniteDiff.finite_difference_jacobian(q -> distance(collision, xp, qp, xc, Quaternion(q..., false)), vector(qc))
    # end

    # return FD

    # @assert norm(D - FD, Inf) < 1.0e-5

    return D
end

# contact point in world frame
function contact_point(relative::Symbol, collision::SphereBoxCollision, xp, qp, xc, qc)
    # contact origin points
    cop = contact_point_origin(xp, qp, collision.origin_sphere)
    coc = contact_point_box(cop, xc, qc, collision.width, collision.depth, collision.height)

    # direction of minimum distance (child to parent)
    d = cop - coc
    dir = normalize(d)

    # contact point
    if relative == :parent
        return cop - collision.radius_sphere * dir
    elseif relative == :child
        return coc
    end
end

function ∂contact_point∂x(relative::Symbol, jacobian::Symbol, collision::SphereBoxCollision, xp, qp, xc, qc)
    # contact origin points
    cop = contact_point_origin(xp, qp, collision.origin_sphere)
    coc = contact_point_box(cop, xc, qc, collision.width, collision.depth, collision.height)

    # direction of minimum distance (child to parent)
    d = cop - coc
    dir = normalize(d)

    if relative == :parent
        # cop - collision.radius_sphere * dir
        if jacobian == :parent
            ∂c∂x = ∂contact_point_origin∂x(xp, qp, collision.origin_sphere)
            X = ∂c∂x
            X -= collision.radius_sphere * ∂normalize∂x(d) * (∂c∂x - ∂contact_point_box∂p(cop, xc, qc, collision.width, collision.depth, collision.height) * ∂c∂x)
        elseif jacobian == :child
            X = -1.0 * collision.radius_sphere * ∂normalize∂x(d) * -1.0 * ∂contact_point_box∂x(cop, xc, qc, collision.width, collision.depth, collision.height)
        end
    elseif relative == :child
        # coc
        if jacobian == :parent
            X = ∂contact_point_box∂p(cop, xc, qc, collision.width, collision.depth, collision.height) * ∂contact_point_origin∂x(xp, qp, collision.origin_sphere)
        elseif jacobian == :child
            X = ∂contact_point_box∂x(cop, xc, qc, collision.width, collision.depth, collision.height)
        end
    end

    # if jacobian == :parent
    #     FD =  FiniteDiff.finite_difference_jacobian(x -> contact_point(relative, collision, x, qp, xc, qc), xp)
    # elseif jacobian == :child
    #     FD = FiniteDiff.finite_difference_jacobian(x -> contact_point(relative, collision, xp, qp, x, qc), xc)
    # end

    # return FD

    # @assert norm(X - FD, Inf) < 1.0e-5

    return X
end

function ∂contact_point∂q(relative::Symbol, jacobian::Symbol, collision::SphereBoxCollision, xp, qp, xc, qc)
    # contact origin points
    cop = contact_point_origin(xp, qp, collision.origin_sphere)
    coc = contact_point_box(cop, xc, qc, collision.width, collision.depth, collision.height)

    # direction of minimum distance (child to parent)
    d = cop - coc
    dir = normalize(d)

    if relative == :parent
        # cop - collision.radius_sphere * dir
        if jacobian == :parent
            ∂c∂q = ∂contact_point_origin∂q(xp, qp, collision.origin_sphere)
            Q = ∂c∂q
            Q -= collision.radius_sphere * ∂normalize∂x(d) * (∂c∂q - ∂contact_point_box∂p(cop, xc, qc, collision.width, collision.depth, collision.height) * ∂c∂q)
        elseif jacobian == :child
            Q = -1.0 * collision.radius_sphere * ∂normalize∂x(d) * -1.0 * ∂contact_point_box∂q(cop, xc, qc, collision.width, collision.depth, collision.height)
        end
    elseif relative == :child
        # coc
        if jacobian == :parent
            Q = ∂contact_point_box∂p(cop, xc, qc, collision.width, collision.depth, collision.height) * ∂contact_point_origin∂q(xp, qp, collision.origin_sphere)
        elseif jacobian == :child
            Q = ∂contact_point_box∂q(cop, xc, qc, collision.width, collision.depth, collision.height)
        end
    end

    # if jacobian == :parent
    #     FD = FiniteDiff.finite_difference_jacobian(q -> contact_point(relative, collision, xp, Quaternion(q..., false), xc, qc), vector(qp))
    # elseif jacobian == :child
    #     FD = FiniteDiff.finite_difference_jacobian(q -> contact_point(relative, collision, xp, qp, xc, Quaternion(q..., false)), vector(qc))
    # end

    # return FD

    # @assert norm(Q - FD, Inf) < 1.0e-5

    return Q
end

function contact_normal(collision::SphereBoxCollision, xp, qp, xc, qc)
    # contact points
    cop = contact_point(:parent, collision, xp, qp, xc, qc)
    coc = contact_point(:child,  collision, xp, qp, xc, qc)

    # unnormalized direction
    dir = cop - coc

    # distance
    dis = distance(collision, xp, qp, xc, qc)

    # normalized direction
    if dis >= 0.0
        return normalize(dir)'
    else
        return -1.0 * normalize(dir)'
    end
end

function ∂contact_normal_transpose∂x(jacobian::Symbol, collision::SphereBoxCollision, xp, qp, xc, qc)
    # contact origin points
    cop = contact_point(:parent, collision, xp, qp, xc, qc)
    coc = contact_point(:child,  collision, xp, qp, xc, qc)

    # unnormalized direction
    dir = cop - coc

    # Jacobians
    X = ∂normalize∂x(dir) * (∂contact_point∂x(:parent, jacobian, collision, xp, qp, xc, qc) - ∂contact_point∂x(:child, jacobian, collision, xp, qp, xc, qc))

    # distance
    dis = distance(collision, xp, qp, xc, qc)

    # normalized direction
    if dis >= 0.0
        return X
    else
        return -1.0 * X
    end
end

function ∂contact_normal_transpose∂q(jacobian::Symbol, collision::SphereBoxCollision, xp, qp, xc, qc)
    # contact origin points
    cop = contact_point(:parent, collision, xp, qp, xc, qc)
    coc = contact_point(:child,  collision, xp, qp, xc, qc)

    # unnormalized direction
    dir = cop - coc

    Q = ∂normalize∂x(dir) * (∂contact_point∂q(:parent, jacobian, collision, xp, qp, xc, qc) - ∂contact_point∂q(:child, jacobian, collision, xp, qp, xc, qc))

    # Jacobians
    # if jacobian == :parent
    #     FD = FiniteDiff.finite_difference_jacobian(q -> contact_normal(collision, xp, Quaternion(q..., false), xc, qc)', vector(qp))
    # elseif jacobian == :child
    #     FD = FiniteDiff.finite_difference_jacobian(q -> contact_normal(collision, xp, qp, xc, Quaternion(q..., false))', vector(qc))
    # end

    # distance
    dis = distance(collision, xp, qp, xc, qc)

    # @assert norm((dis >= 0.0 ? 1.0 : -1.0) * Q - FD, Inf) < 1.0e-2

    # normalized direction
    if dis >= 0.0
        return Q
    else
        return -1.0 * Q
    end
end


text:
# copied from https://github.com/fhs/CRC32.jl/blob/master/src/CRC32.jl to reduce
# dependencies

function maketable(poly::UInt32)
    tab = zeros(UInt32, 256)
    for i in 0:255
        crc = UInt32(i)
        for _ in 1:8
            if (crc&1) == 1
                crc = (crc >> 1) ⊻ poly
            else
                crc >>= 1
            end
        end
        tab[i+1] = crc
    end
    tab
end

const table = maketable(0xedb88320)

function crc32(data::Vector{UInt8}, crc::Integer=0)
    crc = ~UInt32(crc)
    for b in data
        crc = table[(UInt8(crc&0xff) ⊻ b) + 1] ⊻ (crc >> 8)
    end
    ~crc
end

crc32(data::AbstractString, crc::Integer=0) =
    crc32(convert(Vector{UInt8}, data), crc)


text:
# Autogenerated wrapper script for BVLS_jll for aarch64-linux-musl-libgfortran4
export libbvls

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("BVLS")
JLLWrappers.@declare_library_product(libbvls, "libbvls.so")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libbvls,
        "lib/libbvls.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Base.Test

function f22938(a, b, x...)
    d = 1
    a = d
    for i in 1:b
        d += i
    end
    return i * a
end

msig = Tuple{typeof(f22938),Int,Int,Int,Int}
world = typemax(UInt)
_, msp, m = Base._methods_by_ftype(msig, -1, world)[]
mi = Core.Inference.code_for_method(m, msig, msp, world, false)
c0 = Core.Inference.retrieve_code_info(mi)

@test isempty(Core.Inference.validate_code(mi))
@test isempty(Core.Inference.validate_code(c0))

# INVALID_EXPR_HEAD
c = Core.Inference.copy_code_info(c0)
insert!(c.code, 4, Expr(:(=), SlotNumber(2), Expr(:invalid, 1)))
errors = Core.Inference.validate_code(c)
@test length(errors) == 1
@test errors[1].kind === Core.Inference.INVALID_EXPR_HEAD

# INVALID_LVALUE
c = Core.Inference.copy_code_info(c0)
insert!(c.code, 4, Expr(:(=), LabelNode(1), 1))
insert!(c.code, 2, Expr(:(=), :x, 1))
insert!(c.code, 10, Expr(:(=), 3, 1))
errors = Core.Inference.validate_code(c)
@test length(errors) == 3
@test all(e.kind === Core.Inference.INVALID_LVALUE for e in errors)

# INVALID_RVALUE
c = Core.Inference.copy_code_info(c0)
insert!(c.code, 2, Expr(:(=), SlotNumber(2), GotoNode(1)))
insert!(c.code, 4, Expr(:(=), SlotNumber(2), LabelNode(2)))
insert!(c.code, 10, Expr(:(=), SlotNumber(2), LineNumberNode(2)))
for h in (:gotoifnot, :line, :const, :meta)
    push!(c.code, Expr(:(=), SlotNumber(2), Expr(h)))
end
errors = Core.Inference.validate_code(c)
@test length(errors) == 10
@test count(e.kind === Core.Inference.INVALID_RVALUE for e in errors) == 7
@test count(e.kind === Core.Inference.INVALID_EXPR_NARGS for e in errors) == 3

# INVALID_CALL_ARG/INVALID_EXPR_NARGS
c = Core.Inference.copy_code_info(c0)
insert!(c.code, 2, Expr(:(=), SlotNumber(2), Expr(:call, :+, SlotNumber(2), GotoNode(1))))
insert!(c.code, 4, Expr(:call, :-, Expr(:call, :sin, LabelNode(2)), 3))
insert!(c.code, 10, Expr(:call, LineNumberNode(2)))
for h in (:gotoifnot, :line, :const, :meta)
    push!(c.code, Expr(:call, :f, Expr(h)))
end
errors = Core.Inference.validate_code(c)
@test length(errors) == 10
@test count(e.kind === Core.Inference.INVALID_CALL_ARG for e in errors) == 7
@test count(e.kind === Core.Inference.INVALID_EXPR_NARGS for e in errors) == 3

# EMPTY_SLOTNAMES
c = Core.Inference.copy_code_info(c0)
empty!(c.slotnames)
errors = Core.Inference.validate_code(c)
@test length(errors) == 2
@test any(e.kind === Core.Inference.EMPTY_SLOTNAMES for e in errors)
@test any(e.kind === Core.Inference.SLOTFLAGS_MISMATCH for e in errors)

# SLOTFLAGS_MISMATCH
c = Core.Inference.copy_code_info(c0)
push!(c.slotnames, :dummy)
errors = Core.Inference.validate_code(c)
@test length(errors) == 1
@test errors[1].kind === Core.Inference.SLOTFLAGS_MISMATCH

# SLOTTYPES_MISMATCH
c = @code_typed(f22938(1,2,3,4))[1]
pop!(c.slottypes)
errors = Core.Inference.validate_code(c)
@test length(errors) == 1
@test errors[1].kind === Core.Inference.SLOTTYPES_MISMATCH

# SLOTTYPES_MISMATCH_UNINFERRED
c = Core.Inference.copy_code_info(c0)
c.slottypes = 1
errors = Core.Inference.validate_code(c)
@test length(errors) == 1
@test errors[1].kind === Core.Inference.SLOTTYPES_MISMATCH_UNINFERRED

# SSAVALUETYPES_MISMATCH
c = @code_typed(f22938(1,2,3,4))[1]
empty!(c.ssavaluetypes)
errors = Core.Inference.validate_code(c)
@test length(errors) == 1
@test errors[1].kind === Core.Inference.SSAVALUETYPES_MISMATCH

# SSAVALUETYPES_MISMATCH_UNINFERRED
c = Core.Inference.copy_code_info(c0)
c.ssavaluetypes -= 1
errors = Core.Inference.validate_code(c)
@test length(errors) == 1
@test errors[1].kind === Core.Inference.SSAVALUETYPES_MISMATCH_UNINFERRED

# INVALID_ASSIGNMENT_SLOTFLAG
c = Core.Inference.copy_code_info(c0)
c.slotflags[8] = 0x00
errors = Core.Inference.validate_code(c)
@test length(errors) == 1
@test errors[1].kind === Core.Inference.INVALID_ASSIGNMENT_SLOTFLAG

# SIGNATURE_NARGS_MISMATCH
old_sig = mi.def.sig
mi.def.sig = Tuple{1,2}
errors = Core.Inference.validate_code(mi)
mi.def.sig = old_sig
@test length(errors) == 1
@test errors[1].kind === Core.Inference.SIGNATURE_NARGS_MISMATCH

# NON_TOP_LEVEL_METHOD
c = Core.Inference.copy_code_info(c0)
push!(c.code, Expr(:method, :dummy))
errors = Core.Inference.validate_code(c)
@test length(errors) == 1
@test errors[1].kind === Core.Inference.NON_TOP_LEVEL_METHOD

# SLOTNAMES_NARGS_MISMATCH
mi.def.nargs += 20
errors = Core.Inference.validate_code(mi)
mi.def.nargs -= 20
@test length(errors) == 2
@test count(e.kind === Core.Inference.SLOTNAMES_NARGS_MISMATCH for e in errors) == 1
@test count(e.kind === Core.Inference.SIGNATURE_NARGS_MISMATCH for e in errors) == 1

max_stars_repo_name: Mikewl/julia

max_issues_repo_name: Mikewl/julia

max_forks_repo_name: Mikewl/julia


text:
using LazySets

function fiveDimSys_specification()
    # initial set
    X0 = BallInf([1.0, 0.0, 0.0, 0.0, 0.0], 0.1)

    # time horizon: 5 time units
    time_horizon = 5.0

    # specification
    O = Dict{Symbol, Any}(:T => time_horizon)

    return X0, O
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
function always_true(_...)
    return true
end

function always_false(_...)
    return false
end

function always_nothing(_...)
    return nothing
end

function passthrough(args...)
    return (args, )
end


text:
using CreditMetrics
using Test

using StatsFuns

@testset "ProbabilityVector" begin
    # Test ProbabilityVector basics
    @test probabilities(ProbabilityVector([1.0])) == [1.0]
    @test probabilities(ProbabilityVector([0.25, 0.0, 0.5, 0.25])) == [0.25, 0.0, 0.5, 0.25]
    @test probabilities(ProbabilityVector([1])) == [1]
    @test probabilities(ProbabilityVector([Float32(1.0)])) == [Float32(1.0)]
    @test probabilities(ProbabilityVector([Float32(0.25), 0.75])) == [Float32(0.25), 0.75]

    # Test ProbabilityVector from AbstractProbabilityVector using ProbabilityVector
    @test probabilities(ProbabilityVector(ProbabilityVector([1.0]))) == [1.0]
    @test probabilities(ProbabilityVector(ProbabilityVector([0.25, 0.0, 0.5, 0.25]))) == [0.25, 0.0, 0.5, 0.25]
    @test probabilities(ProbabilityVector(ProbabilityVector([1]))) == [1]
    @test probabilities(ProbabilityVector(ProbabilityVector([Float32(1.0)]))) == [Float32(1.0)]
    @test probabilities(ProbabilityVector(ProbabilityVector([Float32(0.25), 0.75]))) == [Float32(0.25), 0.75]

    # Test ProbabilityVector errors
    @test_throws DomainError ProbabilityVector([])
    @test_throws DomainError ProbabilityVector([0.2])
    @test_throws DomainError ProbabilityVector([1.2])
    @test_throws DomainError ProbabilityVector([0.3, 0.6])
    @test_throws DomainError ProbabilityVector([0.7, 0.6])
    @test_throws DomainError ProbabilityVector([1.75, -0.75])

    # Test conversion from AbstractUnitIntervalThresholds to ProbabilityVector
    @test probabilities(ProbabilityVector(UnitIntervalThresholds(Float64[]))) == [1.0]
    @test probabilities(ProbabilityVector(UnitIntervalThresholds([0.0]))) == [0.0, 1.0]
    @test probabilities(ProbabilityVector(UnitIntervalThresholds([1.0]))) == [1.0, 0.0]
    @test probabilities(ProbabilityVector(UnitIntervalThresholds([0.25]))) == [0.25, 0.75]
    @test probabilities(ProbabilityVector(UnitIntervalThresholds([0.25, 0.5]))) == [0.25, 0.25, 0.5]
end

@testset "UnitIntervalThresholds" begin
    # Test UnitIntervalThresholds basics
    @test thresholds(UnitIntervalThresholds([])) == []
    @test thresholds(UnitIntervalThresholds([0.25])) == [0.25]
    @test thresholds(UnitIntervalThresholds([0.25, 0.75, 0.875])) == [0.25, 0.75, 0.875]
    @test thresholds(UnitIntervalThresholds([0.25, 0.25, 0.875])) == [0.25, 0.25, 0.875]
    @test thresholds(UnitIntervalThresholds([0, 0, 1, 1])) == [0, 0, 1, 1]
    @test thresholds(UnitIntervalThresholds([0.2, 1])) == [0.2, 1]

    # Test UnitIntervalThresholds from AbstractUnitIntervalThresholds via UnitIntervalThresholds
    @test thresholds(UnitIntervalThresholds(UnitIntervalThresholds([]))) == []
    @test thresholds(UnitIntervalThresholds(UnitIntervalThresholds([0.25]))) == [0.25]
    @test thresholds(UnitIntervalThresholds(UnitIntervalThresholds([0.25, 0.75, 0.875]))) == [0.25, 0.75, 0.875]
    @test thresholds(UnitIntervalThresholds(UnitIntervalThresholds([0.25, 0.25, 0.875]))) == [0.25, 0.25, 0.875]
    @test thresholds(UnitIntervalThresholds(UnitIntervalThresholds([0, 0, 1, 1]))) == [0, 0, 1, 1]
    @test thresholds(UnitIntervalThresholds(UnitIntervalThresholds([0.2, 1]))) == [0.2, 1]

    # Test UnitIntervalThresholds errors
    @test_throws DomainError UnitIntervalThresholds([-0.5])
    @test_throws DomainError UnitIntervalThresholds([-0.5, 0.5])
    @test_throws DomainError UnitIntervalThresholds([1.5])
    @test_throws DomainError UnitIntervalThresholds([0.5, 1.5])
    @test_throws DomainError UnitIntervalThresholds([0.75, 0.25])

    # Test conversions from AbstractProbabilityVector to UnitIntervalThresholds
    @test thresholds(UnitIntervalThresholds(ProbabilityVector([1.0]))) == Float64[]
    @test thresholds(UnitIntervalThresholds(ProbabilityVector([0.25, 0.75]))) == [0.25]
    @test thresholds(UnitIntervalThresholds(ProbabilityVector([0.25, 0.0, 0.75]))) == [0.25, 0.25]
end

@testset "RealLineThresholds" begin
    # Test RealLineTrhesholds basics
    @test thresholds(RealLineThresholds([])) == []
    @test thresholds(RealLineThresholds([-1.3])) == [-1.3]
    @test thresholds(RealLineThresholds([-1.3, 1.5])) == [-1.3, 1.5]
end

@testset "DeformedRealLineThresholds" begin
    # Test DeformedRealLineThresholds
    @test thresholds(DeformedRealLineThresholds(RealLineThresholds([-0.5, 0.75]), 1.75, 1.0)) == [1.25, 2.5]
end

@testset "NormalThresholdLink" begin
    # NormalThresholdLink basic
    @test cdf(NormalThresholdLink{Float64}())(1.5) == normcdf(1.5)
    @test cdf(NormalThresholdLink{Float32}())(Float32(1.5)) == normcdf(Float32(1.5))
    @test invcdf(NormalThresholdLink{Float64}())(0.3) == norminvcdf(0.3)
    @test invcdf(NormalThresholdLink{Float32}())(Float32(0.3)) == norminvcdf(Float32(0.3))

    # NormalThresholdLink errors
    @test_throws TypeError NormalThresholdLink{String}()
    @test_throws TypeError NormalThresholdLink{Int}()
end