text:
module logging
using Dates: Dates
using Logging
using InteractiveUtils
using ProgressLogging: asprogress

import Term: Theme,
            theme,
            textlen,
            escape_brackets,
            unescape_brackets,
            reshape_text,
            has_markup,
            int

import ..box: ROUNDED
import ..style: apply_style
import ..renderables: AbstractRenderable
import ..Tprint: tprintln
import ..progress: ProgressBar,
                ProgressJob,
                DescriptionColumn,
                ProgressColumn,
                render,
                start!,
                stop!,
                addjob!,
                removejob!,
                update!,
                getjob,
                SeparatorColumn,
                PercentageColumn

import ..console: console_width, console_height, change_scroll_region, move_to_line

export TermLogger, install_term_logger

DEFAULT_LOGGER = global_logger()

# ---------------------------- create term logger ---------------------------- #
"""
    TermLogger

Custom logger type.
"""
struct TermLogger <: Logging.AbstractLogger
    io::IO
    theme::Theme
    pbar::ProgressBar

end
TermLogger(theme::Theme) = TermLogger(stderr, theme)

TermLogger(io::IO, theme::Theme) = TermLogger(
    io,
    theme,
    ProgressBar(;
        transient=true,
        columns=[DescriptionColumn,
        SeparatorColumn,
        ProgressColumn,
        SeparatorColumn,
        PercentageColumn
        ])
)

# set logger beavior
Logging.min_enabled_level(logger::TermLogger) = Logging.Debug
function Logging.shouldlog(logger::TermLogger, level, _module, group, id)
    return true
end
Logging.catch_exceptions(logger::TermLogger) = true

# --------------------------- handle logger message -------------------------- #

"""
    print_closing_line(color::String, width::Int)

Print the final line of a log message with style and date info
"""
function print_closing_line(color::String, width::Int = 48)
    tprintln(
        "  [$color bold dim]$(ROUNDED.bottom.left)" *
        "$(ROUNDED.row.mid)"^(width) *
        "[/$color bold dim]",
    )
    _date = Dates.format(Dates.now(), "e, dd u yyyy")
    _time = Dates.format(Dates.now(), "HH:MM:SS")
    pad = width - textlen(_date * _time) - 2
    return tprintln(
        " "^pad * "[dim]$_date[/dim] [bold dim underline]$_time[/bold dim underline]"
    )
end

"""
    handle_progress(logger::TermLogger, prog)

Handle progress information passed by `ProgressLogging`

It creates/adds/removes `ProgressJob`s to the logger's
`ProgressBar` to create progress visualizations.
"""
function handle_progress(logger::TermLogger, prog)
    pbar = logger.pbar
    pbar.running || start!(pbar)

    # if progress not started yet, ignore
    if isnothing(prog.fraction) || prog.fraction == 0
        return
    end

    # check if a job exists for this progress, or add one
    job = getjob(pbar, prog.id)
    job = isnothing(job) ? addjob!(
                pbar;
                description=length(prog.name) > 0 ? prog.name : "running...",
                N = 100,
                transient=false,
                id = prog.id
            ) : job

    # if done, remove job.
    if prog.done || prog.fraction == 1.0
        job.i = 100
        stop!(job)
    else
        update!(job; i = (Int ∘ floor)(prog.fraction * 100))
    end

    # render
    if all(map(j -> j.finished, pbar.jobs))
        map(j -> removejob!(pbar, j), pbar.jobs)
        stop!(pbar)
    else
        render(pbar)
    end
end

"""
    Logging.handle_message(logger::TermLogger,

Handle printing of log messages, with style!.

In addition to the log message and info such as file/line and time of log,
it prints kwargs styled by their type.
"""
function Logging.handle_message(
    logger::TermLogger, lvl, msg, _mod, group, id, file, line; kwargs...
)
    _progress = asprogress(lvl, msg, _mod, group, id, file, line; kwargs...)
    isnothing(_progress) || return handle_progress(logger, _progress)

    # get name of function where logging message was called from
    fname = ""
    for frame in stacktrace()
        if "$(frame.file)" == file && frame.line == line
            fname = ".[$(logger.theme.func) underline]$(frame.func)[/$(logger.theme.func) underline]"
        end
    end

    # prepare styles
    if lvl == Logging.Info
        color = logger.theme.info
    elseif lvl == Logging.Debug
        color = logger.theme.debug
    elseif lvl == Logging.Warn
        color = logger.theme.warn
    elseif lvl == Logging.Error
        color = logger.theme.error
    else
        color = "#90CAF9"
    end

    outline_markup = "$color dim"
    hor = "[$outline_markup]▶[/$outline_markup]"
    vert = "[$outline_markup]" * ROUNDED.mid.left * "[/$outline_markup]"

    # style message
    if msg isa AbstractString
        msg = has_markup(msg) ? msg : "[#8abeff]$msg[/#8abeff]"
        msg = reshape_text(msg, 64)
    else
        msg = string(msg)
    end

    # get the first line of information
    content = "[$color underline bold]@$(string(lvl))[/$color underline bold] [#edb46f ]($(_mod)$fname):[/#edb46f ]"

    # for multi-lines message, print each line separately.
    msg_lines::Vector{AbstractString} = split(msg, "\n")
    for n in 2:length(msg_lines)
        msg_lines[n] = "  $vert   " * " "^textlen(content) * "[#8abeff]" * msg_lines[n]
    end
    content *= "  " * msg_lines[1]
    tprintln(content)
    tprintln.(msg_lines[2:end])

    # if no kwargs we can just quit
    if length(kwargs) == 0
        print_closing_line(color)
        return nothing
    end

    # get padding width
    _types = string.(typeof.([v for v in values(kwargs)]))
    for (i, _type) in enumerate(values(kwargs))
        if typeof(_type) <: Function
            _types[i] = string(Function)
        end
    end

    wpad = max(textlen.(_types)...) + 2
    namepad = max(textlen.(string.([v for v in keys(kwargs)]))...)

    # print all kwargs
    tprintln("  $vert")
    for ((k, v), _type) in zip(kwargs, _types)
        # get line stub
        pad = wpad - textlen(_type)
        line =
            "  $vert [$(theme.type) dim]($(_type))[/$(theme.type) dim]" *
            " "^pad *
            hor *
            "  [bold #ff9959]$k[/bold #ff9959]"

        epad = namepad - textlen(string(k))
        line *= " "^epad * " [bold red]=[/bold red] "
        lpad = textlen(line) - 4

        # get value style
        if v isa Number
            _style = logger.theme.number
        elseif v isa Symbol
            _style = logger.theme.symbol
        elseif v isa AbstractString
            _style = logger.theme.string
        elseif v isa AbstractVector
            _style = logger.theme.number
            _size = size(v)
            v = escape_brackets(string(v))
            v = textlen(v) > 33 ? v[1:30] * "..." : v
            v *= "\n [dim]size: $(_size)[/dim]"
        elseif v isa AbstractArray || v isa AbstractMatrix
            _style = logger.theme.number
            v =
                "$(typeof(v)) [dim]<: $(supertypes(typeof(v))[end-1])[/dim]" *
                "\n [dim]size: $(size(v))[/dim]"
        elseif v isa Function
            _style = logger.theme.func
        elseif v isa AbstractRenderable
            _style = "default"
            v = "$(typeof(v)) \e[2m(size: $(v.measure))\e[0m"
        else
            _style = nothing
        end

        # print value lines
        vlines = split(string(v), "\n")

        if !isnothing(_style)
            vlines = map(
                ln -> "["*_style*"]"*ln*"[/"*_style*"]", vlines
            )
        end

        if length(vlines) == 1
            tprintln(line * vlines[1])
        else
            tprintln(line * vlines[1])
            for ln in vlines[2:end]
                tprintln("  $vert " * " "^lpad * ln)
            end
        end
    end
    return print_closing_line(color)
end

# ---------------------------- install term logger --------------------------- #

"""
    install_term_logger(theme::Theme=theme)

Install `TermLogger` as the global logging system.

`theme::Theme` can be passed to specify the theme to use for styling objects.
"""
function install_term_logger(theme::Theme = theme)
    _logger = TermLogger(theme)
    return global_logger(_logger)
end

function uninstall_term_logger()
    _lg = global_logger(DEFAULT_LOGGER)
    logger = global_logger(DEFAULT_LOGGER)
    return logger
end
end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations3D(1.4)

initial_condition = initial_condition_convergence_test
source_terms = source_terms_convergence_test

volume_flux  = flux_ranocha
solver = DGMulti(element_type = Hex(),
                 approximation_type = periodic_derivative_operator(
                   derivative_order=1, accuracy_order=4, xmin=0.0, xmax=1.0, N=20),
                 surface_flux = flux_lax_friedrichs,
                 volume_integral = VolumeIntegralFluxDifferencing(volume_flux))

mesh = DGMultiMesh(solver, coordinates_min=(-1.0, -1.0, -1.0),
                           coordinates_max=( 1.0,  1.0,  1.0))

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
                                    source_terms=source_terms)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 0.4)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(solver))

alive_callback = AliveCallback(analysis_interval=analysis_interval)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback)

###############################################################################
# run the simulation

sol = solve(ode, RDPK3SpFSAL49(), abstol = 1.0e-7, reltol = 1.0e-7,
            save_everystep = false, callback = callbacks)
summary_callback() # print the timer summary

ext: jl


text:
some_elm_v1(x) = zero(x)
some_elm_v1(::Type{Missing}) = missing
some_elm_v1(::Type{String}) = ""

compress_then_write_v1(b, io) = compress_then_write(eltype(b), b, io)

compress_then_write_v1(::Type{Missing}, b, io) = (len = 0, type = Missing)

compress_then_write_v1(T, b, io) = begin
    bbc = Blosc.compress(b)
    res = length(bbc)
    write(io, bbc)
    return (len = res, type = T)
end

"""
    save a DataFrames to the outdir
"""
savejdf_v1(outdir, df) =begin
    if VERSION < v"1.3.0-rc1"
        return ssavejdf(outdir, df)
    end

    pmetadatas = Any[missing for i = 1:length(names(df))]

    if !isdir(outdir)
        mkpath(outdir)
    end

    # use a bounded channel to limit
    c1 = Channel{Bool}(Threads.nthreads())
    atexit(()->close(c1))

    for (i, n) in enumerate(names(df))
        put!(c1, true)
        pmetadatas[i] = @spawn begin
            io = BufferedOutputStream(open(joinpath(outdir, string(n)) ,"w"))
            res = compress_then_write(Array(df[!, i]), io)
            close(io)
            res
        end
        take!(c1)
    end

    metadatas = fetch.(pmetadatas)

    fnl_metadata = (
        names = names(df),
        rows = size(df, 1),
        metadatas = metadatas
    )

    serialize(joinpath(outdir, "metadata.jls"), fnl_metadata)
    fnl_metadata
end

"""
    serially save a DataFrames to the outdir
"""
ssavejdf_v1(outdir, df::DataFrame) = begin
    pmetadatas = Any[missing for i = 1:length(names(df))]

    if !isdir(outdir)
        mkpath(outdir)
    end

    for i = 1:length(names(df))
        io = BufferedOutputStream(open(joinpath(outdir, string(names(df)[i])), "w"))
        pmetadatas[i] = compress_then_write(Array(df[!, i]), io)
        close(io)
    end

    fnl_metadata = (
        names = names(df),
        rows = size(df, 1),
        metadatas = pmetadatas
    )

    serialize(joinpath(outdir, "metadata.jls"), fnl_metadata)
    fnl_metadata
end

# figure out from metadata how much space is allocated
get_bytes_v1(metadata) = begin
    if metadata.type == String
        return max(
            metadata.string_compressed_bytes,
            metadata.string_len_bytes,
        )
    elseif metadata.type == Missing
        return 0
    elseif metadata.type >: Missing
        return max(
            get_bytes(metadata.Tmeta),
            get_bytes(metadata.missingmeta),
        )
    else
        return metadata.len
    end
end

hasfieldnames(::Type{T}) where {T} = fieldnames(T) >= 1

column_loader_v1(t::Type{T}, io, metadata) where {T} = begin
    buffer = Vector{UInt8}(undef, metadata.len)
    column_loader!(buffer, t, io, metadata)
end

# load bytes bytes from io decompress into type
column_loader!_v1(buffer, ::Type{T}, io, metadata) where {T} = begin
    readbytes!(io, buffer, metadata.len)
    return Blosc.decompress(T, buffer)
end

compress_then_write_v1(::Type{Bool}, b, io) = begin
    b8 = UInt8.(b)
    bbc = Blosc.compress(b8)
    write(io, bbc)
    return (len = length(bbc), type = Bool)
end

column_loader_v1(T::Type{Bool}, io, metadata) = begin
# Bool are saved as UInt8
    buffer = Vector{UInt8}(undef, metadata.len)
    readbytes!(io, buffer, metadata.len)
    Bool.(Blosc.decompress(UInt8, buffer))
end

column_loader!_v1(buffer, T::Type{Bool}, io, metadata) = begin
    # Bool are saved as UInt8
    read!(io, buffer)
    res = Blosc.decompress(UInt8, buffer)
    Bool.(res)
end

compress_then_write_v1(::Type{T}, b, io) where {T>:Missing} = begin
    S = nonmissingtype(eltype(b))
    b_S = coalesce.(b, some_elm(S))

    metadata = compress_then_write(S, b_S, io)

    b_m = ismissing.(b)
    metadata2 = compress_then_write(eltype(b_m), b_m, io)

    (
     Tmeta = metadata,
     missingmeta = metadata2,
     type = eltype(b),
     len = max(metadata.len, metadata2.len),
    )
end

column_loader!_v1(buffer, ::Type{Union{Missing,T}}, io, metadata) where {T} = begin
    # read the content
    Tmeta = metadata.Tmeta

    t_pre = column_loader!(buffer, Tmeta.type, io, Tmeta) |> allowmissing
    #t = t_pre
    # read the missings as bool
    m = column_loader(Bool, io, metadata.missingmeta)
    #return t_pre
    t_pre[m] .= missing
    t_pre
end

compress_then_write_v1(_, b::StringVector{T}, io) where {T} = compress_then_write(b, io)

compress_then_write_v1(b::StringVector{T}, io) where {T} = begin
    fields = (getfield(b, n) for n in fieldnames(typeof(b)))
    (
     metadata = [(eltype(f), write(io, Blosc.compress(f))) for f in fields],
     type = typeof(b),
    )
end

column_loader_v1(::Type{StringVector{T}}, io, metadata) where {T} = begin
    # uncompress
    args = Vector[]

    # assign the buffer once
    buffer = Vector{UInt8}(undef, maximum(x -> x[2], metadata.metadata))

    for (elm_type, compressed_bytes) in metadata.metadata
        readbytes!(io, buffer, compressed_bytes)
        push!(args, Blosc.decompress(elm_type, buffer))
    end
    metadata.type(args...)
end

# perform a RLE
compress_then_write_v1(::Type{String}, b::Array{String}, io) = begin
    # write the string one by one
    # do a Run-length encoding (RLE)
    previous_b = b[1]
    cnt = 1
    lens = Int[]
    str_lens = Int[]
    for i = 2:length(b)
        if b[i] != previous_b
            push!(str_lens, write(io, previous_b))
            push!(lens, cnt)
            #push!(str_lens, sizeof(previous_b))
            cnt = 0
            previous_b = b[i]
        end
        cnt += 1
    end

    # reach the end: two situation
    # 1) it's a new element, so write it
    # 2) it's an existing element. Also write it
    push!(str_lens, write(io, previous_b))
    push!(lens, cnt)
    #push!(str_lens, sizeof(previous_b))

    @assert sum(lens) == length(b)

    str_lens_compressed = Blosc.compress(Vector{UInt32}(str_lens))
    str_lens_bytes = write(io, str_lens_compressed)

    lens_compressed = Blosc.compress(Vector{UInt64}(lens))
    rle_bytes = write(io, lens_compressed)

    # return metadata
    return (
        string_compressed_bytes = sum(str_lens),
        string_len_bytes = str_lens_bytes,
        rle_bytes = rle_bytes,
        rle_len = length(str_lens),
        type = String,
        len = max(sum(str_lens), str_lens_bytes, rle_bytes),
    )
end

# load a string column
"""
    metadata should consists of length, compressed byte size of string-lengths,
    string content lengths
"""
column_loader_v1!(_, ::Type{String}, io, metadata) = begin
    buffer = Vector{UInt8}(undef, metadata.string_compressed_bytes)
    readbytes!(io, buffer, metadata.string_compressed_bytes)
    #return String(buffer)

    # read the string-lengths
    buffer2 = Vector{UInt8}(undef, metadata.string_len_bytes)
    readbytes!(io, buffer2, metadata.string_len_bytes)

    buffer3 = Vector{UInt8}(undef, metadata.rle_bytes)
    readbytes!(io, buffer3, metadata.rle_bytes)

    counts = Blosc.decompress(UInt64, buffer3)

    str_lens = Blosc.decompress(UInt32, buffer2)

    #return (String(buffer), str_lens, counts)

    lengths = inverse_rle(str_lens, counts)
    offsets = inverse_rle(vcat(0, cumsum(str_lens[1:end-1])), counts)

    #res = StringArray{String, 1}(buffer, vcat(1, cumsum(Blosc.decompress(UInt64, buffer3))[1:end-1]) .-1,  )
    res = StringArray{String,1}(buffer, offsets, lengths)
end

column_loader_v1!(buffer, ::Type{Missing}, io, metadata) = nothing


text:
using GLMakie

width = 3
height = 3

dir = [225 180 135 270 0 90 315 360 45]

Wy = -sind.(mod.(-dir .+ 90, 360))
Wx = -cosd.(mod.(-dir .+ 90, 360))

data1 = Array{Float32}(undef, round(Int, width), round(Int, height))

let k = 1
    for j::Int in 1:height
        for i::Int in 1:width
            data1[i, j] = k
            k = k + 1
        end
    end
end
fig, ax, hm = heatmap(data1)
Colorbar(fig[:, end+1], hm)

xs = 1:width
ys = 1:height
arrows!(fig[:, 1], xs, ys, Wx, Wy, arrowsize = 30, lengthscale = 0.2, linewidth = 3)
display(fig)


text:
include("07_two_site_operator.jl")

println("
#######################################################
# Tutorial 8                                          #
#                                                     #
# 2-site custom operator/gate definition              #
#######################################################
")

import ITensors: op

# Controlled-Ry rotation gate
function op(::OpName"CRy", ::SiteType"S=1/2"; θ)
  return [
    1 0 0           0
    0 1 0           0
    0 0 cos(θ / 2) -sin(θ / 2)
    0 0 sin(θ / 2) cos(θ / 2)
  ]
end

# CRy(θ=π/2)
θ = π / 2
CRyθ = op("CRy", i1, i2; θ)

# |10⟩ = |Z-Z+⟩
ZmZp = Zm1 * Zp2

# CRyθ|Z-Z+⟩ = |Z-X+⟩
@show apply(CRyθ, ZmZp) ≈ Zm1 * Xp2


text:
function readinput(filename)
    map(group -> filter(!isempty, split(replace(group, "\n" => " "), " ")), split(read(filename, String), "\n\n"))
end

passports = readinput("d04.in")

P1, P2 = 0, 0

function lenientvalid(passport)
    valid = true
    for field in ["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"]
        if !haskey(passport, field)
            valid = false
        end
    end
    return valid
end

function strictvalid(passport)
    valid = true
    if !(1920 <= parse(Int64, get(passport, "byr", "")) <= 2002)
        valid = false
    end
    if !(2010 <= parse(Int64, get(passport, "iyr", "")) <= 2020)
        valid = false
    end
    if !(2020 <= parse(Int64, get(passport, "eyr", "")) <= 2030)
        valid = false
    end
    hgt = get(passport, "hgt", "999cm")
    if endswith(hgt,"cm")
        if !(150 <= parse(Int64, hgt[1:end-2]) <= 193)
            valid = false
        end
    elseif endswith(hgt,"in")
        if !(59 <= parse(Int64, hgt[1:end-2]) <= 76)
            valid = false
        end
    else
        valid = false
    end
    if match(r"#[0-9a-f]{6}$", get(passport, "hcl", "NOT")) === nothing
        valid = false
    end
    if !(get(passport, "ecl", "NOT") in ["amb", "blu", "brn", "gry", "grn", "hzl", "oth"])
        valid = false
    end
    if match(r"^[0-9]{9}$", get(passport, "pid", "NOT")) === nothing
        valid = false
    end
    return valid
end

for passportlist in passports
    passport = Dict()
    for field in passportlist
        k, v = split(field, ':')
        passport[k] = v
    end
    if lenientvalid(passport)
        global P1 += 1
        if strictvalid(passport)
            global P2 += 1
        end
    end
end

println("P1: ", P1)
println("P2: ", P2)

max_stars_repo_name: eduellery/adventofcode

max_issues_repo_name: eduellery/adventofcode

max_forks_repo_name: eduellery/adventofcode


text:
struct NotAlmostZeroException <: Exception
    msg::AbstractString
    value::Union{AbstractArray, Number}
    rtol::Real
    atol::Real
end

function Base.showerror(io::IO, e::NotAlmostZeroException)
    n = @sprintf "%.4g" norm(e.value)
    print(e.msg, " :",
          " norm = ", n,
          " rtol = ", e.rtol,
          " atol = ", e.atol)
    if ! get(io, :compact, false)
        println(io)
        print(io, "value = ", e.value)
    end
end

function find_errors(sweep::ContinuationSweep,
                     prob::AbstractContinuationProblem,
                     opts::ContinuationOptions;
                     sweep_index = '?')
    if isempty(sweep.u)
        return []
    end
    cache = get_prob_cache(prob)
    H = _similar(sweep.u[1], length(sweep.u[1]) - 1)
    errors = []
    for (i, u) in enumerate(sweep.u)
        H = residual!(H, u, cache)
        if ! isalmostzero(H, opts.atol)
            push!(errors, NotAlmostZeroException(
                "H(sweeps[$sweep_index].u[$i]) not almost zero",
                H,
                opts.rtol,
                opts.atol,
            ))
        end
    end
    return errors
end

function find_errors(sol::ContinuationSolution, args...)
    errors = []
    for (k, sweep) in enumerate(sol.sweeps)
        append!(errors, find_errors(sweep, args...; sweep_index=k))
    end
    return errors
end

function find_errors(solver::ContinuationSolver)
    errors = find_errors(solver.sol, solver.prob, solver.opts)

    # Create a problem cache since solver.cache.prob_cache must not be
    # modified.  `find_errors` has to be a pure function!
    prob_cache = get_prob_cache(solver.prob)

    for (k, sweep) in enumerate(solver.sol.sweeps)
        if all(isindomain(u, prob_cache) for u in sweep.u)
            push!(errors, ErrorException(
                "$k-th sweep terminates before touching the boundary.",
            ))
        end
        # TODO: allow unbounded domain
    end
    return errors
end

print_errors(errors) = print_errors(stdout, errors)

function print_errors(io, errors)
    for e in errors
        showerror(io, e)
        println(io)
    end
end

# TODO: define `diagnosis` function for user-friendly API


text:
# ____________________________________________________________________
# powerlaw
#
mutable struct powerlaw <: AbstractComponent
    norm::Parameter
    x0::Parameter
    alpha::Parameter

    function powerlaw(x0::Number)
        out = new(Parameter(1),
                  Parameter(x0),
                  Parameter(-1))
        out.norm.low = 0
        out.x0.low = 0
        out.x0.fixed = true
        out.alpha.low = -5; out.alpha.high = 5
        return out
    end
end

function evaluate!(buffer, comp::powerlaw, x::Domain{1},
                   norm, x0, alpha)
    buffer .= norm .* (x ./ x0).^alpha
end


text:
using HOODESolver
using LinearAlgebra
using Test

function testcoeffloat()
    tab=[ (big"0.1", big"0.0001"), (big"0.00001", big"0.01")]
    @testset "test coefexp float vs bigfloat" begin
        for (bigeps, bigdt) in tab
            epsilon = convert(Float64,bigeps)
            dt = convert(Float64, bigeps)
            bigceab = CoefExpAB(10,bigeps,32,bigdt)
            ceab = CoefExpAB(10,epsilon,32,bigdt)
            @test isapprox(bigceab.tab_coef, ceab.tab_coef, atol=1e-15, rtol=1e-14)
    #       println("epsilon=$epsilon norm=$(norm(bigceab.tab_coef-ceab.tab_coef, Inf))")
        end
    end
 end
function testcomp()
    @testset "test solve float vs bigfloat" begin

        bigeps = big"0.0001"
        epsilon = convert(Float64,bigeps)
        bigu0=rand(BigFloat,4)
        u0 = convert.(Float64, bigu0)
        A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
        fct = (u,p,t) -> [0, u[4], -2u[1]*u[2], -u[2]-u[1]^2+u[2]^2] # Henon-Heiles
        bigprob = HOODEProblem(fct, bigu0, (big"0.0",big"3.0"), missing, A, bigeps)
        prob = HOODEProblem(fct, u0, (0.0, 3.0), missing, A, epsilon)
        bigsol = solve(bigprob)
        sol = solve(prob)

        for i=1:101
            @test isapprox(
    bigsol.interp.u_caret[i],
    sol.interp.u_caret[i],
    atol=1e-13,
    rtol=1e-12
)
        end
    end
end

testcoeffloat()
testcomp()


text:
# Local compute version

using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
# Pkg.precompile()

Base.cd(ENV["HOME"]*"/learning-odometry")

# load the necessary Python bindings
include(joinpath(@__DIR__,"PyTensorFlowUsage.jl"))

Base.cd(@__DIR__)

## Load all required packages
using Distributed
# addprocs(5) # make sure there are 4 processes waiting before loading packages

@everywhere begin
  using Pkg
  Pkg.activate(@__DIR__)
end

WP = WorkerPool(2:nprocs() |> collect )

# @show ARGS
include("parsecommands.jl")

using DelimitedFiles
using Dates, Statistics
using YAML
using CoordinateTransformations, Rotations, StaticArrays
using ImageCore
using Images, ImageDraw
# using ImageView
using AprilTags
using DataInterpolations

using RoME
using Caesar
0
@everywhere using Caesar
@everywhere using JLD2

@everywhere begin
# using Fontconfig
# using Compose
using Gadfly
using RoMEPlotting
using Cairo
# using FileIO
# using GeometryTypes # using MeshCat
end

# setup configuration
include(joinpath(@__DIR__,"configParameters.jl") )

# Figure export folder
global currdirtime = now()
if parsed_args["previous"] != ""
  # currdirtime = "2018-11-07T01:36:52.274"
  currdirtime = parsed_args["previous"]
end

resultsparentdir = joinpath(datadir, "results")
resultsdir = joinpath(resultsparentdir, "$(currdirtime)")

if parsed_args["previous"] == ""
  # When running fresh from new data
  include(joinpath(@__DIR__,"createResultsDirs.jl"))
end

# Utils required for this processing script
include(joinpath(@__DIR__,"racecarUtils.jl") )
include(joinpath(@__DIR__,"cameraUtils.jl") )
# include(joinpath(dirname(@__FILE__),"visualizationUtils.jl") )

@everywhere function plotRacecarInterm(fgl::G, resultsdirl, psyml::Symbol)::Nothing where G <: AbstractDFG
  @show ls(fgl)
  pl = drawPosesLandms(fgl, spscale=0.1, drawhist=false, meanmax=:mean) #,xmin=-3,xmax=6,ymin=-5,ymax=2);
  Gadfly.draw(PNG(joinpath(resultsdirl, "images", "$(psyml).png"),15cm, 10cm),pl)
  pl = drawPosesLandms(fgl, spscale=0.1, meanmax=:mean) # ,xmin=-3,xmax=3,ymin=-2,ymax=2);
  Gadfly.draw(PNG(joinpath(resultsdirl, "images", "hist_$(psyml).png"),15cm, 10cm),pl)
  # pl = plotPose2Vels(fgl, Symbol("$(psyml)"), coord=Coord.Cartesian(xmin=-1.0, xmax=1.0))
  # Gadfly.draw(PNG(joinpath(resultsdirl, "images", "vels_$(psyml).png"),15cm, 10cm),pl)
  # save combined image with tags
  nothing
end

global tag_bag = Dict()
if parsed_args["previous"] == ""
  tag_bag = detectTagsInImgs(datafolder, imgfolder, resultsdir, camidxs, iterposes=parsed_args["iterposes"])
  # save the tag bag file for future use
  @save resultsdir*"/tag_det_per_pose.jld2" tag_bag
else
  @load resultsdir*"/tag_det_per_pose.jld2" tag_bag
end

## Load joystick time series data

runnr = parse(Int, parsed_args["folder_name"][end])
joyTimeseries = readdlm(joinpath(datadir,parsed_args["folder_name"],"labRun$(runnr)_joy.csv"), ',')
joyTimeseries = joyTimeseries[2:end,[1,6,8]]
joyTimeseries[:,1] .= joyTimeseries[:,1]*1e-9 # .|> unix2datetime;

# load the detections file to get timestamps
detcData = readdlm(joinpath(datadir,parsed_args["folder_name"],"labRun$(runnr)Detections.csv"), ',')
detcData = detcData[2:end,:]
detcPoseTs = detcData[:,4]*1e-9 .+ detcData[:,3] .|> unix2datetime

## find timeseries segments that go with each pose
joyTsDict = Dict{Symbol, Array{Float64,2}}()
poseTimes = Dict{Symbol, DateTime}(:x0 => detcPoseTs[1])
mask = joyTimeseries[:,1] .< datetime2unix(detcPoseTs[1])
joyTsDict[:x0] = joyTimeseries[mask,:]

for i in 2:length(detcPoseTs)
  mask = datetime2unix(detcPoseTs[i-1]) .<= joyTimeseries[:,1] .< datetime2unix(detcPoseTs[i])
  joyTsDict[Symbol("x$(i-1)")] = joyTimeseries[mask,:]
  poseTimes[Symbol("x$(i-1)")] = detcPoseTs[i-1] # later helper
end

## interpolate up to PyQuest values...

intJoyDict = Dict{Symbol,Vector{Vector{Float64}}}()
# for sym, lclJD = :x1, joyTsDict[:x1]
for (sym, lclJD) in joyTsDict
  if 1 < size(lclJD,1)
    tsLcl = range(lclJD[1,1],lclJD[end,1],length=25)
    intrTrTemp = DataInterpolations.LinearInterpolation(lclJD[:,2],lclJD[:,1])
    intrStTemp = DataInterpolations.LinearInterpolation(lclJD[:,3],lclJD[:,1])
    newVec = Vector{Vector{Float64}}()
    for tsL in tsLcl
      newVal = zeros(4)
      newVal[1] = intrTrTemp(tsL)
      newVal[2] = intrStTemp(tsL)
      push!(newVec, newVal)
    end
    # currently have no velocity values
    intJoyDict[sym] = newVec
  else
    intJoyDict[sym] = [zeros(4) for i in 1:25]
  end
end

# add the NeuralPose2Pose2 factor in Main workspace
include( joinpath(dirname(pathof(Caesar)), "..", "examples", "learning", "hybrid", "NeuralPose2Pose2", "PyNeuralPose2Pose2.jl") )
# include(joinpath(@__DIR__, "NeuralPose2Pose2/PyNeuralPose2Pose2.jl"))

## run the solution

fg = main(WP, resultsdir, camidxs, tag_bag, jldfile=parsed_args["jldfile"], failsafe=parsed_args["failsafe"], show=parsed_args["show"], odopredfnc=PyTFOdoPredictorPoint2, joyvel=intJoyDict, poseTimes=poseTimes, multiproc=false  )

## development

#
# 0
#
# fg = initfg()
#
# addVariable!(fg, :x0, Pose2, timestamp=poseTimes[:x0])
# addVariable!(fg, :x1, Pose2, timestamp=poseTimes[:x1])
#
# pp2 = PriorPose2(MvNormal(zeros(3),LinearAlgebra.diagm([0.01;0.01;0.005].^2)))
# addFactor!(fg, [:x0], pp2)
#
# # the neural factor
# DXmvn = MvNormal(zeros(3),LinearAlgebra.diagm([0.01;0.01;0.005].^2))
# odopredfnc=PyTFOdoPredictorPoint2
# joyvel=intJoyDict
# pp = PyNeuralPose2Pose2(odopredfnc,joyvel[:x0],DXmvn,0.4)
# addFactor!(fg, [:x0;:x1], pp)

# fails
# key 1 not found
# julia101 -p 4 apriltag_and_zed_slam.jl --previous "2018-11-09T01:42:33.279" --jldfile "racecar_fg_x299.jld2" --folder_name "labrun7" --failsafe

#


text:
@testset "WorldClim BioClim" begin
    using RasterDataSources: rastername, zipurl, zipname, zippath

    zip_url = URI(scheme="https", host="biogeo.ucdavis.edu", path="/data/worldclim/v2.1/base/wc2.1_10m_bio.zip")
    @test zipurl(WorldClim{BioClim}, 2; res="10m") == zip_url
    @test zipname(WorldClim{BioClim}, 2; res="10m") == "wc2.1_10m_bio.zip"

    raster_file = joinpath(ENV["RASTERDATASOURCES_PATH"], "WorldClim", "BioClim", "wc2.1_10m_bio_2.tif")
    @test rasterpath(WorldClim{BioClim}, 2; res="10m") == raster_file
    @test getraster(WorldClim{BioClim}, 2; res="10m") == raster_file
    @test getraster(WorldClim{BioClim}, (2,), "10m") == (raster_file,)
    @test isfile(raster_file)
end


text:
using QuasiArrays, Base64, Test

@testset "QuasiSubArray" begin
    @testset "basics" begin
        A = QuasiArray(reshape(1:60, 5,4,3), (range(0;stop=1,length=5), Base.OneTo(4), [2,3,6]))
        sA = view(A, 0.25:0.25, 1:4, :)
        @test @inferred(size(sA)) == (1, 4, 3)
        @test parent(sA) == A
        @test parentindices(sA) == (0.25:0.25, 1:4, Inclusion([2,3,6]))
        @test size(sA) == (1, 4, 3)
        @test axes(sA) == (Base.OneTo(1), Base.OneTo(4), Inclusion([2,3,6]))
        @test sA[1, 2, [2,3,6]][:] == [7,27,47]
        sA = view(A, 0.0:0.25:1, 1:4, 6)
        sA[1:5,1:4] .= -2
        @test all(A[:,:,6] .== -2)
        fill!(sA, -3)
        @test all(A[:,:,6] .== -3)
        sA = view(A, 0.0:0.25:1, 3:3, 2:3)
        @test size(sA) == (5,1,2)
        @test axes(sA) === (Base.OneTo(5), Base.OneTo(1), Base.OneTo(2))
        @test convert(AbstractArray,sA) == convert(AbstractArray{Int},sA) ==
                Array(sA) == Array{Int}(sA) == A[0.0:0.25:1,3:3,2:3]

        @test Matrix(view(sA,1:1,1:1,1))[1] == sA[1,1,1] == 11
        # Test with mixed types
        @test sA[:, Int16(1), [big(1),2]] == [11:15 31:35]
        sA = view(A, 0.25:0.25, 1:4, [2 3; 2 6])
        @test ndims(sA) == 4
        @test axes(sA) === (Base.OneTo(1), Base.OneTo(4), Base.OneTo(2), Base.OneTo(2))
        @test_throws BoundsError view(A, 0:0.25:0.25, 3, [1 3; 4 2])
        sA = view(A, 0:0.25:0.25, 3, [2 3; 2 6])
        @test ndims(sA) == 3
        @test axes(sA) === (Base.OneTo(2), Base.OneTo(2), Base.OneTo(2))
    end

    @testset "FastSubQuasiArray" begin
        a = QuasiArray(randn(5), (1:0.5:3,))
        v = view(a, Inclusion(1.5:0.5:2))
        @test axes(v) == (Inclusion(1.5:0.5:2),)
        @test v[1.5] == a[1.5]
        @test_throws BoundsError v[1]
        v[1.5] = 3
        @test a[1.5] === 3.0

        a = QuasiArray(randn(2), ([[1,2],[3,4]],))
        v = view(a, Inclusion([[1,2],[3,4]]))
        @test v[[1,2]] == a[[1,2]]
        v[[1,2]] = 3
        @test a[[1,2]] == 3.0
    end

    # @testset "logical indexing #4763" begin
    #     A = view(QuasiVector(1:5,range(0,2;length=5)), 0.5:0.5:1.5)
    #     A[Inclusion(0.5:0.5:1.0)]
    #     @test A[A.<3] == view(A, A.<7) == [5, 6]
    #     @test Base.unsafe_getindex(A, A.<7) == [5, 6]
    #     B = reshape(1:16, 4, 4)
    #     sB = view(B, 2:3, 2:3)
    #     @test sB[sB.>8] == view(sB, sB.>8) == [10, 11]
    #     @test Base.unsafe_getindex(sB, sB.>8) == [10, 11]
    # end
    #
    # # Tests where dimensions are dropped
    # A = copy(reshape(1:120, 3, 5, 8))
    # sA = view(A, 2, :, 1:8)
    # @test parent(sA) == A
    # @test parentindices(sA) == (2, Base.Slice(1:5), 1:8)
    # @test Base.parentdims(sA) == [2:3;]
    # @test size(sA) == (5, 8)
    # @test axes(sA) === (Base.OneTo(5), Base.OneTo(8))
    # @test @inferred(strides(sA)) == (3,15)
    # @test sA[2, 1:8][:] == [5:15:120;]
    # @test sA[:,1] == [2:3:14;]
    # @test sA[2:5:end] == [5:15:110;]
    # sA[2:5:end] .= -1
    # @test all(sA[2:5:end] .== -1)
    # @test all(A[5:15:120] .== -1)
    # test_bounds(sA)
    # sA = view(A, 1:3, 1:5, 5)
    # @test Base.parentdims(sA) == [1:2;]
    # @test size(sA) == (3,5)
    # @test axes(sA) === (Base.OneTo(3),Base.OneTo(5))
    # @test @inferred(strides(sA)) == (1,3)
    # test_bounds(sA)
    # sA = view(A, 1:2:3, 3, 1:2:8)
    # @test Base.parentdims(sA) == [1,3]
    # @test size(sA) == (2,4)
    # @test axes(sA) === (Base.OneTo(2), Base.OneTo(4))
    # @test @inferred(strides(sA)) == (2,30)
    # @test sA[:] == A[sA.indices...][:]
    # test_bounds(sA)
    #
    # let a = [5:8;]
    #     @test parent(a) == a
    #     @test parentindices(a) == (1:4,)
    # end
    #
    # # issue #6218 - logical indexing
    # A = rand(2, 2, 3)
    # msk = fill(true, 2, 2)
    # msk[2,1] = false
    # sA = view(A, :, :, 1)
    # sA[msk] .= 1.0
    # @test sA[msk] == fill(1, count(msk))
    #
    # # bounds checking upon construction; see #4044, #10296
    # @test_throws BoundsError view(1:10, 8:11)
    # A = reshape(1:20, 5, 4)
    # sA = view(A, 1:2, 1:3)
    # @test_throws BoundsError view(sA, 1:3, 1:3)
    # @test_throws BoundsError view(sA, 1:2, 1:4)
    # view(sA, 1:2, 1:2)
    # @test_throws BoundsError view(A, 17:23)
    # view(A, 17:20)
    #
    # # Linear indexing by one multidimensional array:
    # A = reshape(1:120, 3, 5, 8)
    # sA = view(A, :, :, :)
    # @test sA[[72 17; 107 117]] == [72 17; 107 117]
    # @test sA[[99 38 119 14 76 81]] == [99 38 119 14 76 81]
    # @test sA[[fill(1, (2, 2, 2)); fill(2, (2, 2, 2))]] == [fill(1, (2, 2, 2)); fill(2, (2, 2, 2))]
    # sA = view(A, 1:2, 2:3, 3:4)
    # @test sA[(1:8)'] == [34 35 37 38 49 50 52 53]
    # @test sA[[1 2 4 4; 6 1 1 4]] == [34 35 38 38; 50 34 34 38]
    #
    # # issue #11871
    # let a = fill(1., (2,2)),
    #     b = view(a, 1:2, 1:2)
    #     b[2] = 2
    #     @test b[2] === 2.0
    # end
    #
    # # issue #15138
    # let a = [1,2,3],
    #     b = view(a, UInt(1):UInt(2))
    #     @test b == view(a, UInt(1):UInt(2)) == view(view(a, :), UInt(1):UInt(2)) == [1,2]
    # end
    #
    # let A = reshape(1:4, 2, 2),
    #     B = view(A, :, :)
    #     @test parent(B) === A
    #     @test parent(view(B, 0x1, :)) === parent(view(B, 0x1, :)) === A
    # end
    #
    # # issue #15168
    # let A = rand(10), sA = view(copy(A), :)
    #     @test sA[Int16(1)] === sA[Int32(1)] === sA[Int64(1)] === A[1]
    #     permute!(sA, Vector{Int16}(1:10))
    #     @test A == sA
    # end
    #
    # # the following segfaults with LLVM 3.8 on Windows, ref #15417
    # @test Array(view(view(reshape(1:13^3, 13, 13, 13), 3:7, 6:6, :), 1:2:5, :, 1:2:5)) ==
    #     cat([68,70,72],[406,408,410],[744,746,748]; dims=3)
    #
    # # tests @view (and replace_ref_end!)
    # X = reshape(1:24,2,3,4)
    # Y = 4:-1:1
    #
    # @test isa(@view(X[1:3]), SubArray)
    #
    # @test X[1:end] == @.(@view X[1:end]) # test compatibility of @. and @view
    # @test X[1:end-3] == @view X[1:end-3]
    # @test X[1:end,2,2] == @view X[1:end,2,2]
    # @test X[1,1:end-2,1] == @view X[1,1:end-2,1]
    # @test X[1,2,1:end-2] == @view X[1,2,1:end-2]
    # @test X[1,2,Y[2:end]] == @view X[1,2,Y[2:end]]
    # @test X[1:end,2,Y[2:end]] == @view X[1:end,2,Y[2:end]]
    #
    # u = (1,2:3)
    # @test X[u...,2:end] == @view X[u...,2:end]
    # @test X[(1,)...,(2,)...,2:end] == @view X[(1,)...,(2,)...,2:end]
    #
    # # test macro hygiene
    # let size=(x,y)-> error("should not happen"), Base=nothing
    #     @test X[1:end,2,2] == @view X[1:end,2,2]
    # end
    #
    # # test that side effects occur only once
    # let foo = [X]
    #     @test X[2:end-1] == @view (push!(foo,X)[1])[2:end-1]
    #     @test foo == [X, X]
    # end
    #
    # # test @views macro
    # @views let f!(x) = x[1:end-1] .+= x[2:end].^2
    #     x = [1,2,3,4]
    #     f!(x)
    #     @test x == [5,11,19,4]
    #     @test x[1:3] isa SubArray
    #     @test x[2] === 11
    #     @test Dict((1:3) => 4)[1:3] === 4
    #     x[1:2] .= 0
    #     @test x == [0,0,19,4]
    #     x[1:2] .= 5:6
    #     @test x == [5,6,19,4]
    #     f!(x[3:end])
    #     @test x == [5,6,35,4]
    #     x[Y[2:3]] .= 7:8
    #     @test x == [5,8,7,4]
    #     x[(3,)..., ()...] += 3
    #     @test x == [5,8,10,4]
    #     i = Int[]
    #     # test that lhs expressions in update operations are evaluated only once:
    #     x[push!(i,4)[1]] += 5
    #     @test x == [5,8,10,9] && i == [4]
    #     x[push!(i,3)[end]] += 2
    #     @test x == [5,8,12,9] && i == [4,3]
    #     @. x[3:end] = 0       # make sure @. works with end expressions in @views
    #     @test x == [5,8,0,0]
    # end
    # @views @test isa(X[1:3], SubArray)
    # @test X[1:end] == @views X[1:end]
    # @test X[1:end-3] == @views X[1:end-3]
    # @test X[1:end,2,2] == @views X[1:end,2,2]
    # @test X[1,2,1:end-2] == @views X[1,2,1:end-2]
    # @test X[1,2,Y[2:end]] == @views X[1,2,Y[2:end]]
    # @test X[1:end,2,Y[2:end]] == @views X[1:end,2,Y[2:end]]
    # @test X[u...,2:end] == @views X[u...,2:end]
    # @test X[(1,)...,(2,)...,2:end] == @views X[(1,)...,(2,)...,2:end]
    # # test macro hygiene
    # let size=(x,y)-> error("should not happen"), Base=nothing
    #     @test X[1:end,2,2] == @views X[1:end,2,2]
    # end
    #
    # # issue #18034
    # # ensure that it is possible to create an isbits, IndexLinear view of an immutable Array
    # let
    #     struct ImmutableTestArray{T, N} <: Base.DenseArray{T, N}
    #     end
    #     Base.size(::Union{ImmutableTestArray, Type{ImmutableTestArray}}) = (0, 0)
    #     Base.IndexStyle(::Union{ImmutableTestArray, Type{ImmutableTestArray}}) = Base.IndexLinear()
    #     a = ImmutableTestArray{Float64, 2}()
    #     @test Base.IndexStyle(view(a, :, :)) == Base.IndexLinear()
    #     @test isbits(view(a, :, :))
    # end
    #
    # # ref issue #17351
    # @test @inferred(reverse(view([1 2; 3 4], :, 1), dims=1)) == [3, 1]
    #
    # let
    #     s = view(reshape(1:6, 2, 3), 1:2, 1:2)
    #     @test @inferred(s[2,2,1]) === 4
    # end
    #
    # # issue #18581: slices with OneTo axes can be linear
    # let
    #     A18581 = rand(5, 5)
    #     B18581 = view(A18581, :, axes(A18581,2))
    #     @test IndexStyle(B18581) === IndexLinear()
    # end
    #
    # @test sizeof(view(zeros(UInt8, 10), 1:4)) == 4
    # @test sizeof(view(zeros(UInt8, 10), 1:3)) == 3
    # @test sizeof(view(zeros(Float64, 10, 10), 1:3, 2:6)) == 120
    #
    # # PR #25321
    # # checks that issue in type inference is resolved
    # A = rand(5,5,5,5)
    # V = view(A, 1:1 ,:, 1:3, :)
    # @test @inferred(strides(V)) == (1, 5, 25, 125)
    #
    # # Issue #26263 — ensure that unaliascopy properly trims the array
    # A = rand(5,5,5,5)
    # V = view(A, 2:5, :, 2:5, 1:2:5)
    # @test @inferred(Base.unaliascopy(V)) == V == A[2:5, :, 2:5, 1:2:5]
    # @test @inferred(sum(Base.unaliascopy(V))) ≈ sum(V) ≈ sum(A[2:5, :, 2:5, 1:2:5])
    #
    # # issue #27632
    # function _test_27632(A)
    #     for J in CartesianIndices(size(A)[2:end])
    #         A[1, J]
    #     end
    #     nothing
    # end
    # # check that this doesn't crash
    # _test_27632(view(ones(Int64, (1, 1, 1)), 1, 1, 1))
    #
    # # issue #29608 - views of single values can be considered contiguous
    # @test Base.iscontiguous(view(ones(1), 1))

    @testset "array subviews" begin
        A = QuasiArray(reshape(1:60, 5,4,3), (range(0;stop=1,length=5), Base.OneTo(4), [2,3,6]))
        @test A[0.0,:,3] isa Array
        @test view(A,0.0,:,3) isa SubArray
        @test A[0.0,:,3] == view(A,0.0,:,3)
    end

    @testset "sub-of-sub" begin
        A = QuasiArray(randn(3,3),(0:0.5:1,1:0.5:2))
        V = view(A,0:0.5:0.5,:)
        V2 = view(V,:,1:0.5:1.5)
        @test V2 isa SubArray
        @test parent(V2) isa QuasiArray
        @test V2 == A[0:0.5:0.5,1:0.5:1.5]

        @test QuasiArrays.sub_materialize(A) == A
    end

    @testset "Show" begin
        a = QuasiVector(randn(3), 0:0.5:1)
        v = view(a, :)
        if VERSION < v"1.7-"
            @test stringmime("text/plain", v) == "view(QuasiVector{Float64, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}, Inclusion(0.0:0.5:1.0)) with eltype Float64"
        else
            @test stringmime("text/plain", v) == "view(QuasiVector{Float64, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, $Int}}}, Inclusion(0.0:0.5:1.0)) with eltype Float64"
        end
    end

    @testset "==" begin
        A = QuasiArray(randn(3,3),(0:0.5:1,1:0.5:2))
        B = QuasiArray(randn(4,3),(0:0.5:1.5,1:0.5:2))
        @test A == view(A,:,:)
        @test view(A,:,:) == A
        @test view(A,:,:) == view(A,:,:)
        @test A ≠ view(A,Inclusion(0:0.5:0.5),:)
        @test view(A,Inclusion(0:0.5:0.5),:) ≠ A
        @test A ≠ view(B,Inclusion(0:0.5:1),:)
        @test view(B,Inclusion(0:0.5:1),:) ≠ A
        @test view(B,Inclusion(0:0.5:1),:) ≠ view(A,Inclusion(0:0.5:0.5),:)
    end
end


text:
using DifferentialEquations
using OrdinaryDiffEq
using MomentClosure
using ModelingToolkit
using DiffEqJump
using DataFrames
using Tables
using Statistics
using Plots
using StatsPlots
using Latexify
using BenchmarkTools

@parameters t β c γ
@variables S(t) I(t);

rxs = [Reaction(β*c, [S,I], [I], [1,1], [2])
       Reaction(γ, [I], nothing)]
rs  = ReactionSystem(rxs, t, [S,I], [β,c,γ]);

δt = 0.1
tmax = 40.0
tspan = (0.0,tmax)
ts = 0:δt:tmax;

u0 = [S => 990.0, I => 10.0];

u0v = [x[2] for x in u0];

p = [β=>0.00005, c=>10.0, γ=>0.25];

central_eqs = generate_central_moment_eqs(rs, 4, combinatoric_ratelaw=false);

closure_methods = ["zero","normal","log-normal","gamma","derivative matching"];

closed_central_eqs = Dict(cm=>moment_closure(central_eqs,cm) for cm in closure_methods);

u0map = Dict(cm=> deterministic_IC(u0v,closed_central_eqs[cm]) for cm in closure_methods);

closed_central_eqs_df = Dict{String,DataFrame}()
for cm in closure_methods
    prob = ODEProblem(closed_central_eqs[cm], u0map[cm], tspan, p)
    sol = solve(prob)
    df = DataFrame(sol(ts)')
    rename!(df,[replace(string(x[1]),"(t)" => "") for x in u0map[cm]])
    df[!,:t] = ts
    closed_central_eqs_df[cm] = df
end;

jumpsys = convert(JumpSystem, rs)
u0i = [S => 990, I => 10]
dprob = DiscreteProblem(jumpsys, u0i, tspan, p)
jprob = JumpProblem(jumpsys, dprob, Direct());

ensemble_jprob = EnsembleProblem(jprob)
ensemble_jsol = solve(ensemble_jprob,SSAStepper(),trajectories=10000)
ensemble_summary = EnsembleSummary(ensemble_jsol,ts);

ensemble_u = DataFrame(ensemble_summary.u')
rename!(ensemble_u,["μ₁₀","μ₀₁"])
ensemble_v = DataFrame(ensemble_summary.v')
rename!(ensemble_v,["M₂₀","M₀₂"])
ensemble_uv = hcat(ensemble_u,ensemble_v)
ensemble_uv[!,:t] = ts;

jplot = @df ensemble_uv plot(:t,[:μ₁₀,:μ₀₁],
     ribbon=[2*sqrt.(:M₂₀),
             2*sqrt.(:M₀₂)],
     label=["S" "I"],
     xlabel="Time",
     ylabel="Number",
     title="Jump process")

pltlist = []
for cm in closure_methods
     plt = @df closed_central_eqs_df[cm] plot(:t,[:μ₁₀,:μ₀₁],
          ribbon=[2*sqrt.(:M₂₀),
                  2*sqrt.(:M₀₂)],
          label=["S" "I"],
          xlabel="Time",
          ylabel="Number",
          title=cm)
    push!(pltlist,plt)
end;

l = @layout [a b c; d e f]
plot(vcat(jplot,pltlist)...,
     layout=l,
     legendfontsize=8,
     titlefontsize=10)

m = []
c = []
v = []
for moment in [:μ₁₀,:μ₀₁,:M₂₀,:M₀₂]
    for cm in closure_methods
        push!(m,moment)
        push!(c,cm)
        push!(v,mean(abs.(closed_central_eqs_df[cm][!,moment] - ensemble_uv[!,moment])))
    end
end
df = DataFrame(Dict(["Moment" => m,"Method" => c,"Normalized L1" => v]))
table = latexify(df,latex=false,env=:mdtable)
table

@benchmark solve(ensemble_jprob,SSAStepper(),trajectories=20000)

prob = ODEProblem(closed_central_eqs["normal"], u0map["normal"], tspan, p)
@benchmark sol = solve(prob)


text:
function test_grad_desc()

    # Vectorize
    for i in 1:10
        target = [10.0 * rand() + 5.0, 10.0 * rand() + 5.0]
        @show target
        x0 = [0.0, 0.0]
        rexps = rand(2:10, 2)
        @show rexps
        maxΔx = 1.0 ./ rexps
        minΔx = 0.0 ./ rexps
        x1 = x0 .+ maxΔx .* 0.1
        gdth = 1e-5
        maxiter = 50000

        function up_fun(gdmodel)
            xs = SimTools.gd_value(gdmodel)
            xs .^ rexps
        end

        gdmodel = SimTools.grad_desc_vec(up_fun;
            target, x0, x1, minΔx, maxΔx, gdth, maxiter,
        )
        @show gdmodel.ϵi
        @show gdmodel.iter
        @show up_fun(gdmodel)
        @test all(isapprox.(up_fun(gdmodel), target; atol = 1e-3))
        println()
    end

    # Scalar
    @testset begin
        for i in 1:10
            target = 10.0 * rand() + 5.0
            @show target
            x0 = 0.0
            rexp = rand(1:10)
            @show rexp
            maxΔx = 1.0 / rexp
            x1 = x0 + maxΔx * 0.1
            maxiter = 1e4
            gdth = 1e-5

            function up_fun(gdmodel)
                x = SimTools.gd_value(gdmodel)
                x ^ rexp
            end
            gdmodel = SimTools.grad_desc(up_fun;
                target, x0, x1, maxΔx, gdth, maxiter
            )
            @show gdmodel.ϵi
            @show gdmodel.iter
            @show up_fun(gdmodel)
            @test all(isapprox.(up_fun(gdmodel), target; atol = 1e-3))
            println()
        end
    end

end
test_grad_desc()


text:
println("loading mumps lib ...")
using MPI
using MUMPS # needs to be intialized before Ipopt don't understand why
println("mumps lib loaded.")

type linear_solver_MUMPS <: abstract_linear_system_solver
	_factor::MUMPS.Mumps{Float64}
	M::SparseMatrixCSC{Float64,Int64}
	sym::Symbol
	safe_mode::Bool

  function linear_solver_MUMPS(sym::Symbol, safe_mode::Bool)
      this = new();
			this.sym = sym
			this.safe_mode = safe_mode

			return this
  end
end

# intialize and finalize function???
# MPI ???
function mumps_sym(sym::Symbol)
	if sym == :symmetric
			return MUMPS.mumps_symmetric
	elseif sym == :unsymmetric
			return MUMPS.mumps_unsymmetric
	elseif sym == :definite
			return MUMPS.mumps_definite
	else
			error("this symmetry symbol is not understood")
	end
end

function create_mumps_factor(solver::linear_solver_MUMPS)
	cntl = MUMPS.default_cntl64[:]; # copy
	cntl[1] = 0.01

	icntl = MUMPS.default_icntl[:]; # copy

	icntl[4] = 1;
	icntl[10] = 1; # 2 iterative refinement steps
	icntl[14] = 200.0 #1000.0;

	return MUMPS.Mumps{Float64}(mumps_sym(solver.sym), icntl, cntl	);  # Real, general unsymmetric
end

function create_mumps_factor2(solver::linear_solver_MUMPS)
	icntl = get_icntl(verbose=false);
	return MUMPS.Mumps{Float64}(mumps_sym(solver.sym), icntl, MUMPS.default_cntl64);  # Real, general unsymmetric
end

function initialize!(solver::linear_solver_MUMPS)
		if ~MPI.Initialized()
			MPI.Init()
		end

    solver._factor = create_mumps_factor2(solver)
end

function finalize!(solver::linear_solver_MUMPS)
    MUMPS.finalize(solver._factor);
end

function ls_factor!(solver::linear_solver_MUMPS, SparseMatrix::SparseMatrixCSC{Float64,Int64}, n::Int64, m::Int64, timer::class_advanced_timer)
		#if
		return ls_factor_LBL!(solver, SparseMatrix, n, m, timer)
end

# A-matrix
function ls_factor_LBL!(solver::linear_solver_MUMPS, SparseMatrix::SparseMatrixCSC{Float64,Int64}, n::Int64, m::Int64, timer::class_advanced_timer)
		start_advanced_timer(timer,"MUMPS")
		@assert(size(SparseMatrix,1) == n + m)
		@assert(size(SparseMatrix,2) == n + m)

		solver.M = deepcopy(SparseMatrix)

		start_advanced_timer(timer,"MUMPS/associate_matrix")
    	associate_matrix!(solver._factor, solver.M);
		pause_advanced_timer(timer,"MUMPS/associate_matrix")

		start_advanced_timer(timer,"MUMPS/factorize")
		factorize!(solver._factor);
		pause_advanced_timer(timer,"MUMPS/factorize")

		zero_eigs = solver._factor.infog[1] == -10 || solver._factor.infog[1] == -2

		pause_advanced_timer(timer,"MUMPS")

		if solver._factor.infog[12] > m
				return -1
		end

		if zero_eigs
				return 0
		end

		return 1
end

function ls_solve(solver::linear_solver_MUMPS, my_rhs::AbstractArray, timer::class_advanced_timer)
		start_advanced_timer(timer,"MUMPS")
		start_advanced_timer(timer,"MUMPS/solve")

		if solver.safe_mode
			MUMPS.associate_matrix!(solver._factor, solver.M);
			MUMPS.factorize!(solver._factor);
		end

		associate_rhs!(solver._factor, reshape(my_rhs,length(my_rhs),1))
		solve!(solver._factor)
		sol = deepcopy(get_solution(solver._factor))
		#sol = solve(solver._factor, my_rhs)

		pause_advanced_timer(timer,"MUMPS/solve")
		pause_advanced_timer(timer,"MUMPS")

    return sol[:]
end


text:
using DelaySSAToolkit, Catalyst
using LaTeXStrings
using Plots
theme(:vibrant)

rn = @reaction_network begin
    C, 0 --> Xₐ
    γ, Xₐ --> 0
    β, Xₐ --> Xᵢ₁ + Xᵢ₂
    γ, Xᵢ₁ --> 0
    γ, Xᵢ₂ --> 0
 end C γ β
 jumpsys = convert(JumpSystem, rn, combinatoric_ratelaws = false)

C, γ, β, τ = [2., 0.1, 0.5, 15.]
p=[C, γ, β]
u0 = [0, 0, 0]
tf = 30.
saveat = .1
de_chan0 = [[],[]]
tspan = (0.,tf)
dprob = DiscreteProblem(u0, tspan,p)

# delay_trigger_affect! = function (integrator, rng)
#    append!(integrator.de_chan[1], τ)
#    append!(integrator.de_chan[2], τ)
# end
# delay_trigger = Dict(3=>delay_trigger_affect!)
delay_trigger = Dict(3=>[1=>τ,2=>τ])
delay_complete = Dict(1=>[2=>-1],2=>[3=>-1])

delay_affect1! = function (integrator, rng)
    i = rand(rng, 1:length(integrator.de_chan[1]))
    deleteat!(integrator.de_chan[1],i)
end
delay_affect2! = function (integrator, rng)
    i = rand(rng, 1:length(integrator.de_chan[2]))
    deleteat!(integrator.de_chan[2],i)
end
delay_interrupt = Dict(4=>delay_affect1!,5=>delay_affect2!)
delayjumpset = DelayJumpSet(delay_trigger,delay_complete,delay_interrupt)

djprob = DelayJumpProblem(jumpsys,dprob, DelayRejection(),  delayjumpset, de_chan0, save_positions=(false,false))
sol =@time solve(djprob, SSAStepper(),seed=10, saveat =.1)
plot(sol, label = [L"X_A" L"X_{I1}" L"X_{I2}"],legend =:topleft,xlabel="time",ylabel="Molecular Number")

timestamps = 0:0.1:tf
a = β + γ
x_A(t) = C/a*(1-exp(-a*t))
x_I(t)= 0<=t<=τ ? C*β/(a-γ)*((1-exp(-γ*t))/γ - (1-exp(-a*t))/a) : C*β/a*((1-exp(-γ*τ))/γ + exp(-a*t)*(1-exp((a-γ)τ))/(a-γ))
using Plots
plot(timestamps,x_A.(timestamps),linewidth=3,label=L"X_A",xlabel="time",ylabel="Mean Value")
plot!(timestamps,x_I.(timestamps),linewidth=3,label=L"X_I")

using StatsBase, DifferentialEquations.EnsembleAnalysis
Sample_size = Int(1e4)
ens_prob = EnsembleProblem(djprob)
ens =@time solve(ens_prob, SSAStepper(),EnsembleThreads(),trajectories = Sample_size, saveat = .1)

tsm(t) = timepoint_mean(ens,t)

mean_A(t) = tsm(t)[1]
mean_I1(t) = tsm(t)[2]
mean_I2(t) = tsm(t)[3]

plot(timestamps,x_A.(timestamps),linewidth=3,label=L"$X_A$ Exact",xlabel="time",ylabel="Mean Value")
plot!(timestamps,mean_A.(timestamps),linewidth=3,line=:dash,label=L"$X_A$ Delay SSA")
plot!(timestamps,x_I.(timestamps),linewidth=3,label=L"$X_I$ Exact")
plot!(timestamps,mean_I1.(timestamps),linewidth=3,line=:dash,label=L"$X_{I1}$ Delay SSA")
plot!(timestamps,mean_I2.(timestamps),linewidth=3,line=:dash, legend =:topleft,label=L"$X_{I2}$ Delay SSA")


text:
using ITensors
using ITensorChemistry
using Test

@testset "ITensorChemistry.jl" begin
  molecule = "H₂"
  basis = "sto-3g"
  (; hamiltonian, state, hartree_fock_energy) = molecular_orbital_hamiltonian_and_state(; molecule, basis)

  s = siteinds("Electron", length(state); conserve_qns=true)
  H = MPO(hamiltonian, s)
  ψhf = MPS(s, state)

  @test inner(ψhf, H, ψhf) ≈ hartree_fock_energy

  sweeps = Sweeps(10)
  setmaxdim!(sweeps, 100, 200)
  setcutoff!(sweeps, 1e-6)
  setnoise!(sweeps, 1e-6, 1e-7, 1e-8, 0.0)
  dmrg_energy, ψ = dmrg(H, ψhf, sweeps; output_level=0)

  @test dmrg_energy < hartree_fock_energy
end


text:
using Documenter, TextSearch

makedocs(;
    modules=[TextSearch],
    authors="Eric S. Tellez",
    repo="https://github.com/sadit/TextSearch.jl/blob/{commit}{path}#L{line}",
    sitename="TextSearch.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", nothing) == "true",
        canonical="https://sadit.github.io/TextSearch.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "API" => "api.md"
    ],
)

deploydocs(;
    repo="github.com/sadit/TextSearch.jl",
    devbranch=nothing,
    branch = "gh-pages",
    versions = ["stable" => "v^", "v#.#.#", "dev" => "dev"]
)

max_stars_repo_name: sadit/TextModel

max_issues_repo_name: sadit/TextModel

max_forks_repo_name: sadit/TextModel


text:
using NIfTI, CodecZlib, TranscodingStreams
using Test

function extractto(gzname, out)
    open(out, "w") do io_out
        open(gzname, "r") do io_in
            gz = GzipDecompressorStream(io_in)
            write(io_out, read(gz))
        end
    end
end

const TEMP_DIR_NAME = mktempdir()

# single file storage
const GZIPPED_NII = joinpath(dirname(@__FILE__), "data/example4d.nii.gz")
const NII = joinpath(TEMP_DIR_NAME, "$(tempname()).nii")
extractto(GZIPPED_NII, NII)

# dual file storage
const GZIPPED_HDR = joinpath(dirname(@__FILE__), "data/example4d.hdr.gz")
hdr_stem = joinpath(TEMP_DIR_NAME, tempname())
const HDR = "$hdr_stem.hdr"
const IMG = "$hdr_stem.img"
extractto(GZIPPED_HDR, HDR)
extractto(joinpath(dirname(@__FILE__), "data/example4d.img.gz"), IMG)

for (fname, mmap) in ((NII, false), (NII, true), (HDR, false), (HDR, true),
                    (GZIPPED_NII, false), (GZIPPED_HDR, false))
    file = niread(fname, mmap=mmap)

    # Header
    @test time_step(file.header) == 2000000 # Actually an error in the file AFAIK
    @test voxel_size(file.header) ≈ Float32[2.0, 2.0, 2.2]
    @test size(file) == (128, 96, 24, 2)

    # Content
    @test file.raw[65, 49, 13, :][:] == [265, 266]
    @test vox(file, 64, 48, 12, :)[:] == [265, 266]
    @test vox(file, 69, 56, 13, :)[:] == [502, 521]

    @assert maximum(file) == maximum(file.raw)
end

@test_throws ErrorException niread(GZIPPED_NII; mmap=true)
@test_throws ErrorException niread(GZIPPED_HDR; mmap=true)

# Test writing
const TEMP_FILE = joinpath(TEMP_DIR_NAME, "$(tempname()).nii")
vol = NIVolume()
niwrite(TEMP_FILE, vol)
niread(TEMP_FILE)

const TEMP_GZIPPED_FILE = joinpath(TEMP_DIR_NAME, "$(tempname()).nii.gz")
niwrite(TEMP_GZIPPED_FILE, vol)
niread(TEMP_GZIPPED_FILE)

# Write and read DT_BINARY
const BOOL_WRITE = joinpath(TEMP_DIR_NAME, "$(tempname()).nii")
const BIT_WRITE = joinpath(TEMP_DIR_NAME, "$(tempname()).nii")
mask = rand(Bool, 3, 5, 7) # Array{Bool}
mask_bitarray = BitArray(mask) # BitArray
niwrite(BOOL_WRITE, NIVolume(mask))
niwrite(BIT_WRITE, NIVolume(mask_bitarray))
@test niread(BOOL_WRITE).raw == mask
@test niread(BIT_WRITE).raw == mask_bitarray

# Open mmaped file for reading and writing
const WRITE = joinpath(TEMP_DIR_NAME, "$(tempname()).nii")
const VERIFY_WRITE = joinpath(TEMP_DIR_NAME, "$(tempname()).nii")
cp(NII, WRITE)
img = niread(WRITE; mmap=true, mode="r+")
img.raw[1,1,1,1] = 5
img.raw[:,2,1,1] = ones(size(img)[1])
cp(WRITE, VERIFY_WRITE)
@test niread(VERIFY_WRITE)[1,1,1,1] == 5
@test niread(VERIFY_WRITE)[:,2,1,1] == ones(size(img)[1])
# Site is currently down TODO: reintroduce this test when site is up
# Big endian
# const BE = "$(tempname()).nii"
# download("https://nifti.nimh.nih.gov/nifti-1/data/avg152T1_LR_nifti.nii.gz", BE)
img = niread(joinpath(dirname(@__FILE__), "data/avg152T1_LR_nifti.nii.gz"))
@test size(img) == (91,109,91)

GC.gc() # closes mmapped files


text:
# --- Coupling Model Creation --- #

"""
    Coupling(::LiftingLine, ::GXBeamAssembly; kwargs...)

Coupling model for coupling a lifting line aerodynamic model and a geometrically exact beam
element assembly model.  This model introduces additional parameters corresponding to
the external point loads ``F_{x,i}, F_{y,i}, F_{z,i}, M_{x,i}, M_{y,i}, M_{z,i}`` or
displacements ``u_{x,i}, u_{y,i}, u_{z,i}, \\theta_{x,i}, \\theta_{y,i}, \\theta_{z,i}``
applied to each node, followed by the distributed loads ``f_{x,i}, f_{y,i}, f_{z,i},
m_{x,i}, m_{y,i}, m_{z,i}`` applied to each beam element (excluding aerodynamic loads),
followed by the properties of point masses attached to each beam element ``m, p, I_{11},
I_{22}, I_{33}, I_{12}, I_{13}, I_{23}``, followed by the gravity vector, linear velocity,
angular velocity, linear acceleration, and angular acceleration of the system, followed by
the freestream air density ``\\rho_\\infty`` and air speed of sound ``c``.

# Keyword Arguments
 - `lifting_elements`: Beam element index corresponding to each aerodynamic section.
    By default, aerodynamic models are assigned to beam elements in order

**NOTE: This model assumes that each beam element is oriented with the x-axis along the
beam's axis, the y-axis forward (into the freestream), and the z-axis normal**
"""
function Coupling(aero::LiftingLine, stru::GXBeamAssembly;
    lifting_elements = 1:length(aero.submodels))

    # section models
    section_models = aero.section_models

    # number of points and elements
    npoint = length(gxbeam.constants.icol_point)
    nelem = length(gxbeam.constants.icol_elem)

    # displacement constraints
    displacement = gxbeam.constants.displacement

    # number of additional parameters
    npc = 6*npoint + 6*nelem + 10*nelem + 17

    # dimensions
    nx = number_of_states(liftingline) + number_of_states(gxbeam)
    ny = number_of_inputs(liftingline) + number_of_inputs(gxbeam)
    np = number_of_parameters(liftingline) + number_of_parameters(gxbeam) + npc

    # coupling function
    g = (y, dx, x, p, t) -> liftingline_gxbeam_inputs!(y, dx, x, p, t;
        aero, stru, lifting_elements, section_models)

    # jacobians
    ratejac = Nonlinear() # TODO: define rate jacobian
    statejac = Nonlinear() # TODO: define state jacobian
    paramjac = Nonlinear() # TODO: define parameter jacobian
    tgrad = Nonlinear() # TODO: define time gradient

    # input/output functions
    setparam = (p; kwargs...) -> liftingline_gxbeam_setparam!(p, displacement,
        npoint, nelem; kwargs...)
    sepparam = (p) -> liftingline_gxbeam_sepparam(p, displacement, npoint, nelem)

    return Coupling{true}(g, nx, ny, np, npc;
        ratejac = ratejac,
        statejac = statejac,
        paramjac = paramjac,
        tgrad = tgrad,
        setparam = setparam,
        sepparam = sepparam)
end

function liftingline_gxbeam_inputs!(y, dx, x, p, t;
    liftingline_submodel, gxbeam_submodel, lifting_elements, section_models)

    # get indices of states, inputs, and parameters for each submodel
    submodels = (liftingline, gxbeam)
    ixa, ixs = state_indices(submodels)
    iya, iys = input_indices(submodels)
    ipa, ips = parameter_indices(submodels)

    # get views of states, inputs, and parameters for each submodel
    dxa, dxs = view(dx, ixa), view(dx, ixs)
    xa, xs = view(x, ixa), view(x, ixs)
    ya, ys = view(y, iya), view(y, iys)
    pa, ps = view(p, ipa), view(p, ips)

    # get indices of states, inputs, and parameters for each lifting line section model
    ixas = state_indices(section_models)
    iyas = input_indices(section_models)
    ipas = parameter_indices(section_models)

    # get views of states, inputs, and parameters for each lifting line section model
    dxas = view.(Ref(dxa), ixas)
    xas = view.(Ref(xa), ixas)
    yas = view.(Ref(ya), iyas)
    pas = view.(Ref(pa), ipas)

    # separate gxbeam model rates, states, and parameters
    gxbeam_rates = separate_rates(gxbeam, dxs)
    gxbeam_states = separate_states(gxbeam, xs)
    gxbeam_parameters = separate_parameters(gxbeam, ps)

    # separate additional parameters
    np = length(p)
    npc = np - sum(number_of_parameters.(submodels))
    ipc = np - npc : np
    pc = view(p, ipc)
    additional_parameters = liftingline_gxbeam_sepparam(pc, gxbeam.constants.displacement,
        length(gxbeam.constants.icol_point), length(gxbeam.constants.icol_elem))

    # initialize distributed loads
    f_d = SVector{3,eltype(y)}.(additional_parameters.f_d)
    m_d = SVector{3,eltype(y)}.(additional_parameters.m_d)

    # calculate aerodynamic inputs and loads
    for (i, ielem) in enumerate(lifting_elements)

        # unpack element properties
        elem = gxbeam_parameters.assembly.elements[ielem]

        # transformation from body to element frame
        C1 = gxbeam_body_to_element(elem, gxbeam_states.theta_e[ielem])

        # element velocities and accelerations (in the element frame)
        vi, ωi = gxbeam_element_velocities(gxbeam_states.V_e[ielem], gxbeam_states.Omega_e[ielem])
        ai, αi = gxbeam_element_accelerations(elem, gxbeam_states.u_e[ielem],
            gxbeam_states.theta_e[ielem], gxbeam_rates.V_e[ielem], gxbeam_rates.Omega_e[ielem],
            additional_parameters.a, additional_parameters.alpha)

        # transformation from element to local aerodynamic frame
        C2 = @SMatrix [0 -1 0; 1 0 0; 0 0 1]

        # freestream velocities and accelerations (in the local aerodynamic frame)
        vi, ωi = -C2*vi, C2*ωi
        ai, αi = -C2*ai, C2*αi

        # calculate lifting line inputs
        yi = liftingline_section_inputs(section_models[i], dxas[i], xas[i], pas[i],
            vi, ωi, ai, αi, additional_parameters.rho, additional_parameters.c, t)

        # separate section aerodynamic and structural inputs
        aerodynamic_inputs, structural_inputs = separate_inputs(section_models[i], yi)

        # save section aerodynamic inputs
        set_inputs!(yas[i], liftingline_section_models[i]; aerodynamic_inputs...)

        # transformation from local aerodynamic frame to body frame
        C3 = C1'*C2'

        # add aerodynamic loads to prescribed distributed loads
        f_d[ielem] += C3*structural_inputs.f
        m_d[ielem] += C3*structural_inputs.m
    end

    # save gxbeam inputs
    set_inputs!(ys, gxbeam;
        u_p = additional_parameters.u_p,
        theta_p = additional_parameters.theta_p,
        F_p = additional_parameters.F_p,
        M_p = additional_parameters.M_p,
        f_d = f_d,
        m_d = m_d,
        m_m = additional_parameters.m_m,
        p_m = additional_parameters.p_m,
        I_m = additional_parameters.I_m,
        gvec = additional_parameters.gvec,
        v = additional_parameters.v,
        omega = additional_parameters.omega,
        a = additional_parameters.a,
        alpha = additional_parameters.alpha)

    return y
end

# convenience function for defining this model's input vector
function liftingline_gxbeam_setparam!(p, displacement, np, ne;
    u_p = nothing, theta_p = nothing, F_p = nothing, M_p = nothing,
    f_d = nothing, m_d = nothing,
    m_m = nothing, p_m = nothing, I_m = nothing,
    gvec = nothing, v = nothing, omega = nothing, a = nothing, alpha = nothing,
    rho = nothing, c = nothing)

    for ip = 1:np
        prescribed_forces = SVector{6}(view(displacement, :, ip)) .== false
        for i = 1:3
            if prescribed_forces[i]
                if !isnothing(F_p)
                    p[6*(ip-1) + i] = F_p[ip][i]
                end
            else
                if !isnothing(u_p)
                    p[6*(ip-1) + i] = u_p[ip][i]
                end
            end
        end
        for i = 1:3
            if prescribed_forces[i]
                if !isnothing(M_p)
                    p[6*(ip-1) + 3 + i] = M_p[ip][i]
                end
            else
                if !isnothing(theta_p)
                    p[6*(ip-1) + 3 + i] = theta_p[ip][i]
                end
            end
        end
    end

    for ie = 1:ne
        if !isnothing(f_d)
            for i = 1:3
                p[6*np + 6*(ie-1) + i] = f_d[ie][i]
            end
        end
        if !isnothing(m_d)
            for i = 1:3
                p[6*np + 6*(ie-1) + 3 + i] = m_d[ie][i]
            end
        end
    end

    for ie = 1:ne
        if !isnothing(m_m)
            p[6*np + 6*ne + 10*(ie-1) + 1] = m_m[ie]
        end
        if !isnothing(p_m)
            p[6*np + 6*ne + 10*(ie-1) + 2] = p_m[ie][1]
            p[6*np + 6*ne + 10*(ie-1) + 3] = p_m[ie][2]
            p[6*np + 6*ne + 10*(ie-1) + 4] = p_m[ie][3]
        end
        if !isnothing(I_m)
            p[6*np + 6*ne + 10*(ie-1) + 5] = I_m[ie][1,1]
            p[6*np + 6*ne + 10*(ie-1) + 6] = I_m[ie][2,2]
            p[6*np + 6*ne + 10*(ie-1) + 7] = I_m[ie][3,3]
            p[6*np + 6*ne + 10*(ie-1) + 8] = I_m[ie][1,2]
            p[6*np + 6*ne + 10*(ie-1) + 9] = I_m[ie][1,3]
            p[6*np + 6*ne + 10*(ie-1) + 10] = I_m[ie][2,3]
        end
    end

    if !isnothing(gvec)
        p[6*np + 6*ne + 10*ne + 1] = gvec[1]
        p[6*np + 6*ne + 10*ne + 2] = gvec[2]
        p[6*np + 6*ne + 10*ne + 3] = gvec[3]
    end

    if !isnothing(v)
        p[6*np + 6*ne + 10*ne + 4] = v[1]
        p[6*np + 6*ne + 10*ne + 5] = v[2]
        p[6*np + 6*ne + 10*ne + 6] = v[3]
    end

    if !isnothing(omega)
        p[6*np + 6*ne + 10*ne + 7] = omega[1]
        p[6*np + 6*ne + 10*ne + 8] = omega[2]
        p[6*np + 6*ne + 10*ne + 9] = omega[3]
    end

    if !isnothing(a)
        p[6*np + 6*ne + 10*ne + 10] = a[1]
        p[6*np + 6*ne + 10*ne + 11] = a[2]
        p[6*np + 6*ne + 10*ne + 12] = a[3]
    end

    if !isnothing(alpha)
        p[6*np + 6*ne + 10*ne + 13] = alpha[1]
        p[6*np + 6*ne + 10*ne + 14] = alpha[2]
        p[6*np + 6*ne + 10*ne + 15] = alpha[3]
    end

    if !isnothing(rho)
        p[6*np + 6*ne + 10*ne + 16] = rho
    end

    if !isnothing(c)
        p[6*np + 6*ne + 10*ne + 17] = c
    end

    return p
end

# convenience function for separating this model's input vector
function liftingline_gxbeam_sepparam(p, displacement, np, ne)

    TF = eltype(p)

    # separate prescribed conditions
    u_p = Vector{SVector{3,TF}}(undef, np)
    theta_p = Vector{SVector{3,TF}}(undef, np)
    F_p = Vector{SVector{3,TF}}(undef, np)
    M_p = Vector{SVector{3,TF}}(undef, np)
    for ip = 1:np
        prescribed_forces = SVector{6}(view(displacement, :, ip)) .== false
        F_p[ip] = SVector(
            ifelse(prescribed_forces[1], p[6*(ip-1) + 1], NaN),
            ifelse(prescribed_forces[2], p[6*(ip-1) + 2], NaN),
            ifelse(prescribed_forces[3], p[6*(ip-1) + 3], NaN))
        M_p[ip] = SVector(
            ifelse(prescribed_forces[4], p[6*(ip-1) + 4], NaN),
            ifelse(prescribed_forces[5], p[6*(ip-1) + 5], NaN),
            ifelse(prescribed_forces[6], p[6*(ip-1) + 6], NaN))
        u_p[ip] = SVector(
            ifelse(prescribed_forces[1], NaN, p[6*(ip-1) + 1]),
            ifelse(prescribed_forces[2], NaN, p[6*(ip-1) + 2]),
            ifelse(prescribed_forces[3], NaN, p[6*(ip-1) + 3]))
        theta_p[ip] = SVector(
            ifelse(prescribed_forces[4], NaN, p[6*(ip-1) + 4]),
            ifelse(prescribed_forces[5], NaN, p[6*(ip-1) + 5]),
            ifelse(prescribed_forces[6], NaN, p[6*(ip-1) + 6]))
    end

    # separate distributed loads
    f_d = Vector{SVector{3,TF}}(undef, ne)
    m_d = Vector{SVector{3,TF}}(undef, ne)
    for ie = 1:ne
        f_d[ie] = SVector(
            p[6*np + 6*(ie-1) + 1],
            p[6*np + 6*(ie-1) + 2],
            p[6*np + 6*(ie-1) + 3])
        m_d[ie] = SVector(
            p[6*np + 6*(ie-1) + 4],
            p[6*np + 6*(ie-1) + 5],
            p[6*np + 6*(ie-1) + 6])
    end

    # separate point masses
    m_m = Vector{TF}(undef, ne)
    p_m = Vector{SVector{3,TF}}(undef, ne)
    I_m = Vector{SMatrix{3,3,TF}}(undef, ne)
    for ie = 1:ne
        m_m[ie] = p[6*np + 6*ne + 10*(ie-1) + 1]
        p_m[ie] = SVector(
            p[6*np + 6*ne + 10*(ie-1) + 2],
            p[6*np + 6*ne + 10*(ie-1) + 3],
            p[6*np + 6*ne + 10*(ie-1) + 4])
        I11 = p[6*np + 6*ne + 10*(ie-1) + 5]
        I22 = p[6*np + 6*ne + 10*(ie-1) + 6]
        I33 = p[6*np + 6*ne + 10*(ie-1) + 7]
        I12 = p[6*np + 6*ne + 10*(ie-1) + 8]
        I13 = p[6*np + 6*ne + 10*(ie-1) + 9]
        I23 = p[6*np + 6*ne + 10*(ie-1) + 10]
        I_m[ie] = @SMatrix [I11 I12 I13; I12 I22 I23; I13 I23 I33]
    end

    # gravity
    gvec =  view(p, 6*np + 6*ne + 10*ne +  1 : 6*np + 6*ne + 10*ne +  3)

    # linear/angular velocity
    v =     view(p, 6*np + 6*ne + 10*ne +  4 : 6*np + 6*ne + 10*ne +  6)
    omega = view(p, 6*np + 6*ne + 10*ne +  7 : 6*np + 6*ne + 10*ne +  9)

    # linear/angular acceleration
    a =     view(p, 6*np + 6*ne + 10*ne + 10 : 6*np + 6*ne + 10*ne + 12)
    alpha = view(p, 6*np + 6*ne + 10*ne + 13 : 6*np + 6*ne + 10*ne + 15)

    # freestream air density
    rho = p[6*np + 6*ne + 10*ne + 16]

    # air speed of sound
    c = p[6*np + 6*ne + 10*ne + 17]

    return (u_p=u_p, theta_p=theta_p, F_p=F_p, M_p=M_p, f_d=f_d, m_d=m_d, m_m=m_m, p_m=p_m,
        I_m=I_m, gvec=gvec, v=v, omega=omega, a=a, alpha=alpha, rho=rho, c=c)
end


text:
@testset "DictTabular" begin
    @testset "DictSeries" begin
        @test @inferred(Series(Dict(:a => 1, :b => 2))) isa Series
        s = Series(Dict(:a => 1, :b => 2))
        @test @inferred(indices(s)) == (keys(get(s)),)
        @test @inferred(s[:a]) === 1
        @test @inferred(s[[:a, :b]]) == s
        @test @inferred(s[:]) == s

        @test (s[:a] = 5; s[:a] === 5)
    end

    @testset "DictTable" begin
        # DictTable with inner DictSeries
        @test @inferred(Table(Dict(:a => Dict(:x=>1,:y=>2,:z=>3), :b => Dict(:x=>4,:y=>5,:z=>6)))) isa Table
        t = Table(Dict(:a => Dict(:x=>1,:y=>2,:z=>3), :b => Dict(:x=>4,:y=>5,:z=>6)))
        @test @inferred(indices(t)) == (keys(get(t)[:a]), keys(get(t)))
        @test @inferred(t[:y, :a]) === 2
        @test @inferred(t[:y, [:a, :b]]) == Series(Dict(:a => 2, :b => 5))
        #if VERSION <= v"0.6.0-rc3.0"
        #    # Inference works at REPL but not in test
        #    @test t[:y, :] == Series(Dict(:a => 2, :b => 5))
        #    @test_broken @inferred(t[:y, :]) == Series(Dict(:a => 2, :b => 5))
        #else
            @test @inferred(t[:y, :]) == Series(Dict(:a => 2, :b => 5))
        #end
        @test @inferred(t[[:x,:y,:z], :a]) == Series(Dict(:x => 1, :y => 2, :z => 3))
        @test @inferred(t[[:x,:y,:z], [:a, :b]]) == t
        if VERSION <= v"0.6.0-rc3.0"
            # Inference works at REPL but not in test
            @test t[[:x,:y,:z], :] == t
            @test_broken @inferred(t[[:x,:y,:z], :]) == t
        else
            @test @inferred(t[[:x,:y,:z], :]) == t
        end
        @test @inferred(t[:, :a]) == Series(Dict(:x => 1, :y => 2, :z => 3))
        @test @inferred(t[:, [:a, :b]]) == t
        #if VERSION <= v"0.6.0-rc3.0"
        #    # Inference works at REPL but not in test
        #    @test t[:, :] == t
        #    @test_broken @inferred(t[:, :]) == t
        #else
            @test @inferred(t[:, :]) == t
        #end

        @test (t[:y, :a] = 5; t[:y, :a] === 5)

        # DictTable with inner ArraySeries
        @test @inferred(Table(Dict(:a => [1,2,3], :b => [4,5,6]))) isa Table
        t2 = Table(Dict(:a => [1,2,3], :b => [4,5,6]))
        @test @inferred(indices(t2)) == (Base.OneTo(3), keys(get(t2)))
        @test @inferred(t2[2, :a]) === 2
        @test @inferred(t2[2, [:a, :b]]) == Series(Dict(:a => 2, :b => 5))
        #if VERSION <= v"0.6.0-rc3.0"
        #    # Inference works at REPL but not in test
        #    @test t2[2, :] == Series(Dict(:a => 2, :b => 5))
        #    @test_broken @inferred(t2[2, :]) == Series(Dict(:a => 2, :b => 5))
        #else
            @test @inferred(t2[2, :]) == Series(Dict(:a => 2, :b => 5))
        #end
        @test @inferred(t2[:, :a]) == Series([1,2,3])
        @test @inferred(t2[:, [:a, :b]]) == t2
        #if VERSION <= v"0.6.0-rc3.0"
        #    # Inference works at REPL but not in test
        #    @test t2[:, :] == t2
        #    @test_broken @inferred(t2[:, :]) == t2
        #else
            @test @inferred(t2[:, :]) == t2
        #end

        @test (t2[2, :a] = 5; t2[2, :a] === 5)

        # DictTable with inner TupleSeries
        @test @inferred(Table(Dict(:a => (l"x"=>1,l"y"=>2,l"z"=>3), :b => (l"x"=>4,l"y"=>5,l"z"=>6)))) isa Table
        t3 = Table(Dict(:a => (l"x"=>1,l"y"=>2,l"z"=>3), :b => (l"x"=>4,l"y"=>5,l"z"=>6)))
        @test @inferred(indices(t3)) == ((l"x", l"y", l"z"), keys(get(t3)))
        @test @inferred(t3[l"y", :a]) === 2
        if VERSION <= v"0.6.0-rc3.0"
            # Inference works at REPL but not in test
            @test t3[l"y", :] == Series(Dict(:a => 2, :b => 5))
            @test_broken @inferred(t3[l"y", :]) == Series(Dict(:a => 2, :b => 5))
        else
            @test @inferred(t3[l"y", :]) == Series(Dict(:a => 2, :b => 5))
        end
        @test @inferred(t3[:, :a]) === Series(l"x" => 1, l"y" => 2, l"z" => 3)
        #if VERSION <= v"0.6.0-rc3.0"
        #    # Inference works at REPL but not in test
        #    @test t3[:, :] == t3
        #    @test_broken @inferred(t3[:, :]) == t3
        #else
            @test @inferred(t3[:, :]) == t3
        #end

        # DictTable with inner StructSeries
        @test @inferred(Table(Dict(:a => 1+2im, :b => 3+4im))) isa Table
        t4 = Table(Dict(:a => 1+2im, :b => 3+4im))
        @test @inferred(indices(t4)) == ((l"re", l"im"), keys(get(t4)))
        @test @inferred(t4[l"im", :a]) === 2
        if VERSION <= v"0.6.0-rc3.0"
            # Inference works at REPL but not in test
            @test t4[l"im", :] == Series(:a => 2, :b => 4)
            @test_broken @inferred(t4[l"im", :]) == Series(:a => 2, :b => 4)
        else
            @test @inferred(t4[l"im", :]) == Series(:a => 2, :b => 4)
        end
        @test @inferred(t4[:, :a]) == Series(1+2im)
        if VERSION <= v"0.6.0-rc3.0"
            # Inference works at REPL but not in test
            @test t4[:, :] == t4
            @test_broken @inferred(t4[:, :]) == t4
        else
            @test @inferred(t4[:, :]) == t4

        end
    end
end


text:
# NOTE Pkgs not intended to be watched by {Revise}, use them before {Revise} section

# General ----------------------------------------

# TODO Which one goes to {PackageCompiler}

# {OhMyREPL} -------------------------------------

using Crayons
using OhMyREPL  # TODO Not working with sysimg?
import OhMyREPL: Passes.SyntaxHighlighter

OhMyREPL.input_prompt!("❯ ", :black)
OhMyREPL.enable_pass!("RainbowBrackets", false)
OhMyREPL.enable_pass!("BracketHighlighter", false)

# TODO
# Markdown syntax
# OhMyREPL.Passes.BracketHighlighter.setcrayon!(
#     Crayon(foreground = :)
# )

scheme = SyntaxHighlighter.ColorScheme()

SyntaxHighlighter.keyword!(scheme,      crayon"light_cyan")
SyntaxHighlighter.function_def!(scheme, crayon"default")  # TODO
SyntaxHighlighter.macro!(scheme,        crayon"light_blue")  # TODO `@`
SyntaxHighlighter.text!(scheme,         crayon"default")  # TOO

SyntaxHighlighter.symbol!(scheme,       crayon"default")
SyntaxHighlighter.comment!(scheme,      crayon"dark_gray")
SyntaxHighlighter.call!(scheme,         crayon"light_blue")

SyntaxHighlighter.op!(scheme,           crayon"light_gray")
# TODO Brackets

SyntaxHighlighter.error!(scheme,        crayon"yellow")
SyntaxHighlighter.argdef!(scheme,       crayon"cyan")

SyntaxHighlighter.number!(scheme,       crayon"yellow")
SyntaxHighlighter.string!(scheme,       crayon"yellow")

# test_colorscheme(scheme)

SyntaxHighlighter.add!("elegant_dark", scheme)
colorscheme!("elegant_dark")

# Colors -----------------------------------------

ENV["JULIA_ERROR_COLOR"] = :red  # TODO Why this color works?
# ENV["JULIA_WARN_COLOR"] = :yellow
# ENV["JULIA_INFO_COLOR"] = :yellow

# TODO Not working
# function customize_colors(repl)
#     repl.help_color = Base.text_colors[:yellow]
# end
# atreplinit(customize_colors)

# REPL: Keybind, prompt --------------------------

import REPL
import REPL.LineEdit

# TODO Simplify with hotkey
function REPL.ends_with_semicolon(line::AbstractString)
    match = findlast(isequal(';'), line)::Union{Nothing,Int}

    # --------------------------------------------
    # CUSTOM Avoid printing for object assignment by testing if the first `=` is inside round brackets `(`
    # --------------------------------------------

    # Test: Manual or assign these to 'line' variable
    # @model Linear(x, y) = begin ; α ~ Normal() ; end
    # function f() ; a = 1 ; end
    # f(x, α, β = 1) = α .+ β*x  # TODO
    # a = 1
    # x = DataFrame(a = 1)
    # DataFrame(a = 1)
    # rand()
    # a == 1
    # a === 1

    # TODO `==` and `(` in comments
    # # a == 1
    # a = 1

    # # a == (1)
    # a = 1

    # Test for e.g. `DataFrame(a = 1)`
    first_equal_sign_idx = findfirst(isequal('='), line)::Union{Nothing, Int}

    if first_equal_sign_idx == nothing ; return false
    else
        # Test for function, macro & struct
        is__func_macro_struct =
            occursin(r"(\bfunction|@|struct|mutable struct\b)", line)

        str_to_first_equal_sign = line[1 : first_equal_sign_idx]
        has__open_bracket = occursin("(", str_to_first_equal_sign)

        if is__func_macro_struct ; return true

        elseif has__open_bracket
            has__close_bracket_followed_by_equal_sign = occursin(r"\)\s?=", line)

            # Test for e.g. `DataFrame(a = 1)`
            if ! has__close_bracket_followed_by_equal_sign ; return false
            # Test for e.g. `f(x, α, β = 1) = α .+ β*x`
            else ; return true  # TODO Validate. Should be correct
            # else
            #     first_close_bracket_idx =
            #         findfirst(isequal(')'), line)::Union{Nothing, Int}
            #     first_equal_sign_after_bracket_idx =
            #         findfirst(isequal('='), line[first_close_bracket_idx : end])
            #     str_from_close_bracket_to_equal_sign =
            #         line[first_close_bracket_idx : first_close_bracket_idx + first_equal_sign_after_bracket_idx]

            #     if occursin('=', str_from_close_bracket_to_equal_sign)
            #         return true
            #     end
            end

        elseif ! has__open_bracket
            # Test for `==` also covers the case for `===`
            str_to_first_equal_sign_plus_one = line[1 : first_equal_sign_idx + 1]
            is__non_comparison_equal_sign =
                occursin(r"(?<!=)=(?!=)", str_to_first_equal_sign_plus_one)

            if is__non_comparison_equal_sign ; return true
            else ; return false
            end
        end
    end
    # --------------------------------------------

    if match !== nothing
        # state for comment parser, assuming that the `;` isn't in a string or comment
        # so input like ";#" will still thwart this to give the wrong (anti-conservative) answer
        comment = false
        comment_start = false
        comment_close = false
        comment_multi = 0
        for c in line[(match + 1):end]
            if comment_multi > 0
                # handle nested multi-line comments
                if comment_close && c == '#'
                    comment_close = false
                    comment_multi -= 1
                elseif comment_start && c == '='
                    comment_start = false
                    comment_multi += 1
                else
                    comment_start = (c == '#')
                    comment_close = (c == '=')
                end
            elseif comment
                # handle line comments
                if c == '\r' || c == '\n'
                    comment = false
                end
            elseif comment_start
                # see what kind of comment this is
                comment_start = false
                if c == '='
                    comment_multi = 1
                else
                    comment = true
                end
            elseif c == '#'
                # start handling for a comment
                comment_start = true
            else
                # outside of a comment, encountering anything but whitespace
                # means the semi-colon was internal to the expression
                isspace(c) || return false
            end
        end
        return true
    end
    return false
end

print_response_custom = """
function REPL.ends_with_semicolon(line::AbstractString)
    match = findlast(isequal(';'), line)::Union{Nothing,Int}

    # --------------------------------------------
    # CUSTOM Avoid printing for object assignment by testing if the first `=` is inside round brackets `(`
    # --------------------------------------------

    # Test for `DataFrame(a = 1)`
    first_equal_sign_idx = findfirst(isequal('='), line)::Union{Nothing, Int}

    if first_equal_sign_idx == nothing ; return false
    else
        # Test for function, macro & struct
        is__func_macro_struct =
            occursin(r"(\bfunction|@|struct|mutable struct\b)", line)

        str_to_first_equal_sign = line[1 : first_equal_sign_idx]
        has__open_bracket = occursin("(", str_to_first_equal_sign)

        if is__func_macro_struct ; return true
        elseif has__open_bracket ; return false
        elseif ! has__open_bracket
            # Test for `==` also covers the case for `===`
            str_to_first_equal_sign_plus_one = line[1 : first_equal_sign_idx + 1]
            is__non_comparison_equal_sign =
                occursin(r"(?<!=)=(?!=)", str_to_first_equal_sign_plus_one)

            if is__non_comparison_equal_sign ; return true
            else ; return false
            end
        end
    end
    # --------------------------------------------

    if match !== nothing
        # state for comment parser, assuming that the `;` isn't in a string or comment
        # so input like ";#" will still thwart this to give the wrong (anti-conservative) answer
        comment = false
        comment_start = false
        comment_close = false
        comment_multi = 0
        for c in line[(match + 1):end]
            if comment_multi > 0
                # handle nested multi-line comments
                if comment_close && c == '#'
                    comment_close = false
                    comment_multi -= 1
                elseif comment_start && c == '='
                    comment_start = false
                    comment_multi += 1
                else
                    comment_start = (c == '#')
                    comment_close = (c == '=')
                end
            elseif comment
                # handle line comments
                if c == '\r' || c == '\n'
                    comment = false
                end
            elseif comment_start
                # see what kind of comment this is
                comment_start = false
                if c == '='
                    comment_multi = 1
                else
                    comment = true
                end
            elseif c == '#'
                # start handling for a comment
                comment_start = true
            else
                # outside of a comment, encountering anything but whitespace
                # means the semi-colon was internal to the expression
                isspace(c) || return false
            end
        end
        return true
    end
    return false
end
"""

print_response_orig = """
function REPL.ends_with_semicolon(line::AbstractString)
    match = findlast(isequal(';'), line)::Union{Nothing,Int}

    # --------------------------------------------
    # CUSTOM Avoid printing for object assignment by testing if the first `=` is inside round brackets `(`
    # --------------------------------------------

    # Test for `DataFrame(a = 1)`
    first_equal_sign_idx = findfirst(isequal('='), line)::Union{Nothing, Int}

    if first_equal_sign_idx == nothing ; return false
    else
        # Test for function, macro & struct
        is__func_macro_struct =
            occursin(r"(\bfunction|@|struct|mutable struct\b)", line)

        str_to_first_equal_sign = line[1 : first_equal_sign_idx]
        has__open_bracket = occursin("(", str_to_first_equal_sign)

        if is__func_macro_struct ; return true
        elseif has__open_bracket ; return false
        elseif ! has__open_bracket
            # Test for `==` also covers the case for `===`
            str_to_first_equal_sign_plus_one = line[1 : first_equal_sign_idx + 1]
            is__non_comparison_equal_sign =
                occursin(r"(?<!=)=(?!=)", str_to_first_equal_sign_plus_one)

            if is__non_comparison_equal_sign ; return true
            else ; return false
            end
        end
    end
    # --------------------------------------------

    if match !== nothing
        # state for comment parser, assuming that the `;` isn't in a string or comment
        # so input like ";#" will still thwart this to give the wrong (anti-conservative) answer
        comment = false
        comment_start = false
        comment_close = false
        comment_multi = 0
        for c in line[(match + 1):end]
            if comment_multi > 0
                # handle nested multi-line comments
                if comment_close && c == '#'
                    comment_close = false
                    comment_multi -= 1
                elseif comment_start && c == '='
                    comment_start = false
                    comment_multi += 1
                else
                    comment_start = (c == '#')
                    comment_close = (c == '=')
                end
            elseif comment
                # handle line comments
                if c == '\r' || c == '\n'
                    comment = false
                end
            elseif comment_start
                # see what kind of comment this is
                comment_start = false
                if c == '='
                    comment_multi = 1
                else
                    comment = true
                end
            elseif c == '#'
                # start handling for a comment
                comment_start = true
            else
                # outside of a comment, encountering anything but whitespace
                # means the semi-colon was internal to the expression
                isspace(c) || return false
            end
        end
        return true
    end
    return false
end
"""

const mykeys = Dict{Any,Any}(
    # Remove '^D'
    # "^D" => (s,o...)-> nothing,  # TODO Not working

    # Move up/dn, history next/prev
    "\\M-k" => (s,o...)->(LineEdit.edit_move_up(s)   || LineEdit.history_prev(s, LineEdit.mode(s).hist)),
    "\\M-j" => (s,o...)->(LineEdit.edit_move_down(s) || LineEdit.history_next(s, LineEdit.mode(s).hist)),

    # Move char left/right
    "\\M-h" => (s,o...)->(LineEdit.edit_move_left(s)),
    "\\M-l" => (s,o...)->(LineEdit.edit_move_right(s)),

    # Kill line forward/backward
    "\e[1;2D" => (s,o...)->(LineEdit.edit_kill_line_backwards(s)),
    "\e[1;2C" => (s,o...)->(LineEdit.edit_kill_line_forwards(s)),

    # Send Julia history
    # TODO
    # "\\M-r"    => (s,o...)->(LineEdit.enter_search(s, p, true)),
    # "..."    => (s,o...)->(enter_search(s, p, false)),
    # "\\M-r"      => (s,data,c)->(history_set_backward(data, true); history_next_result(s, data)),
    # "..."      => (s,data,c)->(history_set_backward(data, false); history_next_result(s, data)),
    "\\C-r" => (s,o...)->(run(`send_jl_history`);),

    # Insert text
    "\\M-p" => (s,o...)->(LineEdit.edit_insert(s, " |> ")),
    # Activate/restore print_response_custom/print_response_orig TODO Enter
    "\\M-m" => (s,o...)->(LineEdit.edit_insert(s, print_response_custom)),
    # TODO Not working See above
    "\\M-n" => (s,o...)->(LineEdit.edit_insert(s, print_response_orig)),
)

function customize_keys(repl)
    # # repl = Base.active_repl;
    # repl.interface = REPL.setup_interface(repl);
    # main_mode = repl.interface.modes[1]
    # hp = main_mode.hist
    # p = LineEdit.HistoryPrompt(hp)

    repl.interface = REPL.setup_interface(repl; extra_repl_keymap = mykeys)

    repl.interface.modes[2].prompt = "shell❯ ";
    repl.interface.modes[3].prompt = "help❯ ";
end

atreplinit(customize_keys)

# Revise -----------------------------------------

# NOTE Call {Revise} here for watching pkgs listed after this section

# For Julia < v1.5
atreplinit() do repl
    try
        @eval using Revise
        @async Revise.wait_steal_repl_backend()
    catch e
        @warn(e.msg)
    end
end

# TODO OhMyREPL workaround -----------------------

# Reinstall keybindings to work around https://github.com/KristofferC/OhMyREPL.jl/issues/166
@async begin  sleep(1) ; OhMyREPL.Prompt.insert_keybindings()  end

# {ClearStacktrace} ------------------------------

# TODO [Stacktrace gripes #33065](https://github.com/JuliaLang/julia/issues/33065)
using ClearStacktrace

# TODO v0.1.0
ClearStacktrace.CRAYON_HEAD[] = crayon"dark_gray";
ClearStacktrace.CRAYON_HEADSEP[] = crayon"dark_gray";
ClearStacktrace.MODULECRAYONS[] = [
    crayon"#E4CAAF", crayon"green", crayon"light_cyan"];
ClearStacktrace.TYPECOLORS[] = Any[0xCCD9DD, 0xE4CAAF, 0xDCE5EC];

# TODO New version
# ClearStacktrace.MODULECOLORS[] = [:light_blue, :light_yellow, :light_red, :light_green, :light_magenta, :light_cyan, :blue, :yellow, :red, :green, :magenta, :cyan]

# try sum(a) catch end  # TODO Not working

# {DrWatson} -------------------------------------

using DrWatson
# TODO
# if isfile("Project.toml")  @quickactivate  end
using Pkg ; if isfile("Project.toml")  Pkg.activate(".")  end

# TODO
# Pkg.UPDATED_REGISTRY_THIS_SESSION[] = false

# Custom printing --------------------------------

using PrettyPrinting: pprint
Base.show(io::IO, ::MIME"text/plain", mytype::NamedTuple) =
    pprint(io, mytype)

# _ Turing ---------------------------------------

# Adapt from {MCMCChains}
# Turing Chains TODO `using Turing` and sysimg
# TODO Rounding

using Turing, StatsBase
import DataAPI.describe

function StatsBase.summarystats(
    chains::Chains;
    sections = MCMCChains._default_sections(chains),
    append_chains::Bool = true,
    method::MCMCChains.AbstractESSMethod = ESSMethod(),
    maxlag = 250,
    etype = :bm,
    kwargs...
)
    # Store everything.
    cskip = MCMCChains.cskip
    funs = [mean∘cskip,
            std∘cskip,
            x -> quantile(cskip(x), .05),
            x -> quantile(cskip(x), .95)]
            # sem∘cskip,
            # x -> MCMCChains.mcse(cskip(x), etype; kwargs...)]
    func_names = [:mean, :std, Symbol("5%"), Symbol("95%")]  # :naive_se, :mcse]

    # Subset the chain.
    _chains = Chains(chains, MCMCChains._clean_sections(chains, sections))

    # Calculate ESS separately.
    ess_df = ess(_chains; sections = nothing, method = method, maxlag = maxlag)

    # Summarize.
    summary_df = summarize(
        _chains, funs...;
        func_names = func_names,
        append_chains = append_chains,
        additional_df = ess_df,
        name = "Summary Statistics",
        sections = nothing
    )
    return summary_df
end

function describe(
    io::IO,
    chains::Chains;
    etype = :bm,
    kwargs...
)
    dfs = vcat(summarystats(chains; etype = etype, kwargs...))
    return dfs
end

# R plot -----------------------------------------

# For quick plots - without comments (same chunks as above)

# using RCall

# reval("library(data.table)")
# @rlibrary forcats
# @rlibrary ggplot2;
# @rlibrary patchwork
# reval("suppressMessages(library(plotly))");
# reval("suppressMessages(library(theme))");
# @rlibrary workflow;
# @rlibrary ggdist;
# try ggplotly(ggplot(mapping = aes(1, 1)) + geom_line())  catch end
# @rlibrary theme;  # for {theme}-related {plotly} e.g. `theme::ggplotly`

# {PrettyTables} ---------------------------------

# using PrettyTables, DataFrames

# hl_lt_(n::Number) = Highlighter(
#     f = (data,i,j)->begin
#         if applicable(<,data[i,j],n) && data[i,j] < n
#             return true
#         else
#             return false
#         end
#     end,
#     crayon = Crayon(foreground = :light_blue)
# )

# function o(df::DataFrame, n::Real = 10)  # TODO Deal with `Inf` ?
#     n_row, n_col = size(df)
#     println("\033[38;2;228;202;175m$(n_row)×$(n_col)\033[0m")

#     half_n = min(div(n, 2), size(df, 1))
#     head_tail_idx = [collect(1:half_n)..., collect(n_row-half_n+1 : n_row)...]
#     head_tail = df[head_tail_idx, :]

#     eltypes_ = eltype.(eachcol(df))
#     num_cols, str_cols = [findall(x -> x <: type, eltypes_)
#                               for type in [Real, String]]

#     pretty_table(head_tail,
#         header_crayon = Crayon(bold = false),
#         columns_width = 8,  # Dynamic width?
#         row_names = head_tail_idx,
#         tf = borderless,
#         body_hlines = [half_n],
#         formatters      = ft_printf("%.2f", num_cols))  # TODO Decimals
#         # TODO Not working for string cols `String`
#         # highlighters    = (hl_lt_(0)))  # hl_gt(.3)
# end

# # # Test
# # using DataFrames
# # df = DataFrame(A = 1:2:200, B = rand(100), C = -rand(100), D = "haha")

# # @ptconf show_row_number = true tf = borderless formatters = ft_printf("%.3f", [2,3]) highlighters = (hl_lt_(0))
# # @pt df

# # Base.show(df::AbstractDataFrame;
# #           allrows::Bool = !get(stdout, :limit, true),
# #           allcols::Bool = !get(stdout, :limit, true),
# #           splitcols = get(stdout, :limit, true),
# #           rowlabel::Symbol = :Row,
# #           summary::Bool = true) =
# #     show(stdout, pretty_table(df, tf = borderless))
# #         #  allrows=allrows, allcols=allcols, splitcols=splitcols,
# #         #  rowlabel=rowlabel, summary=summary)

# # function _show(io::IO,
# #                df::AbstractDataFrame;
# #                allrows::Bool = !get(io, :limit, false),
# #                allcols::Bool = !get(io, :limit, false),
# #                splitcols = get(io, :limit, false),
# #                rowlabel::Symbol = :Row,
# #                summary::Bool = true,
# #                rowid=nothing)
# #     _check_consistency(df)
# #     nrows = size(df, 1)
# #     if rowid !== nothing
# #         if size(df, 2) == 0
# #             rowid = nothing
# #         elseif nrows != 1
# #             throw(ArgumentError("rowid may be passed only with a single row data frame"))
# #         end
# #     end
# #     dsize = displaysize(io)
# #     availableheight = dsize[1] - 7
# #     nrowssubset = fld(availableheight, 2)
# #     bound = min(nrowssubset - 1, nrows)
# #     if allrows || nrows <= availableheight
# #         rowindices1 = 1:nrows
# #         rowindices2 = 1:0
# #     else
# #         rowindices1 = 1:bound
# #         rowindices2 = max(bound + 1, nrows - nrowssubset + 1):nrows
# #     end
# #     maxwidths = getmaxwidths(df, io, rowindices1, rowindices2, rowlabel, rowid)
# #     width = getprintedwidth(maxwidths)
# #     showrows(io,
# #              df,
# #              rowindices1,
# #              rowindices2,
# #              maxwidths,
# #              splitcols,
# #              allcols,
# #              rowlabel,
# #              summary,
# #              rowid)
# #     return
# # end

# # using DataFrames

# # df = DataFrame(A = 1:2:20, B = rand(10), C = rand(10))

# # pretty_table(df,
# #     tf = borderless,
# #     formatters = ft_printf("%.3f", [2,3]), highlighters = (hl_lt(0.2), hl_gt(0.8)))