text:
export
    RBState,
    randbetween
import TrajectoryOptimization.Dynamics: angular_velocity, linear_velocity, orientation

struct RBState{T}
    r::SVector{3,T}
    q::UnitQuaternion{T,TO.DEFAULT_QUATDIFF}
    v::SVector{3,T}
    ω::SVector{3,T}
end

function RBState(r::AbstractVector, q::Rotation, v::AbstractVector, ω::AbstractVector)
    r_ = @SVector [r[1],r[2],r[3]]
    q_ = UnitQuaternion{TO.DEFAULT_QUATDIFF}(q)
    v_ = @SVector [v[1],v[2],v[3]]
    ω_ = @SVector [ω[1],ω[2],ω[3]]
    RBState(r_, q_, v_, ω_)
end

function RBState(r::AbstractVector, q::AbstractVector, v::AbstractVector, ω::AbstractVector)
    @assert length(q) == 4
    q = UnitQuaternion(q...)
    RBState(r, q, v, ω)
end

function RBState(x::SVector{13})
    r_ = @SVector [x[1],x[2],x[3]]
    q_ = UnitQuaternion(x[4], x[5], x[6], x[7])
    v_ = @SVector [x[8],x[9],x[10]]
    ω_ = @SVector [x[11],x[12],x[13]]
    RBState(r_, q_, v_, ω_)
end

function RBState(model::RigidBody, x::SVector)
    r,q,v,ω = Dynamics.parse_state(model, x)
    RBState(r,UnitQuaternion(q),v,ω)
end

function RBState(model::RigidBody, Z::Traj)
    [RBState(model, state(z)) for z in Z]
end

@inline Base.position(x::RBState) = x.r
@inline orientation(x::RBState) = x.q
@inline angular_velocity(x::RBState) = x.ω
@inline linear_velocity(x::RBState) = x.v

function Dynamics.build_state(model::RigidBody{R}, rbs::RBState) where R
    Dynamics.build_state(model, rbs.r, rbs.q, rbs.v, rbs.ω)
end

function Base.:+(s1::RBState, s2::RBState)
    RBState(s1.r+s2.r, s1.q*s2.q, s1.v+s2.v, s1.ω+s2.ω)
end

function Base.:-(s1::RBState, s2::RBState)
    RBState(s1.r-s2.r, s2.q\s1.q, s1.v-s2.v, s1.ω-s2.ω)
end

function TrajectoryOptimization.:⊖(s1::RBState, s2::RBState, rmap=ExponentialMap)
    dx = s1.r-s2.r
    dq = rmap(s2.q\s1.q)
    dv = s1.v-s2.v
    dw = s1.ω-s2.ω
    @SVector [dx[1], dx[2], dx[3], dq[1], dq[2], dq[3],
              dv[1], dv[2], dv[3], dw[1], dw[2], dw[3]]
end

Base.zero(s1::RBState) = zero(RBState)
function Base.zero(::Type{<:RBState})
    RBState((@SVector zeros(3)), UnitQuaternion(I),
        (@SVector zeros(3)), (@SVector zeros(3)))
end

function Base.rand(::Type{RBState})
    RBState(rand(3), rand(UnitQuaternion), rand(3), rand(3))
end

function randbetween(xmin::RBState, xmax::RBState)
    dx = xmax - xmin
    RBState(
        xmin.r .+ rand(3) .* dx.r,
        # rand(UnitQuaternion),
        expm((@SVector randn(3))*rand()*deg2rad(170)),
        xmin.v .+ rand(3) .* dx.v,
        xmin.ω .+ rand(3) .* dx.ω
    )
end

function LinearAlgebra.norm(s::RBState)
    sqrt(s.r's.r + s.v's.v + s.ω's.ω + LinearAlgebra.norm2(s.q))
end

function TrajectoryOptimization.Traj(model::RigidBody,
        X::Vector{<:RBState}, U::Vector{<:AbstractVector}, dt)
    N = length(X)
    equal = N == length(U)
    map(1:length(X)) do k
        x = Dynamics.build_state(model, X[k])
        if k == N
            if equal
                u = U[k]
            else
                try
                    u = Dynamics.trim_controls(model)
                catch
                    u = zeros(model)[2]
                end
            end
        else
            u = U[k]
        end
        KnotPoint(x,u,dt,dt*(k-1))
    end
end


text:
import MathOptInterface
import HiGHS

const MOI = MathOptInterface
const MOIT = MOI.Test

const CONFIG = MOIT.TestConfig()

@testset "Attributes" begin
    o = HiGHS.Optimizer()
    @test MOI.supports(o, MOI.SolverName())
    @test MOI.get(o, MOI.SolverName()) == "HiGHS"
    @test MOI.get(o, MOI.TimeLimitSec()) > 10000
    MOI.set(o, MOI.TimeLimitSec(), 500)
    @test MOI.get(o, MOI.TimeLimitSec()) == 500.0
    @test MOI.supports(o, MOI.RawSolver())
    @test MOI.get(o, MOI.RawSolver()) == o
end

@testset "MOI variable count and empty" begin
    o = HiGHS.Optimizer()
    x1 = MOI.add_variable(o)
    @test x1.value == 0
    @test MOI.supports_constraint(o, MOI.SingleVariable(x1), MOI.Interval(0, 1))
    (x2, _) = MOI.add_constrained_variable(o, MOI.Interval(0, 1))
    @test x2.value == 1
    @test MOI.get(o, MOI.NumberOfVariables()) == 2
    MOI.empty!(o)
    @test MOI.get(o, MOI.NumberOfVariables()) == 0
end

@testset "Objective function and value in box constraints" begin
    o = HiGHS.Optimizer()
    @testset "Getting objective value" begin
        (x, _) = MOI.add_constrained_variable(o, MOI.Interval(-3.0, 6.0))
        HiGHS.Highs_changeColCost(o, Cint(x.value), 1.0)
        @test MOI.get(o, MOI.ObjectiveSense()) == MOI.FEASIBILITY_SENSE
        MOI.set(o, MOI.ObjectiveSense(), MOI.MIN_SENSE)
        @test MOI.get(o, MOI.ObjectiveSense()) == MOI.MIN_SENSE
        @test MOI.get(o, MOI.ResultCount()) == 0
        MOI.optimize!(o)
        @test MOI.get(o, MOI.ResultCount()) == 1
        @test MOI.get(o, MOI.ObjectiveValue()) ≈ -3
    end
    @testset "Max in box" begin
        MOI.empty!(o)
        @test MOI.get(o, MOI.ResultCount()) == 0
        (x, _) = MOI.add_constrained_variable(o, MOI.Interval(-3.0, 6.0))
        MOI.set(o, MOI.ObjectiveSense(), MOI.MAX_SENSE)
        HiGHS.Highs_changeColCost(o, Cint(x.value), 2.0)
        MOI.optimize!(o)
        @test MOI.get(o, MOI.ObjectiveValue()) ≈ 2 * 6
        obj_func = MOI.get(o, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
        @test MOI.get(o, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test obj_func ≈ MOI.ScalarAffineFunction([
                MOI.ScalarAffineTerm(2.0, x),
            ], 0.0,
        )
    end
    @testset "Objective function obtained from model corresponds" begin
        MOI.empty!(o)
        (x1, _) = MOI.add_constrained_variable(o, MOI.Interval(-3.0, 6.0))
        (x2, _) = MOI.add_constrained_variable(o, MOI.Interval(1.0, 2.0))
        MOI.set(o, MOI.ObjectiveSense(), MOI.MIN_SENSE)
        HiGHS.Highs_changeColCost(o, Cint(x1.value), 2.0)
        HiGHS.Highs_changeColCost(o, Cint(x2.value), -1.0)
        F = MOI.get(o, MOI.ObjectiveFunctionType())
        @test F <: MOI.ScalarAffineFunction{Float64}
        obj_func = MOI.get(o, MOI.ObjectiveFunction{F}())
        @test MOI.supports(o, MOI.ObjectiveFunction{F}())
        @test all(MOI.get(o, MOI.ListOfVariableIndices()) .== [x1, x2])
        @test obj_func ≈ MOI.ScalarAffineFunction([
                MOI.ScalarAffineTerm(2.0, x1),
                MOI.ScalarAffineTerm(-1.0, x2),
            ], 0.0,
        )
        MOI.set(o, MOI.ObjectiveFunction{F}(), obj_func)
        obj_func = MOI.get(o, MOI.ObjectiveFunction{F}())
        @test obj_func ≈ MOI.ScalarAffineFunction([
                MOI.ScalarAffineTerm(2.0, x1),
                MOI.ScalarAffineTerm(-1.0, x2),
            ], 0.0,
        )
        obj_func.terms[1] = MOI.ScalarAffineTerm(3.0, x1)
        MOI.set(o, MOI.ObjectiveFunction{F}(), obj_func)
        obj_func = MOI.get(o, MOI.ObjectiveFunction{F}())
        @test obj_func ≈ MOI.ScalarAffineFunction([
                MOI.ScalarAffineTerm(3.0, x1),
                MOI.ScalarAffineTerm(-1.0, x2),
            ], 0.0,
        )
    end

    @testset "Constrained variable equivalent to add constraint" begin
        MOI.empty!(o)
        x = MOI.add_variable(o)
        _ = MOI.add_constraint(o, MOI.SingleVariable(x), MOI.Interval(-3.0, 6.0))
        HiGHS.Highs_changeColCost(o, Cint(x.value), 1.0)
        MOI.set(o, MOI.ObjectiveSense(), MOI.MIN_SENSE)
        @test MOI.get(o, MOI.ResultCount()) == 0
        MOI.optimize!(o)
        @test MOI.get(o, MOI.ResultCount()) == 1
        @test MOI.get(o, MOI.ObjectiveValue()) ≈ -3
    end

    @testset "Constant in objective function" begin
        MOI.empty!(o)
        x = MOI.add_variable(o)
        _ = MOI.add_constraint(o, MOI.SingleVariable(x), MOI.Interval(-3.0, 6.0))
        MOI.set(o, MOI.ObjectiveSense(), MOI.MIN_SENSE)
        obj_func = MOI.ScalarAffineFunction(
            [MOI.ScalarAffineTerm(1.0, x)], 3.0,
        )
        MOI.set(o, MOI.ObjectiveFunction{typeof(obj_func)}(), obj_func)
        MOI.optimize!(o)
        @test MOI.get(o, MOI.ResultCount()) == 1
        @test MOI.get(o, MOI.ObjectiveValue()) ≈ 0
        obj_func = MOI.get(o, MOI.ObjectiveFunction{typeof(obj_func)}())
        @test MOI.constant(obj_func) ≈ 3
        MOI.set(o, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
        obj_func = MOI.get(o, MOI.ObjectiveFunction{typeof(obj_func)}())
        @test MOI.constant(obj_func) ≈ 0
        @test isempty(obj_func.terms)
    end
end

@testset "Linear constraints" begin
    # max x1 + 2x2
    # st 0 <= x{1,2} <= 5
    # 0 <= x1 + x2 <= 7.5
    o = HiGHS.Optimizer()
    (x1, _) = MOI.add_constrained_variable(o, MOI.Interval(0.0, 5.0))
    (x2, _) = MOI.add_constrained_variable(o, MOI.Interval(0.0, 5.0))
    MOI.set(o, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    func = MOI.ScalarAffineFunction(
        [
            MOI.ScalarAffineTerm(1.0, x1),
            MOI.ScalarAffineTerm(2.0, x2),
        ], 0.0,
    )
    @test MOI.supports_constraint(o, func, MOI.Interval(0, 1))
    MOI.set(o, MOI.ObjectiveFunction{typeof(func)}(), func)
    @test MOI.get(o, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64}}()) == 0
    MOI.add_constraint(o,
        MOI.ScalarAffineFunction(
            [
                MOI.ScalarAffineTerm(1.0, x1),
                MOI.ScalarAffineTerm(1.0, x2),
            ], 0.0,
        ), MOI.Interval(0.0, 7.5),
    )
    @test MOI.get(o, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64}}()) == 1
    MOI.optimize!(o)
    @test MOI.get(o, MOI.ObjectiveValue()) ≈ 12.5
    @test MOI.get(o, MOI.SimplexIterations()) > 0
    @test MOI.get(o, MOI.BarrierIterations()) == 0
end

@testset "Variable names" begin
    o = HiGHS.Optimizer()
    MOIT.variablenames(o, CONFIG)
    MOI.empty!(o)
    y = MOI.add_variable(o)
    MOI.set(o, MOI.VariableName(), y, "y")
    y2 = MOI.get(o, MOI.VariableIndex, "y")
    @test y == y2
    @test MOI.get(o, MOI.VariableIndex, "y0") === nothing
end

@testset "HiGHS custom options" begin
    o = HiGHS.Optimizer()
    @test MOI.supports(o, MOI.RawParameter("solver"))
    @test MOI.get(o, MOI.RawParameter("solver")) == "choose"
    MOI.set(o, MOI.RawParameter("solver"), "simplex")
    @test MOI.get(o, MOI.RawParameter("solver")) == "simplex"
    @test MOI.get(o, MOI.RawParameter("message_level")) == 4
    MOI.set(o, MOI.RawParameter("message_level"), 1)
    @test MOI.get(o, MOI.RawParameter("message_level")) == 1
    @test MOI.get(o, MOI.RawParameter("time_limit")) > 1000
    MOI.set(o, MOI.RawParameter("time_limit"), 1000.0)
    @test MOI.get(o, MOI.RawParameter("time_limit")) == 1000.0
    # unsupported test
    @test_throws ArgumentError MOI.get(o, MOI.RawParameter("wrong_param"))
end

@testset "Model empty" begin
    o = HiGHS.Optimizer()
    @test MOI.is_empty(o)
    MOI.add_variable(o)
    @test !MOI.is_empty(o)
    MOI.empty!(o)
    MOI.set(
        o,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction{Float64}([], 0.0)
    )
    @test MOI.is_empty(o)
    MOI.set(
        o,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction{Float64}([], 3.0)
    )
    @test !MOI.is_empty(o)
    MOI.set(o, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    @test MOI.is_empty(o)
    @test MOI.get(o, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}()) ≈ MOI.ScalarAffineFunction{Float64}([], 0.0)
    x = MOI.add_variable(o)
    MOI.set(
        o,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction{Float64}([MOI.ScalarAffineTerm(1.0, x)], 0.0)
    )
    @test_throws ErrorException MOI.optimize!(o)
end

ext: jl


text:
module ClusterValidityIndices

# Package dependencies
using Statistics: mean
using LinearAlgebra

# Common structures
include("common.jl")

# CVI utility functions
include("utils.jl")

# All of the CVI modules
include("CVI/CVI.jl")

# Export all public names
export

    # Abstract types
    CVI,

    # CVI constructors
    XB,
    DB,
    PS,
    CH,
    cSIL,
    GD53,
    GD43,
    WB,
    rCIP,
    # CONN,

    # CVI functions
    param_inc!,
    param_batch!,
    evaluate!,
    get_icvi!,
    get_cvi!,

    # CVI utilities
    sort_cvi_data,
    relabel_cvi_data

    # Not exported
    # get_cvi_data
    # get_bernoulli_subset
    # showtypetree

end


text:
using RCall

# R"require(OpenMx)"
R"ToPSD <- function(distance, frobenius = FALSE){
  DistanceMatrix <- as.matrix(distance)
  N <- dim(DistanceMatrix)[1]
  D2 <- DistanceMatrix^2
  I <- diag(N)
  one <- as.matrix(rep(1, N))
  one.t <- t(one)
  K <- -0.5*(I - one %*% one.t/N) %*% D2 %*% (I - one%*%one.t/N)
	if (frobenius){
  	f = norm(K, type = 'F')
  	K = K/f
  }
  evecK = as.matrix(eigen(K)$vector)
  evalK = as.vector(abs(eigen(K)$values))
  evalK[evalK < 1e-4] = 1e-4
  diagK = diag(evalK)
  psdK = evecK %*% diagK %*% t(evecK)
  return(psdK)
}"

function KernelMatrix(arg...; tableFile::String = "", treeFile::String = "",
  out::String = "", format::String = "nwk", alpha::Array{Float64,1} = [1.0],
  Dtype::String = "d_UW", Kernel::Bool = true, frobenius::Bool = true)

  tree = readNewick(filepath = treeFile, format = format)
  if Dtype == "d_VAW"
    VAW = true
  else
    VAW = false
  end

  Distance = GUniFrac(otuTableFile = tableFile, phy = tree, alpha = alpha, VAW = VAW)

  if Dtype == "d_UW"
    D = Distance[Dtype]
  elseif Dtype == "d_alpha"
    D = Distance[[string("d_", alpha)]]
  elseif Dtype == "d_VAW"
    D = Distance[Dtype]
  end

  if Kernel
    @rput D
    @rput out

    if frobenius
      R"K <- ToPSD(distance = D, frobenius = TRUE)"
    else
      R"K <- ToPSD(distance = D)"
    end
    @rget K
    writecsv(string("K_",out,".csv"), K)
  #  R"write.table(K, paste0('Kernel_', out '.csv'), row.names = F, col.names = F, sep=',')"
    return K
  else
    writecsv(string("D_",out,".csv"), D)
    return D
  end
end


text:
# ---
# title: 899. Orderly Queue
# id: problem899
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math, String
# link: <https://leetcode.com/problems/orderly-queue/description/>
# hidden: true
# ---
#
# A string `S` of lowercase letters is given.  Then, we may make any number of
# _moves_.
#
# In each move, we choose one of the first `K` letters (starting from the left),
# remove it, and place it at the end of the string.
#
# Return the lexicographically smallest string we could have after any number of
# moves.
#
#
#
# **Example 1:**
#
#
#
#     Input: S = "cba", K = 1
#     Output: "acb"
#     Explanation:
#     In the first move, we move the 1st character ("c") to the end, obtaining the string "bac".
#     In the second move, we move the 1st character ("b") to the end, obtaining the final result "acb".
#
#
# **Example 2:**
#
#
#
#     Input: S = "baaca", K = 3
#     Output: "aaabc"
#     Explanation:
#     In the first move, we move the 1st character ("b") to the end, obtaining the string "aacab".
#     In the second move, we move the 3rd character ("c") to the end, obtaining the final result "aaabc".
#
#
#
#
# **Note:**
#
#   1. `1 <= K <= S.length <= 1000`
#   2. `S` consists of lowercase letters only.
#
#
## @lc code=start
using LeetCode

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


text:
using Photometry
using BenchmarkTools
using CSV
using DataFrames
using ProgressLogging

data = randn(512, 512) .+ 10

rows = []

@progress for N in [1, 10, 50, 100, 200, 400, 500, 1000, 2000]
    apers = CircularAperture.(256, 256, 10 .* rand(N))
    time = @belapsed aperture_photometry($apers, $data)
    push!(rows, (N=N, time=time))
end

results = DataFrame(rows)
outfile = "julia_circle_apertures_nt=$(Threads.nthreads()).csv"
@info "saving to $outfile"
CSV.write(joinpath(@__DIR__, outfile), results)

nothing


text:
using Gadfly, DataFrames

set_default_plot_size(6inch, 3inch)

xs = 0:0.1:20

df = DataFrame(
    x=xs,
    y=cos.(xs),
    ymin=cos.(xs) .- 0.5,
    ymax=cos.(xs) .+ 0.5,
)

plot(df, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon)


text:
#
# IndexedVector tests
#

srand(1)
a = DataArray(rand(1:5,20))
a[1:2] = NA
ia = IndexedVector(a)
ia2 = IndexedVector(rand(1:4,20))

ia .== 4
v = [1:20]
@assert v[ia .== 4] == v[a .== 4]
@assert sort(v[(ia .== 4) | (ia .== 5)]) == v[(a .== 4) | (a .== 5)]
@assert sort(v[(ia .>= 4) & (ia .== 5)]) == v[(a .>= 4) & (a .== 5)]
@assert sort(v[!(ia .== 4)]) == v[!(a .== 4)]

df = DataFrame(quote
         x1 = IndexedVector(vcat(fill([1:5],4)...))
         x2 = IndexedVector(vcat(fill(letters[1:10],2)...))
     end)

df[:(x2 .== "a"), :]
df[:( (x2 .== "a") | (x1 .== 2) ), :]
df[:( ("b" .<= x2 .<= "c") | (x1 .== 5) ), :]
df[:( (x1 .== 1) & (x2 .== "a") ), :]

df[:( in(x2, ["c","e"]) ), :]


text:
# scalar-valued kernel
Base.maximum(kernel::CalibrationErrors.ExponentialKernel) = exp(zero(kernel.γ))
Base.maximum(kernel::CalibrationErrors.SquaredExponentialKernel) = exp(zero(kernel.γ))

# matrix-valued kernels
Base.maximum(kernel::CalibrationErrors.UniformScalingKernel) =
    kernel.λ * maximum(kernel.kernel)
Base.maximum(kernel::CalibrationErrors.DiagonalKernel) =
    maximum(kernel.diag) * maximum(kernel.kernel)

# by default assume p = q
Base.maximum(kce::SKCE) = 2 * maximum(CalibrationErrors.kernel(kce))


text:
@testset "218.the-skyline-problem.jl" begin
    @test get_skyline([[2, 9, 10], [3, 7, 15], [5, 12, 12], [15, 20, 10], [19, 24, 8]]) ==
          [(2, 10), (3, 15), (7, 12), (12, 0), (15, 10), (20, 8), (24, 0)]
    @test get_skyline([[0, 2, 3], [2, 5, 3]]) == [(0, 3), (5, 0)]
end


text:
# TODO use 0.7 logging capabilities!
using Printf
# info, debug, gc, signals, performance
const log_level = RefValue((false, false, false, false, false, false))

function enable_ith(i, value::Bool)
    log_level[] = ntuple(j-> j == i ? value : log_level[][j], Val(5))
    value
end
log_info(value::Bool = true) = enable_ith(1, value)

const logging_io = RefValue(stdout)

macro log_gc(args...)
    quote
        if log_level[][3]
            save_print(logging_io[], $(esc.(args)...))
        end
    end
end
macro log_signals(args...)
    quote
        if log_level[][4]
            save_print(logging_io[], $(esc.(args)...))
        end
    end
end

function print_stats(io::IO, elapsedtime, bytes, gctime, allocs)
    @printf(io, "%10.6f seconds", elapsedtime/1e9)
    if bytes != 0 || allocs != 0
        bytes, mb = Base.prettyprint_getunits(bytes, length(Base._mem_units), Int64(1024))
        allocs, ma = Base.prettyprint_getunits(allocs, length(Base._cnt_units), Int64(1000))
        if ma == 1
            @printf(io, " (%d%s allocation%s: ", allocs, Base._cnt_units[ma], allocs==1 ? "" : "s")
        else
            @printf(io, " (%.2f%s allocations: ", allocs, Base._cnt_units[ma])
        end
        if mb == 1
            @printf(io, "%d %s%s", bytes, Base._mem_units[mb], bytes==1 ? "" : "s")
        else
            @printf(io, "%.3f %s", bytes, Base._mem_units[mb])
        end
        if gctime > 0
            @printf(io, ", %.2f%% gc time", 100*gctime/elapsedtime)
        end
        print(io, ")")
    elseif gctime > 0
        @printf(io, ", %.2f%% gc time", 100*gctime/elapsedtime)
    end
    println(io)
end

macro log_performance(name, expr)
    quote
        if log_level[][5]
            gc_stats = Base.gc_num()
            t_stats = Base.time_ns()
            val = $(esc(expr))
            time_diff = Base.time_ns() - t_stats
            gc_diff = Base.GC_Diff(Base.gc_num(), gc_stats)
            print(logging_io[], $name, ": ")
            print_stats(
                logging_io[], time_diff, gc_diff.allocd,
                gc_diff.total_time, Base.gc_alloc_count(gc_diff)
            )
        else
            val = $(esc(expr))
        end
        return val
    end
end

const io_lock = ReentrantLock()
save_print(args...) = save_print(io, args...)
function save_print(io::IO, args...)
    @async begin
        try
            lock(io)
            lock(io_lock)
            print(io, string(args..., "\n"))
        finally
            unlock(io_lock)
            unlock(io)
        end
    end
end


text:
include("../src/CleverAlgorithms.jl")

using CA

val = CA.Stochastic.randdom_search((x) -> x^2, (-5,5))
println(val)
val = CA.Stochastic.randdom_search((x) -> x^2, (-5,5); op = >)
println(val)


text:
using BenchmarkTools, MicroscopyTools

function bin_3(arr::AbstractArray{T, N}, factors=ntuple(x -> 2, Val(N))) where {T, N}
    # old size
    sz = size(arr)
    # new size
    sn = sz .÷ factors
    s2 = ntuple(i -> isodd(i) ? factors[(i+1)/2] : sn[i/2], Val(2 * N))

    dims = ntuple(i -> 2 * i - 1, Val(N))
    reshape(sum(reshape(arr,s2),dims=dims),sn)
end

function f(x)
    @info "Random data with size $(size(x))"
    print("Tullio: ")
    @btime bin2($x)
    print("CartesianIndices: ")
    @btime bin($x, (2,2))
    print("RH: ")
    @btime bin_3($x, (2,2))
end

function main()
    @info "Binning by (2,2)"
    x = randn(Float32, (256, 512))
    f(x)

    x = randn(Float32, (8192, 7113))
    f(x)

    x = randn((97, 32))
    f(x)

    x = randn((3097, 113, 32))
    f(x)

    x = randn((23, 31, 8))
    f(x)

    x = randn(Float32, (3097, 4096, 32))
    f(x)

    return
end

main()

max_stars_repo_name: JuliaMicroscopy/MicroscopyTools

max_issues_repo_name: JuliaMicroscopy/MicroscopyTools


text:
module PolygonMinkowskiSum
export minkowski_sum

using GeometryTypes: Point
using LinearAlgebra: normalize, dot

"""
Polygon with no holes
"""
struct SimplePolygon{VP}
    points::VP
end

function diff(p::SimplePolygon)
    [p.points[mod1(n + 1, end)] - p.points[n] for n=1:length(p.points)]
end

function crosses(p::SimplePolygon)
    d = diff(p)
    [cross(d[i], d[mod1(i+1, end)]) for i=1:length(d)]
end

function find_ccw_violations(p::SimplePolygon)
    # TODO have option for strict or not (allow co-linear points?)
    cs = crosses(p)
    bad_i = findall(cs .< 0) # TODO generic zero
    return bad_i
end

is_ccw(p1, p2, p3) = cross(p1, p2) ≥ 0 && cross(p2, p3) ≥ 0

# @recipe macro requires RecipesBase to be defined in scope of macro use
using RecipesBase: RecipesBase, @recipe

@recipe function f(p::SimplePolygon)
    cycle = p.points[[1:end...,1]]
    seriestype --> :shape
    (x, y) = ([p[1] for p in cycle], [p[2] for p in cycle])
    x, y
end

# TODO doesn't work for `Array`. What is a generic constructor?
rot_quarter(v::P) where {P} = P(-v[2], v[1]) # TODO constrain P

# asin(angle from v1 to v2) * |v1| * |v2|
cross(v1, v2) = v1[1] * v2[2] - v1[2] * v2[1]

"""
find the vector that points in the same direction (positive dot product),
is ccw from v (positive cross product), that is next (smallest cross product)
"""
function arg_next_direction(vs, v)
    I = eachindex(vs)
    i_best = nothing
    c_best = typemax(Float64) # TODO generic
    for i in I
        v_i = normalize(vs[i])
        if dot(v_i, v) > 0
            c = cross(v_i, v)
            if c > 0 && c < c_best
                c_best = c
                i_best = i
            end
        end
    end
    return i_best
end

"""
Assume all point sequences are arranged counter-clockwise
Assume both shapes are convex

Effectively https://en.wikipedia.org/wiki/Minkowski_addition#Two_convex_polygons_in_the_plane
"""
function minkowski_sum(this::SimplePolygon, alongthis::SimplePolygon, plotdebug=nothing)
    IJ = find_minkowski_sum(this, alongthis)
    r = SimplePolygon([this.points[i] + alongthis.points[j] for (i,j) in IJ])
    if plotdebug !== nothing
        ss = [string((k,i,j)) for (k, (i, j)) in enumerate(IJ)]
        plotdebug.ns.scatter!(plotdebug.plots["result"], r.points, series_annotation=ss)
        plotdebug.ns.scatter!(plotdebug.plots["poly1"], this.points, series_annotation=map(string,1:length(this.points)))
        plotdebug.ns.scatter!(plotdebug.plots["poly2"], alongthis.points, series_annotation=map(string,1:length(alongthis.points)))
    end
    return r
end

"""
only works if there are no co-linear edges
"""
function sort_find_sort_polar(p1::SimplePolygon, p2::SimplePolygon)
    polar_i(p) = [(angle=atan(y,x), index=i) for (i,(x,y)) in enumerate(diff(p))]

    (s1, s2) = ([(vertex=pi, index_polygon=1) for pi in polar_i(p)] for p in (p1, p2))
    # TODO these sorts need to be stable with repsect to cyclic order, and they
    s1 = sort(s1)
    s2 = sort(s2)
    s1, s2, sort(vcat(s1, s2))
end

function find_minkowski_sum2(p1::SimplePolygon, p2::SimplePolygon)
    ps = (p1, p2)
    s1, s2, sboth = sort_find_sort_polar(ps...)

    if first(sboth) == first(s2)
        return [(i2, i1) for (i1, i2) in find_minkowski_sum(p2, p1)]
    end
    @assert first(sboth) == first(s1)

    r = Tuple{Int64, Int64}[]

    is = [s1[1].vertex.index, s2[end].vertex.index]

    j = 2
    while true
        push!(r, tuple(is...))
        if !(j <= length(sboth)); break; end
        s = sboth[j]
        is[s.index_polygon] = s.vertex.index
        j += 1
    end
    is = [s1[end].vertex.index, s2[1].vertex.index]
    push!(r, tuple(is...))
    return r
end

function find_minkowski_sum(this::SimplePolygon, alongthis::SimplePolygon)
    result = Tuple{Int64,Int64}[]
    n = 1

    d = alongthis.points[n + 1] - alongthis.points[n]
    # TODO don't create an array for these differences.
    this_diff = diff(this)

    # find two edges in `this` that straddle `d`
    is_ccws = [is_ccw(this_diff[mod1(i, end)], d, this_diff[mod1(i+1, end)]) for i = 1:length(this_diff)]
    # index of `this` vertex that is incident to both edges
    m = mod1(findfirst(is_ccws) + 1, length(is_ccws))
    push!(result, (m, n))

    while length(result) < length(this.points) + length(alongthis.points)
        # TODO avoid recomputing
        # TODO check for repeated points (causes div-by-zero in `normalize`)
        d_path = normalize(alongthis.points[mod1(n+1, end)] - alongthis.points[n])
        d_poly = normalize(this.points[mod1(m+1, end)] - this.points[mod1(m, end)])
        # TODO these assertions are being violated, understand why they are wrong or why the algorithm is wrong
        #@assert cross(d, d_path) >= 0
        #@show cross(d, d_poly)
        #@assert cross(d, d_poly) >= 0
        if cross(d, d_path) < cross(d, d_poly)
            n += 1
            d = d_path
        else
            m += 1
            d = d_poly
        end
        push!(result,
            (mod1(m, length(this.points)), mod1(n, length(alongthis.points)))
        )
    end
    return result
end
end # module


text:
using LinearAlgebra
function sfvf(x, y)
    sign(x) * max(abs(x) - y, 0)
end
function fit(x, t; alpha = 0.1, tol = 0.0001, mi = 1000000)
    function update(n, d, x, t, w, alpha)
        l = length(w)
        w[1] = sum(t - x * w[2:l]) / n
        wvec = ones(n) * w[1]
        for k = 1:d
            ww = w[2:l]
            ww[k] = 0
            q = dot(vec(t) - wvec - x * ww, x[:, k])
            r = dot(x[:, k], x[:, k])
            w[k+1] = sfvf(q / r, alpha)
        end
    end
    if size(x)[1] != size(t)[1]
        x = x'
        t = t'
    end
    n, d = size(x)
    w = zeros(d + 1)
    e = 0.0
    for i = 1:mi
        eb = e
        update(n, d, x, t, w, alpha)
        e = sum(broadcast(abs, w)) / size(w)[1]
        if abs(e - eb) <= tol
            return w
        end
    end
end
function predict(x, w)
    x2 = ones(1, size(x)[2])
    x = vcat(x2, x)
    return w' * x
end


text:
struct DeivVecTag end

# J(f(x))*v
function auto_jacvec!(dy, f, x, v,
                      cache1 = ForwardDiff.Dual{DeivVecTag}.(x, v),
                      cache2 = ForwardDiff.Dual{DeivVecTag}.(x, v))
    cache1 .= Dual{DeivVecTag}.(x, v)
    f(cache2,cache1)
    dy .= partials.(cache2, 1)
end
function auto_jacvec(f, x, v)
    partials.(f(Dual{DeivVecTag}.(x, v)), 1)
end

function num_jacvec!(dy,f,x,v,cache1 = similar(v),
                     cache2 = similar(v);
                     compute_f0 = true)
    compute_f0 && (f(cache1,x))
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    @. x += ϵ*v
    f(cache2,x)
    @. x -= ϵ*v
    @. dy = (cache2 - cache1)/ϵ
end

function num_jacvec(f,x,v,f0=nothing)
    f0 === nothing ? _f0 = f(x) : _f0 = f0
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(minimum(x)))
    (f(x.+ϵ.*v) .- _f0)./ϵ
end

function num_hesvec!(dy,f,x,v,
                     cache1 = similar(v),
                     cache2 = similar(v),
                     cache3 = similar(v))
    cache = FiniteDiff.GradientCache(v[1],cache1,Val{:central})
    g = let f=f,cache=cache
        (dx,x) -> FiniteDiff.finite_difference_gradient!(dx,f,x,cache)
    end
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    @. x += ϵ*v
    g(cache2,x)
    @. x -= 2ϵ*v
    g(cache3,x)
    @. dy = (cache2 - cache3)/(2ϵ)
end

function num_hesvec(f,x,v)
    g = (x) -> FiniteDiff.finite_difference_gradient(f,x)
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    x += ϵ*v
    gxp = g(x)
    x -= 2ϵ*v
    gxm = g(x)
    (gxp - gxm)/(2ϵ)
end

function numauto_hesvec!(dy,f,x,v,
                     cache = ForwardDiff.GradientConfig(f,v),
                     cache1 = similar(v),
                     cache2 = similar(v))
    g = let f=f,x=x,cache=cache
        g = (dx,x) -> ForwardDiff.gradient!(dx,f,x,cache)
    end
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    @. x += ϵ*v
    g(cache1,x)
    @. x -= 2ϵ*v
    g(cache2,x)
    @. dy = (cache1 - cache2)/(2ϵ)
end

function numauto_hesvec(f,x,v)
    g = (x) -> ForwardDiff.gradient(f,x)
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    x += ϵ*v
    gxp = g(x)
    x -= 2ϵ*v
    gxm = g(x)
    (gxp - gxm)/(2ϵ)
end

function autonum_hesvec!(dy,f,x,v,
                     cache1 = ForwardDiff.Dual{DeivVecTag}.(x, v),
                     cache2 = ForwardDiff.Dual{DeivVecTag}.(x, v))
    cache = FiniteDiff.GradientCache(v[1],cache1,Val{:central})
    g = (dx,x) -> FiniteDiff.finite_difference_gradient!(dx,f,x,cache)
    cache1 .= Dual{DeivVecTag}.(x, v)
    g(cache2,cache1)
    dy .= partials.(cache2, 1)
end

function autonum_hesvec(f,x,v)
    g = (x) -> FiniteDiff.finite_difference_gradient(f,x)
    partials.(g(Dual{DeivVecTag}.(x, v)), 1)
end

function num_hesvecgrad!(dy,g,x,v,
                     cache2 = similar(v),
                     cache3 = similar(v))
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    @. x += ϵ*v
    g(cache2,x)
    @. x -= 2ϵ*v
    g(cache3,x)
    @. dy = (cache2 - cache3)/(2ϵ)
end

function num_hesvecgrad(g,x,v)
    T = eltype(x)
    # Should it be min? max? mean?
    ϵ = sqrt(eps(real(T))) * max(one(real(T)), abs(norm(x)))
    x += ϵ*v
    gxp = g(x)
    x -= 2ϵ*v
    gxm = g(x)
    (gxp - gxm)/(2ϵ)
end

function auto_hesvecgrad!(dy,g,x,v,
                     cache2 = ForwardDiff.Dual{DeivVecTag}.(x, v),
                     cache3 = ForwardDiff.Dual{DeivVecTag}.(x, v))
    cache2 .= Dual{DeivVecTag}.(x, v)
    g(cache3,cache2)
    dy .= partials.(cache3, 1)
end

function auto_hesvecgrad(g,x,v)
    partials.(g(Dual{DeivVecTag}.(x, v)), 1)
end

### Operator Forms

struct JacVec{F,T1,T2,xType}
    f::F
    cache1::T1
    cache2::T2
    x::xType
    autodiff::Bool
end

function JacVec(f,x::AbstractArray;autodiff=true)
    if autodiff
        cache1 = ForwardDiff.Dual{DeivVecTag}.(x, x)
        cache2 = ForwardDiff.Dual{DeivVecTag}.(x, x)
    else
        cache1 = similar(x)
        cache2 = similar(x)
    end
    JacVec(f,cache1,cache2,x,autodiff)
end

Base.size(L::JacVec) = (length(L.cache1),length(L.cache1))
Base.size(L::JacVec,i::Int) = length(L.cache1)
Base.:*(L::JacVec,v::AbstractVector) = L.autodiff ? auto_jacvec(_x->L.f(_x),L.x,v) : num_jacvec(_x->L.f(_x),L.x,v)

function LinearAlgebra.mul!(dy::AbstractVector,L::JacVec,v::AbstractVector)
    if L.autodiff
        auto_jacvec!(dy,(_y,_x)->L.f(_y,_x),L.x,v,L.cache1,L.cache2)
    else
        num_jacvec!(dy,(_y,_x)->L.f(_y,_x),L.x,v,L.cache1,L.cache2)
    end
end

struct HesVec{F,T1,T2,xType}
    f::F
    cache1::T1
    cache2::T2
    cache3::T2
    x::xType
    autodiff::Bool
end

function HesVec(f,x::AbstractArray;autodiff=true)
    if autodiff
        cache1 = ForwardDiff.GradientConfig(f,x)
        cache2 = similar(x)
        cache3 = similar(x)
    else
        cache1 = similar(x)
        cache2 = similar(x)
        cache3 = similar(x)
    end
    HesVec(f,cache1,cache2,cache3,x,autodiff)
end

Base.size(L::HesVec) = (length(L.cache2),length(L.cache2))
Base.size(L::HesVec,i::Int) = length(L.cache2)
Base.:*(L::HesVec,v::AbstractVector) = L.autodiff ? numauto_hesvec(L.f,L.x,v) : num_hesvec(L.f,L.x,v)

function LinearAlgebra.mul!(dy::AbstractVector,L::HesVec,v::AbstractVector)
    if L.autodiff
        numauto_hesvec!(dy,L.f,L.x,v,L.cache1,L.cache2,L.cache3)
    else
        num_hesvec!(dy,L.f,L.x,v,L.cache1,L.cache2,L.cache3)
    end
end

struct HesVecGrad{G,T1,T2,uType}
    g::G
    cache1::T1
    cache2::T2
    x::uType
    autodiff::Bool
end

function HesVecGrad(g,x::AbstractArray;autodiff=false)
    if autodiff
        cache1 = ForwardDiff.Dual{DeivVecTag}.(x, x)
        cache2 = ForwardDiff.Dual{DeivVecTag}.(x, x)
    else
        cache1 = similar(x)
        cache2 = similar(x)
    end
    HesVecGrad(g,cache1,cache2,x,autodiff)
end

Base.size(L::HesVecGrad) = (length(L.cache2),length(L.cache2))
Base.size(L::HesVecGrad,i::Int) = length(L.cache2)
Base.:*(L::HesVecGrad,v::AbstractVector) = L.autodiff ? auto_hesvecgrad(L.g,L.x,v) : num_hesvecgrad(L.g,L.x,v)

function LinearAlgebra.mul!(dy::AbstractVector,L::HesVecGrad,v::AbstractVector)
    if L.autodiff
        auto_hesvecgrad!(dy,L.g,L.x,v,L.cache1,L.cache2)
    else
        num_hesvecgrad!(dy,L.g,L.x,v,L.cache1,L.cache2)
    end
end


text:
function skiptofield!(buf, pos, len, options, row, header)
    while row < header
        while pos <= len
            _, code, _, _, tlen = Parsers.xparse(String, buf, pos, len, options)
            pos += tlen
            Parsers.delimited(code) && break
        end
        row += 1
    end
    return row, pos
end

function countfields(buf, pos, len, options)
    rows = 0
    while pos <= len
        _, code, _, _, tlen = Parsers.xparse(String, buf, pos, len, options)
        pos += tlen
        rows += 1
        Parsers.delimited(code) && continue
        (Parsers.newline(code) || pos > len) && break
    end
    return rows, pos
end

function columnname(buf, vpos, vlen, code, options, i)
    if Parsers.sentinel(code)
        return "Column$i"
    elseif Parsers.escapedstring(code)
        return unescape(PointerString(pointer(buf, vpos), vlen), options.e)
    else
        return unsafe_string(pointer(buf, vpos), vlen)
    end
end

function datalayout_transpose(header, buf, pos, len, options, datarow, normalizenames)
    if isa(header, Integer) && header > 0
        # skip to header column to read column names
        row, pos = skiptofield!(buf, pos, len, options, 1, header)
        # io now at start of 1st header cell
        _, code, vpos, vlen, tlen = Parsers.xparse(String, buf, pos, len, options)
        columnnames = [columnname(buf, vpos, vlen, code, options, 1)]
        pos += tlen
        row, pos = skiptofield!(buf, pos, len, options, header+1, datarow)
        columnpositions = Int64[pos]
        datapos = pos
        rows, pos = countfields(buf, pos, len, options)

        # we're now done w/ column 1, if EOF we're done, otherwise, parse column 2's column name
        cols = 1
        while pos <= len
            # skip to header column to read column names
            row, pos = skiptofield!(buf, pos, len, options, 1, header)
            cols += 1
            _, code, vpos, vlen, tlen = Parsers.xparse(String, buf, pos, len, options)
            push!(columnnames, columnname(buf, vpos, vlen, code, options, cols))
            pos += tlen
            row, pos = skiptofield!(buf, pos, len, options, header+1, datarow)
            push!(columnpositions, pos)
            pos = readline!(buf, pos, len, options)
        end
    elseif isa(header, AbstractRange)
        # column names span several columns
        columnpositions = Int64[]
        columnnames = String[]
        throw(ArgumentError("not implemented for transposed csv files"))
    elseif pos > len
        # emtpy file, use column names if provided
        datapos = pos
        columnpositions = Int64[]
        columnnames = header isa Vector && !isempty(header) ? String[string(x) for x in header] : []
        rows = 0
    else
        # column names provided explicitly or should be generated, they don't exist in data
        # skip to datarow
        row, pos = skiptofield!(buf, pos, len, options, 1, datarow)
        # io now at start of 1st data cell
        columnnames = [isa(header, Integer) || isempty(header) ? "Column1" : string(header[1])]
        columnpositions = Int64[pos]
        datapos = pos
        rows, pos = countfields(buf, pos, len, options)
        # we're now done w/ column 1, if EOF we're done, otherwise, parse column 2's column name
        cols = 1
        while pos <= len
            # skip to datarow column
            row, pos = skiptofield!(buf, pos, len, options, 1, datarow)
            cols += 1
            push!(columnnames, isa(header, Integer) || isempty(header) ? "Column$cols" : string(header[cols]))
            push!(columnpositions, pos)
            pos = readline!(buf, pos, len, options)
        end
    end
    return rows, makeunique(map(x->normalizenames ? normalizename(x) : Symbol(x), columnnames))::Vector{Symbol}, columnpositions
end

function datalayout(header::Integer, buf, pos, len, options, datarow, normalizenames, cmt, ignoreemptylines)
    # default header = 1
    if header <= 0
        # no header row in dataset; skip to data to figure out # of columns
        pos = skipto!(buf, pos, len, options, 1, datarow)
        datapos = Int64(pos)
        fields, pos = readsplitline(buf, pos, len, options, cmt, ignoreemptylines)
        columnnames = [Symbol(:Column, i) for i = eachindex(fields)]
    else
        pos = skipto!(buf, pos, len, options, 1, header)
        fields, pos = readsplitline(buf, pos, len, options, cmt, ignoreemptylines)
        columnnames = makeunique([normalizenames ? normalizename(x) : Symbol(x) for (i, x) in enumerate(fields)])
        if datarow != header+1
            pos = skipto!(buf, pos, len, options, header+1, datarow)
        end
        datapos = Int64(pos)
    end
    return columnnames, datapos::Int64
end

function datalayout(header::AbstractVector{<:Integer}, buf, pos, len, options, datarow, normalizenames, cmt, ignoreemptylines)
    pos = skipto!(buf, pos, len, options, 1, header[1])
    columnnames, pos = readsplitline(buf, pos, len, options, cmt, ignoreemptylines)
    for row = 2:length(header)
        pos = skipto!(buf, pos, len, options, 1, header[row] - header[row-1])
        fields, pos = readsplitline(buf, pos, len, options, cmt, ignoreemptylines)
        for (i, x) in enumerate(fields)
            columnnames[i] *= "_" * x
        end
    end
    if datarow != last(header)+1
        pos = skipto!(buf, pos, len, options, last(header)+1, datarow)
    end
    datapos = pos
    return makeunique([normalizenames ? normalizename(nm) : Symbol(nm) for nm in columnnames]), datapos
end

function datalayout(header::Union{Vector{Symbol}, Vector{String}}, buf, pos, len, options, datarow, normalizenames, cmt, ignoreemptylines)
    pos = skipto!(buf, pos, len, options, 1, datarow)
    datapos = pos
    if pos > len
        columnnames = makeunique([normalizenames ? normalizename(nm) : Symbol(nm) for nm in header])
    else
        fields, pos = readsplitline(buf, pos, len, options, cmt, ignoreemptylines)
        if isempty(header)
            columnnames = [Symbol("Column$i") for i in eachindex(fields)]
        else
            length(header) == length(fields) || throw(ArgumentError("The length of provided header ($(length(header))) doesn't match the number of columns at row $datarow ($(length(fields)))"))
            columnnames = makeunique([normalizenames ? normalizename(nm) : Symbol(nm) for nm in header])
        end
    end
    return columnnames, datapos
end

# readline! is used for implementation of skipto!
function readline!(buf, pos, len, options)
    while pos <= len
        _, code, _, _, tlen = Parsers.xparse(String, buf, pos, len, options)
        pos += tlen
        (Parsers.newline(code) || pos > len) && break
    end
    return pos
end

function skipto!(buf, pos, len, options, cur, dest)
    cur >= dest && return pos
    for _ = 1:(dest-cur)
        pos = readline!(buf, pos, len, options)
    end
    return pos
end

function readsplitline(buf, pos, len, options::Parsers.Options{ignorerepeated}, cmt, ignoreemptylines) where {ignorerepeated}
    vals = String[]
    pos > len && return vals, pos
    col = 1
    while true
        pos = consumecommentedline!(buf, pos, len, cmt, ignoreemptylines)
        if ignorerepeated
            pos = Parsers.checkdelim!(buf, pos, len, options)
        end
        _, code, vpos, vlen, tlen = Parsers.xparse(String, buf, pos, len, options)
        push!(vals, columnname(buf, vpos, vlen, code, options, col))
        pos += tlen
        col += 1
        Parsers.delimited(code) && continue
        (Parsers.newline(code) || pos > len) && break
    end
    return vals, pos
end

@inline function skipemptyrow!(buf, pos, len)
    @inbounds b = buf[pos]
    if b == UInt8('\n')
        return pos + 1
    elseif b == UInt8('\r')
        if pos + 1 < len && buf[pos + 1] == UInt8('\n')
            return pos + 2
        else
            return pos + 1
        end
    end
    return pos
end

function consumecommentedline!(buf, pos, len, ::Nothing, ::Val{true})
    while pos <= len
        newpos = skipemptyrow!(buf, pos, len)
        newpos == pos && break
        pos = newpos
    end
    return pos
end

function consumecommentedline!(buf, pos, len, (cmtptr, cmtlen), ::Val{true})
    ptr = pointer(buf, pos)
    while (pos + cmtlen - 1) <= len
        match = Parsers.memcmp(ptr, cmtptr, cmtlen)
        if match
            pos += cmtlen
            pos > len && break
            @inbounds b = buf[pos]
            while b != UInt8('\n') && b != UInt8('\r')
                pos += 1
                pos > len && break
                @inbounds b = buf[pos]
            end
            pos += 1
        else
            newpos = skipemptyrow!(buf, pos, len)
            if newpos > pos
                pos = newpos
            else
                break
            end
        end
        ptr = pointer(buf, pos)
    end
    return pos
end

consumecommentedline!(buf, pos, len, ::Nothing, ::Val{false}) = pos
function consumecommentedline!(buf, pos, len, (cmtptr, cmtlen), ::Val{false})
    ptr = pointer(buf, pos)
    while (pos + cmtlen - 1) <= len
        match = Parsers.memcmp(ptr, cmtptr, cmtlen)
        if match
            pos += cmtlen
            pos > len && break
            @inbounds b = buf[pos]
            while b != UInt8('\n') && b != UInt8('\r')
                pos += 1
                pos > len && break
                @inbounds b = buf[pos]
            end
            pos += 1
        else
            break
        end
        ptr = pointer(buf, pos)
    end
    return pos
end

struct ByteValueCounter
    counts::Vector{Int64}
    ByteValueCounter() = new(zeros(Int64, 256))
end

function incr!(c::ByteValueCounter, b::UInt8)
    @inbounds c.counts[b] += 1
    return
end

function guessnrows(buf, oq::UInt8, cq::UInt8, eq::UInt8, source, delim, comment, ignoreemptylines, debug)
    len = fs = length(buf)
    pos = 1
    nbytes = 0
    smallestline = nothing
    lastbytenewline = false
    nlines = 0
    bvc = ByteValueCounter()
    b = 0x00
    pos = consumecommentedline!(buf, pos, len, comment, ignoreemptylines)
    while pos <= len && nlines < 10
        @inbounds b = buf[pos]
        pos += 1
        nbytes += 1
        if b == oq
            while pos <= len
                @inbounds b = buf[pos]
                pos += 1
                nbytes += 1
                if b == eq
                    if pos > len
                        break
                    elseif eq == cq && buf[pos] != cq
                        break
                    end
                    @inbounds b = buf[pos]
                    pos += 1
                    nbytes += 1
                elseif b == cq
                    break
                end
            end
        elseif b == UInt8('\n')
            consumecommentedline!(buf, pos, len, comment, ignoreemptylines)
            nlines += 1
            lastbytenewline = true
            smallestline = min(nbytes, something(smallestline, typemax(Int64)))
            nbytes = 0
        elseif b == UInt8('\r')
            pos <= len && buf[pos] == UInt8('\n') && (pos += 1)
            consumecommentedline!(buf, pos, len, comment, ignoreemptylines)
            nlines += 1
            lastbytenewline = true
            smallestline = min(nbytes, something(smallestline, typemax(Int64)))
            nbytes = 0
        else
            lastbytenewline = false
            incr!(bvc, b)
        end
    end
    nlines += !lastbytenewline
    if nbytes > 0
        smallestline = min(nbytes, something(smallestline, typemax(Int64)))
    end

     if delim === nothing
        if isa(source, AbstractString) && endswith(source, ".tsv")
            d = UInt8('\t')
        elseif isa(source, AbstractString) && endswith(source, ".wsv")
            d = UInt8(' ')
        elseif nlines > 1
            d = nothing
            for attempted_delim in (UInt8(','), UInt8('\t'), UInt8(' '), UInt8('|'), UInt8(';'), UInt8(':'))
                debug && @show Char(attempted_delim)
                debug && @show bvc.counts[Int(attempted_delim)]
                debug && @show nlines
                cnt = bvc.counts[Int(attempted_delim)]
                if cnt > 0 && cnt % nlines == 0
                    d = attempted_delim
                    break
                end
            end
            d = something(d, UInt8(','))
        else
            d = UInt8(',')
        end
    else
        d = (delim isa Char && isascii(delim)) ? delim % UInt8 :
            (sizeof(delim) == 1 && isascii(delim)) ? delim[1] % UInt8 : delim
    end
    smallestline = something(smallestline, 0)
    debug && println("smallest line detected with $smallestline bytes")
    guess = fs / smallestline
    rowsguess = isfinite(guess) ? ceil(Int, guess) : 0
    return rowsguess, d
end


text:
export TabularDeterministicPolicy

"""
    TabularDeterministicPolicy(table::Vector{Int}, nactions::Int) -> TabularDeterministicPolicy

Always select the specific action according to the state defined in `table`.
`nactions` is used to define the number of possible actions.
"""
Base.@kwdef struct TabularDeterministicPolicy <: AbstractPolicy
    table::Vector{Int}
    nactions::Int
end

(p::TabularDeterministicPolicy)(obs) = p(get_state(obs))
(p::TabularDeterministicPolicy)(s::Int) = p.table[s]

RLBase.get_prob(p::TabularDeterministicPolicy, s::Int, a::Int) = p.table[s] == a ? 1.0 : 0.0

function RLBase.update!(p::TabularDeterministicPolicy, correction)
    s, a = correction
    p.table[s] = a
end