text:
isempty(ARGS) && error("Please pass a config file as command line argument.")
config_file = ARGS[1]
isfile(config_file) || error("Cannot read '$config_file'")
println("Config supplied: '$config_file'")
include(config_file)

using SVDD, OneClassActiveLearning, OneClassSampling
using JuMP
using Gurobi
using MLKernels
using Memento
using Random
using Serialization

function setup_experiment_folder(config_file)
    isdir(exp_dir) || mkpath(exp_dir)
    mkpath(joinpath(exp_dir, "log", "experiment"))
    mkpath(joinpath(exp_dir, "log", "worker"))
    mkpath(joinpath(exp_dir, "results"))
    cp(config_file, joinpath(exp_dir, basename(config_file)))
    cp(joinpath(first(splitdir(config_file)), "config.jl"), joinpath(exp_dir, "config.jl"))
end

function create_experiments(sampling_strategies::Vector{Dict{Symbol, Any}}, solver::JuMP.OptimizerFactory, quality_metrics, model_init_strategy::Dict{Symbol, Any})::Vector{Dict{Symbol, Any}}
    experiment_configurations = [
        (data_file, model_params_by_filename[basename(data_file)], s[:method], s[:seed])
        for s in sampling_strategies
        for data_file in data_files
    ]
    all_experiments = []
    for (data_file, model_params, sampling_strategy, seed) in experiment_configurations
        @show (data_file, model_params, sampling_strategy, seed)
        num_samples = model_params[:n_observations]

        data_set_name = splitdir(splitdir(data_file)[1])[2]
        local_data_file_name = splitext(basename(data_file))[1]

        output_path = joinpath(exp_dir, "results", data_set_name)
        isdir(output_path) || mkpath(output_path)

        experiment = Dict{Symbol, Any}(
            :data_file => data_file,
            :data_set_name => data_set_name,
            :log_dir => joinpath(exp_dir, "log"),
            :models => Dict{Symbol, Any}(
                :svdd => Dict{Symbol, Any}(
                    :type => model_params[:svdd][:type],
                    :init_strategies => Dict{Symbol, Any}(
                        :svdd => SVDD.SimpleCombinedStrategy(SVDD.FixedGammaStrategy(MLKernels.GaussianKernel(model_params[:svdd][:gamma])), SVDD.FixedCStrategy(model_params[:svdd][:C])),
                        [Symbol("svdd_reinit_$k") => v for (k, v) in model_init_strategy[:svdd]]...
                    ),
                    :solver => solver),
                :pwc => model_params[:pwc],
                :sample_qe => Dict{Symbol, Any}(
                    :gamma => model_params[:pwc][:gamma],
                    :threshold_strategies => model_init_strategy[:sample_qe][:threshold_strategies]
                )
            ),
            :sampling_strategy => deepcopy(sampling_strategy),
            :quality_metrics => quality_metrics,
            :seed => seed
        )
        exp_hash = hash(experiment)
        experiment[:hash] = "$(exp_hash)"
        experiment[:output_file] = joinpath(output_path, "$(local_data_file_name)_$(typeof(sampling_strategy))_VanillaSVDD_$(exp_hash).json")
        push!(all_experiments, experiment)
    end
    @info "Created $(realpath(exp_dir)) with $(length(all_experiments)) experiment settings."
    return all_experiments
end

function save_experiments(all_experiments::Array{Dict{Symbol, Any},1})
    exp_hashes_file = joinpath(exp_dir, "experiment_hashes")
    @info "Saving experiment hashes at $exp_hashes_file..."
    open(exp_hashes_file, "a") do f
        for e in all_experiments
            write(f, "$(e[:hash])\n")
        end
    end
    @info "... done."
    exp_filename = joinpath(exp_dir, "experiments.jser")
    serialize(exp_filename, all_experiments)
    @info "Experiments written to file $(exp_filename)."
end

# setup experiment directory
exp_dir = joinpath(data_output_root, experiment_name)
if isdir(exp_dir)
    print("Type 'yes' or 'y' to delete and overwrite experiment $(exp_dir): ")
    argin = readline()
    if argin == "yes" || argin == "y"
        rm(exp_dir, recursive=true)
    else
        error("Aborting...")
    end
end
mkpath(exp_dir)
exp_dir = realpath(exp_dir)

# find all data files
all(isdir.(joinpath.(data_input_root, data_dirs))) || error("Not all data dirs are valid.")
data_files = vcat(([joinpath.(data_input_root, x, readdir(joinpath(data_input_root, x))) for x in data_dirs])...)
@info "Found $(length(data_files)) data files."

# parse precompute parameters
isfile(data_parameter_file) || error("Precomputed parameters not found. Please follow the instructions in the Readme.")
@info "Using precompute parameters from file \"$data_parameter_file\"."
model_params_by_filename = deserialize(data_parameter_file)

# generate experiments
setup_experiment_folder(config_file)
Random.seed!(0)
experiments = create_experiments(sampling_strategies, solver, quality_metrics, model_init_strategy)
save_experiments(experiments)

max_stars_repo_name: englhardt/ocs-evaluation

max_issues_repo_name: englhardt/ocs-evaluation

max_forks_repo_name: englhardt/ocs-evaluation


text:
type UnassignedError   <: Exception; end
type AssignedError     <: Exception; end
type SizeMismatchError <: Exception; end

#wire.jl - defines the Wire type.  Wires are stored under the hood as bitvectors.
doc"""
  `Wire{R}`

  is the basic type for Verilog operations.  R specifies a "unit range" of integers.
  Use the "v" suffix to enable verilog-style ranging.

  Wire{3:0v} declares a four-digit verilog wire with indices spanning from 0->3;
  Wire{6:2v} declares a five-bit verilog wire with indices spanning from
  2->6.
"""

immutable Wire{R}
  values::BitVector
  assigned::BitVector

  function Wire(bv1, bv2)
    isa(R, VerilogRange) || throw(TypeError(:Wire, "specifier must be a VerilogRange", VerilogRange, typeof(R)))
    (R.start <= R.stop) || throw(TypeError(:Wire, "constructor range direction failed", R, "backwards"))
    (length(bv1) == length(bv2) == length(R)) || throw(SizeMismatchError())
    new(bv1, bv2)
  end
end

################################################################################
## Aliased naked constructors.
(::Type{Wire})(bv::Bool)               = Wire{0:0v}(BitArray([bv]),trues(1))
(::Type{Wire})(N::Signed)              = Wire{(N-1):0v}(BitVector(N), falses(N))
(::Type{Wire})(R::VerilogRange)        = Wire{R}(BitVector(length(R)), falses(length(R)))
(::Type{Wire})(bv::BitVector)          = Wire{(length(bv)-1):0v}(bv, trues(length(bv)))
#allow initialization of wire with an array, but remember to reverse it.
(::Type{Wire})(wa::Vector{Wire})       = Wire(vcat(map((w) -> w.values, reverse(wa))...))
(::Type{Wire}){N}(wa::Vector{Wire{N}}) = Wire(vcat(map((w) -> w.values, reverse(wa))...))
(::Type{Wire})(ws::Wire...)            = Wire(collect(Wire, ws))

#declaration with an unsigned integer
function (::Type{Wire})(N::Unsigned, l::Integer = 0)
  #override an unspecified length.
  l = (l == 0) ? sizeof(N) * 8 : l
  #mask out crap we don't want.
  Wire(N, range(l))
end
function (::Type{Wire})(N::Unsigned, r::VerilogRange)
  Wire{r}(N)
end
function (::Type{Wire{R}}){R}(N::Unsigned)
  #instantiate a bitarray.
  l = length(R)
  ba = BitVector(length(R))

  N = (l < 64) ? (UInt64(N) & ((1 << l) - 1)) : UInt64(N)
  ba.chunks[1] = N
  #pass this to the bitarray-based constructor.
  Wire{R}(ba, trues(l))
end

function (::Type{Wire{R}}){R}()
  Wire{R}(falses(length(R)), falses(length(R)))
end

type UnsignedBigInt <: Unsigned
  value::BitArray
end

################################################################################
# conversion away from wires - necessary for integer reintepretation of verilog
# wire definitons

function Base.convert{R}(::Type{Unsigned}, w::Wire{R})
  if length(R) <= 64
    return w.values.chunks[1]
  else
    return UnsignedBigInt(w.values)
  end
end

# also allow conversion from wire tuples - will be automatically triggered during
# integer reinterpretations.

function Base.convert(::Type{Unsigned}, wt::Tuple)::Unsigned
  for el in wt
    !isa(el, Wire) && throw(ArgumentError("Tuple must be a wire tuple"))
  end
  Wire(wt...)  #expand the tuple.  will be autoconverted to unsigned.
end

################################################################################

#useful helper functions
import Base: length, range

length{R}(w::Wire{R}) = length(R)
range{R}(w::Wire{R}) = R
assigned(w::Wire) = (&)(w.assigned...)

################################################################################
# getters and setters
import Base: getindex, setindex!

function getindex{R}(w::Wire{R}, n::Integer)
  (n in R) || throw(BoundsError(w, n))
  #adjust for array indexing.
  access_idx = n + 1 - R.start
  #gets the relevant index, if it's been defined.
  w.assigned[access_idx] || throw(UnassignedError())
  Wire(w.values[access_idx])
end

getindex{R}(w::Wire{R}, ::Type{msb}) = getindex(w, R.stop)
getindex{R}(w::Wire{R}, ridx::msb)   = getindex(w, R.stop - ridx.value)

function getindex{R}(w::Wire{R}, r::VerilogRange)
  #returns a wire with the relevant selected values.
  issubset(r, R) || throw(BoundsError(w, r))
  rr = ((r.stop >= r.start) ? (r.start:r.stop) : (r.start:-1:r.stop))
  (&)(w.assigned[rr + 1 - R.start]...) || throw(UnassignedError())
  Wire(w.values[rr + 1 - R.start])
end

getindex{R}(w::Wire{R}, r::RelativeRange) = getindex(w, parse_msb(r, R))

################################################################################
## setters

function setindex!{R}(dst::Wire{R}, src::Wire{0:0v}, n::Integer)
  (n in R) || throw(BoundsError(dst, n))
  offset_idx = n - R.start + 1

  #chcek that the src value exists.
  src.assigned[1] || throw(UnassignedError())
  dst.assigned[offset_idx] && throw(AssignedError())
  dst.assigned[offset_idx] = true
  dst.values[offset_idx] = src.values[1]
  nothing
end

setindex!{R}(dst::Wire{R}, src::Wire{0:0v}, ::Type{msb}) = setindex!(dst, src, R.stop)
setindex!{R}(dst::Wire{R}, src::Wire{0:0v}, m::msb) = setindex!(dst, src, R.stop - m.value)

#you can dereference things as stepranges, but you can't dereference things
#as stepranges.
function Base.setindex!{RD, RS}(dst::Wire{RD}, src::Wire{RS}, r::VerilogRange)
  #check for size mismatch.
  (r.stop >= r.start) || throw(ArgumentError("only forward VerilogRanges allowed for setting"))
  (length(r) == length(RS)) || throw(SizeMismatchError())
  (issubset(r, RD)) || throw(BoundsError(dst, r))

  #the range offset to where they're actually stored in the destination array
  offset_range = r - RD.start + 1
  for idx in 1:length(r)
    dst.assigned[offset_range[idx]] && throw(AssignedError())
    src.assigned[idx]               || throw(UnassignedError())
  end

  for idx in 1:length(r)
    dst.assigned[offset_range[idx]] = true
    dst.values[offset_range[idx]]   = src.values[idx]
  end

  nothing
end

Base.setindex!{RD, RS}(dst::Wire{RD}, src::Wire{RS}, r::RelativeRange) = setindex!(dst, src, parse_msb(r, RD))

#it's useful to declare a single wire shorthand
typealias SingleWire Wire{0:0v}

#=
typealias OptionalWire{R}    Union{Void, Wire{R}}
typealias OptionalSingleWire Union{Void, SingleWire}
=#

export Wire, SingleWire, OptionalWire, OptionalSingleWire

ext: jl


text:
function test_excite()
    println("test excitation generation")
    srand(98765)
    dummy_input = rand(1024)

    println("-- test_excite")
    for hopsize in [40, 80, 160, 320]
        for g in [true, false]
            ex = excite(dummy_input, hopsize, gaussian=g)
            @test all(isfinite.(ex))
        end
    end
end

test_excite()


text:
###########################################################################
#  Riemannian manifold Hamiltonian Monte Carlo (RMHMC)
#  The RMHMC sampler is work in progress (not running yet)
#
#  Reference: Girolami M, Calderhead B. Riemann Manifold Langevin and Hamiltonian Monte Carlo Methods. Journal of the
#  Royal Statistical Society: Series B (Statistical Methodology), 2011, 73 (2), pp 123–214
#
#  Parameters :
#    - nLeaps : number of leapfrog steps
#    - leapStep : leapfrog step size
#    - nNewton: number of Newton steps
#    - tuner: for tuning the RMHMC parameters
#
###########################################################################

export RMHMC

println("Loading RMHMC(nLeaps, leapStep, nNewton, tuner) sampler")

###########################################################################
# RMHMC specific tuners
###########################################################################
abstract RMHMCTuner <: MCMCTuner

###########################################################################
# RMHMC type
###########################################################################
immutable RMHMC <: MCMCSampler
  nLeaps::Int
  leapStep::Float64
  nNewton::Int
  tuner::Union(Nothing, MCMCTuner)

  function RMHMC(nLeaps::Int, leapStep::Float64, nNewton::Int, tuner::Union(Nothing, MCMCTuner))
    @assert nLeaps>0 "Number of leapfrog steps should be > 0"
    @assert leapStep>0 "Leapfrog step size should be > 0"
    @assert nNewton>0 "Number of Newton steps should be > 0"
    new(nLeaps, leapStep, nNewton, tuner)
  end
end

RMHMC(tuner::Union(Nothing, MCMCTuner)=nothing) = RMHMC(6, 0.5, 4, tuner)
RMHMC(nLeaps::Int, tuner::Union(Nothing, MCMCTuner)=nothing) = RMHMC(nLeaps, 3/nLeaps, 4, tuner)
RMHMC(nLeaps::Int, leapStep::Float64, tuner::Union(Nothing, MCMCTuner)=nothing) = RMHMC(nLeaps, leapStep, 4, tuner)
RMHMC(nLeaps::Int, nNewton::Int, tuner::Union(Nothing, MCMCTuner)=nothing) = RMHMC(nLeaps, 3/nLeaps, nNewton, tuner)
RMHMC(leapStep::Float64, tuner::Union(Nothing, MCMCTuner)=nothing) = RMHMC(int(floor(3/leapStep)), leapStep, 4, tuner)
RMHMC(leapStep::Float64, nNewton::Int, tuner::Union(Nothing, MCMCTuner)=nothing) =
  RMHMC(int(floor(3/leapStep)), leapStep, nNewton, tuner)

###########################################################################
# RMHMC sampler
###########################################################################
function SamplerTask(model::MCMCModel, sampler::RMHMC, runner::MCMCRunner)
  # Define local variables
  local pars, proposedPars, leapPars
  local logTarget, proposedLogTarget
  local grad, proposedGrad, leapGrad
  local H, proposedH
  local G, invG, cholG, dG, invGxdG, traceInvGxdG
  local momentum, leapMomentum, momentumTerm, invGMomentum01, invGMomentum02
  local timeStep, nRandomLeaps
  local nLeaps, leapStep

  @assert hasgradient(model) "RMHMC sampler requires model with gradient function"
  @assert hastensor(model) "RMHMC sampler requires model with tensor function"
  @assert hasdtensor(model) "RMHMC sampler requires model with function of tensor derivatives"

  # Allocate memory for some of the local variables
  invGxdG = Array(Float64, model.size, model.size, model.size)
  traceInvGxdG = Array(Float64, model.size)
  momentumTerm = Array(Float64, model.size)

   #  Task reset function
  function reset(resetPars::Vector{Float64})
    pars = copy(resetPars)
    logTarget, grad, G, dG = model.evalalldt(pars)
  end
  # Hook inside Task to allow remote resetting
  task_local_storage(:reset, reset)

  # Initialization
  pars = copy(model.init)
  logTarget, grad = model.evalallg(pars)
  @assert isfinite(logTarget) "Initial values out of model support, try other values"

  if isa(sampler.tuner, EmpMCTuner); tune = EmpiricalHMCTune(sampler.nLeaps, sampler.leapStep, 0, 0); end

  i = 1
  while true
    if isa(sampler.tuner, EmpMCTuner)
      tune.proposed += 1
      nLeaps, leapStep = tune.nLeaps, tune.leapStep
    else
      nLeaps, leapStep = sampler.nLeaps, sampler.leapStep
    end

    proposedPars = copy(pars)

    # Calculate metric tensor G, its inverse invG and its Cholesky factor cholG
    G = model.evalt(proposedPars)
    invG = inv(G)
    cholG = chol(G)

    # Calculate Hamiltonian
    momentum = cholG'*randn(model.size)
    # 0.5*(log(2)+model.size*log(pi)+2*sum(log(diag(cholG)))) is the log determinant of metric tensor G
    H = -logTarget+0.5*(log(2)+model.size*log(pi)+2*sum(log(diag(cholG))))+(momentum'*(invG*momentum))/2

    # Calculate derivative dG of metric tensor G, as well ass invG*dG and trace(invG*dG)
    dG = model.evaldt(proposedPars)
    for j = 1:model.size
      invGxdG[:, :, j] = invG*dG[:, :, j]
      traceInvGxdG[j] = trace(invGxdG[:, :, j])
    end

    # Perform leapfrog steps
    timeStep = (randn() > 0.5 ? 1. : -1.)
    nRandomLeaps = ceil(rand()*nLeaps)

    for j = 1:nRandomLeaps
      leapGrad = model.evalg(proposedPars)
      leapMomentum = copy(momentum)
      for k = 1:sampler.nNewton
        invGMomentum01 = invG*leapMomentum
        for r = 1:model.size
          momentumTerm[r] = (0.5*(leapMomentum'*invGxdG[:, :, r]*invGMomentum01))[1]
        end
        leapMomentum = (momentum+timeStep*(leapStep/2)*(leapGrad-0.5*traceInvGxdG+momentumTerm))
      end

      momentum = copy(leapMomentum)
      invGMomentum02 = invG*momentum
      leapPars = copy(proposedPars)
      for k = 1:sampler.nNewton
        G = model.evalt(leapPars)
        invGMomentum01 = G\momentum
        leapPars = proposedPars+(timeStep*(leapStep/2))*(invGMomentum01+invGMomentum02)
      end

      proposedPars = copy(leapPars)
      G = model.evalt(proposedPars)
      invG = inv(G)
      dG = model.evaldt(proposedPars)
      for k = 1:model.size
        invGxdG[:, :, k] = invG*dG[:, :, k]
        traceInvGxdG[k] = trace(invGxdG[:, :, k])
      end
      invGMomentum01 = invG*momentum
      for k = 1:model.size
        momentumTerm[k] = (0.5*(momentum'*invGxdG[:, :, k]*invGMomentum01))[1]
      end

      proposedGrad = model.evalg(proposedPars)
      momentum = (momentum+timeStep*(leapStep/2)*(proposedGrad-0.5*traceInvGxdG+momentumTerm))
    end

    # Calculate log-target based on the proposed parameters
    proposedLogTarget = model.eval(proposedPars)

    # Calculate Hamiltonian based on the proposed parameters
    proposedH =
      -proposedLogTarget+0.5*(log(2)+model.size*log(pi)+2*sum(log(diag(chol(G)))))+(momentum'*(invG*momentum))/2

    ratio = (H-proposedH)[1]

    if ratio > 0 || (ratio > log(rand()))  # i.e. if accepted
      produce(MCMCSample(proposedPars, proposedLogTarget, proposedGrad, pars, logTarget, grad, {"accept" => true}))
      pars, logTarget, grad = copy(proposedPars), copy(proposedLogTarget), copy(proposedGrad)
      if isa(sampler.tuner, EmpMCTuner); tune.accepted += 1; end
    else
      produce(MCMCSample(pars, logTarget, grad, pars, logTarget, grad, {"accept" => false}))
    end

    if isa(sampler.tuner, EmpMCTuner) && i<= runner.burnin && mod(i, sampler.tuner.adaptStep) == 0
      adapt!(tune, sampler.tuner)
      reset!(tune)
      if sampler.tuner.verbose
        println("Burn-in teration $i of $(runner.burnin): ", round(100*tune.rate, 2), " % acceptance rate")
      end
    end

    i += 1
  end
end


text:
"""
    register!([:datatype|:predicate|:function|:modifier], name, x)

Register `x` as a global datatype, predicate, function or modifier under
the specified `name`.
"""
function register!(category::Symbol, name::Symbol, x)
    if category == :datatype
        register_datatype!(name, x)
    elseif category == :predicate
        register_predicate!(name, x)
    elseif category == :function
        register_function!(name, x)
    elseif category == :modifier
        register_modifier!(name, x)
    else
        error("Category must be :datatype, :predicate, :function or :modifier.")
    end
end

register!(category::Symbol, name::AbstractString, x) =
    register!(category, Symbol(name), x)

"Register global datatype."
function register_datatype!(name::Symbol, ty)
    GLOBAL_DATATYPES[name] = ty
end

"Register global predicate."
function register_predicate!(name::Symbol, f)
    GLOBAL_PREDICATES[name] = f
    GLOBAL_FUNCTIONS[name] = f
end

"Register global function."
function register_function!(name::Symbol, f)
    GLOBAL_FUNCTIONS[name] = f
end

"Register global modifier."
function register_modifier!(name::Symbol, f)
    GLOBAL_MODIFIERS[name] = f
end

"""
    deregister!([:datatype|:predicate|:function|:modifier], name)

Deregister the datatype, predicate, function or modifier specified by `name`.
"""
function deregister!(category::Symbol, name::Symbol)
    if category == :datatype
        deregister_datatype!(name)
    elseif category == :predicate
        deregister_predicate!(name)
    elseif category == :function
        deregister_function!(name)
    elseif category == :modifier
        deregister_modifier!(name)
    else
        error("Category must be :datatype, :predicate, :function or :modifier.")
    end
end

deregister!(category::Symbol, name::AbstractString) =
    deregister!(category, Symbol(name))

"Deregister global datatype."
function deregister_datatype!(name::Symbol)
    delete!(GLOBAL_DATATYPES, name)
end

"Deregister global predicate."
function deregister_predicate!(name::Symbol)
    delete!(GLOBAL_FUNCTIONS, name)
    delete!(GLOBAL_PREDICATES, name)
end

"Deregister global function."
function deregister_function!(name::Symbol)
    delete!(GLOBAL_PREDICATES, name)
    delete!(GLOBAL_FUNCTIONS, name)
end

"Deregister global modifier."
function deregister_modifier!(name::Symbol)
    delete!(GLOBAL_MODIFIERS, name)
end

"""
    attach!(domain, :function, name, f)

Attach the function `f` as the implementation of the functional fluent
specified by `name`.

    attach!(domain, :datatype, name, ty)

Attach the type `ty` as the implementation of the PDDL datatype
specified by `name`.
"""
function attach!(domain::GenericDomain, category::Symbol, name::Symbol, f)
    if category == :datatype
        attach_datatype!(domain, name, f)
    elseif category == :function
        attach_function!(domain, name, f)
    else
        error("Category must be :datatype or :function.")
    end
end
attach!(domain::GenericDomain, category::Symbol, name::AbstractString, f) =
    attach!(domain, category, Symbol(name), f)

"Attach datatype to domain."
function attach_datatype!(domain::GenericDomain, name::Symbol, ty)
    domain.datatypes[name] = ty
end

"Attach function to domain."
function attach_function!(domain::GenericDomain, name::Symbol, f)
    if name in keys(get_functions(domain))
        domain.funcdefs[name] = f
    else
        error("Domain does not have a function named $name.")
    end
end


text:
using Test
using FreeParameters
using Distributions

output = joinpath(@__DIR__,"..","output")
run(`rm -rf $output`)

include("macro.jl")

include("dycomsmodel.jl")
include("define_model.jl")

include("struct_2_dict.jl")
include("extract.jl")
include("export_dict.jl")
include("import_dict.jl")
include("instantiate.jl")


text:
module DaemonMode

export serve, runfile, runargs, runexpr, sendExitCode

using Logging, LoggingExtras

using Crayons
using Crayons.Box

using Sockets

const PORT = 3000

function add_packages(fname::AbstractString)
end

const first_time = Ref(true)
const token_runfile = "DaemonMode::runfile"
const token_runexpr = "DaemonMode::runexpr"
const token_exit = "DaemonMode::exit"
const token_client_end = "DaemonMode::client_end"
const token_server_async_end = "DaemonMode::async_end"
const token_server_sync_end = "DaemonMode::sync_end"
const token_ok_end = "DaemonMode::end_ok"
const token_error_end = "DaemonMode::end_er"
const token_end = "DaemonMode::end"

"""
    serve_sync(port=3000, shared=false, print_stack=true, async=false)

Run the daemon, running all files and expressions sended by the client function.

# Optionals

- port: port to listen (default=3000).
- shared: Share the environment between calls. If it is false (default) each run
  has its own environment, so the variables/functions are not shared.
- print_stack: Print the complete stack when there is an error. By default it is true.
- async: Run the clients in different clients, **Experimental**.
"""
function serve(port=PORT, shared=missing; print_stack=true, async=true)
    if (async)
        return serve_async(port, shared, print_stack=print_stack)
    else
        return serve_sync(port, shared; print_stack=print_stack)
    end

end

function serve_async(port=PORT, shared=missing; print_stack=true, nicely=true)
    server = Sockets.listen(Sockets.localhost, port)
    current = pwd()
    continue_server = Threads.Atomic{Bool}(true)
    tasks = Task[]

    while continue_server[] && isopen(server)
        sock = accept(server)
        task = @async begin
            local mode
            mode = readline(sock)

            if mode == token_runfile
                serverRunFile(sock, coalesce(shared, false), print_stack)
            elseif mode == token_runexpr
                serverRunExpr(sock, coalesce(shared, true), print_stack)
            elseif mode == token_client_end
                println(sock, token_server_async_end)
                continue_server[] = false
            else
                println(sock, "Error, unrecognised mode, expected (\"runFile()\", \"runExpr()\" or \"exit()\", but received \"$mode\"")
                continue_server[] = false
            end
        end

        push!(tasks, task)
    end

    if nicely
        # wait all pending tasks
        for task in tasks
            try
                wait(task)
            catch e
            end
        end
    end

    close(server)
end

function serve_sync(port=PORT, shared=missing; print_stack=true)
    server = Sockets.listen(Sockets.localhost, port)
    current = pwd()
    continue_server = true

    while continue_server
        sock = accept(server)
        mode = readline(sock)

        if mode == token_runfile
            serverRunFile(sock, coalesce(shared, false), print_stack)
        elseif mode == token_runexpr
            serverRunExpr(sock, coalesce(shared, true), print_stack)
        elseif mode == token_client_end
            println(sock, token_server_sync_end)
            continue_server = false
        else
            println(sock, "Error, unrecognised mode, expected (\"runFile()\", \"runExpr()\" or \"exit()\", but received \"$mode\"")
            continue_server = false
        end
    end

    close(server)
end

function myshowerror(sock::IO, e)
    io = IOBuffer()
    showerror(io, e)
    msg = String(take!(io))
    line = first(split(msg, '\n'))
    posi = findfirst(':', line)
    print(sock, RED_FG, BOLD, "ERROR: ")
    print(sock, Crayon(foreground=:white, bold=false), line[1:posi])
    print(sock, Crayon(foreground=:red, bold=false), line[posi+1:end])
    print(sock, Crayon(reset=true))
end

function serverReplyError(sock, e)
    try
        myshowerror(sock, e)
        println(sock)
        println(sock, token_error_end)
    catch e
        if (e isa Base.IOError) && abs(e.code) == abs(Libc.EPIPE)
            # client disconnected early, ignore
        else
            rethrow()
        end
    end
end

function send_backtrace(sock, bt, fname)
    stacks = map(first, Base.process_backtrace(bt))
    fullname = ""

    if !isempty(stacks)
        println(sock)
        println(sock, "Stacktrace:")
    end

    if !isempty(fname)
        fullname = joinpath(pwd(), fname)
    end

    for (i, stack) in enumerate(stacks)
        file = String(stack.file)

        if occursin("string", file) && !isempty(fullname)
            file = fullname
        end

        if occursin("loading.jl", file)
            return
        end

        print(sock, " [$i] ")

        print(sock, BOLD, stack.func)
        print(sock, Crayon(reset=true), " at ")
        print(sock, BOLD, file, ":", stack.line)
        println(sock, Crayon(reset=true))
    end
end

function serverReplyError(sock, e, bt, fname)
    try
        myshowerror(sock, e)
        send_backtrace(sock, bt, fname)
        println(sock)
        println(sock, token_error_end)
    catch e
        if (e isa Base.IOError) && abs(e.code) == abs(Libc.EPIPE)
            # client disconnected early, ignore
        else
            rethrow()
        end
    end
end

function create_mylog(fname)
    function mylog(io, args)
        reset = Crayon(reset=true)
        color = reset
        type = string(args.level)

        if args.level == Logging.Warn
            color = Crayon(foreground=:yellow, bold=true)
            type = "Warning"
        elseif args.level == Logging.Error
            color = Crayon(foreground=:red, bold=true)
        elseif args.level == Logging.Info
            color = Crayon(foreground=:cyan, bold=true)
        end

        print(io, color, "┌ ", type, ": ", reset);
        lines = split(args.message, "\n")
        println(io, lines[1])

        for line in lines[2:end]
            println(io, color, "│ ", reset, line)
        end

        module_str = string(args._module)
        module_str = replace(module_str, ".anonymous" => "")
        file = args.file

        if file == "string"
            file = joinpath(pwd(), fname)
        end

        println(io, color, "└ ", Crayon(foreground=:dark_gray), "@ ", module_str, " ", file, ": ", args.line)
    end
end

function serverRun(run, sock, shared, print_stack, fname, args)
    try
        if shared
            redirect_stdout(sock) do
                redirect_stderr(sock) do
                    run(Main)
                end
            end
        else
            redirect_stdout(sock) do
                redirect_stderr(sock) do
                m = Module()
                Logging.global_logger(MinLevelLogger(FormatLogger(create_mylog(fname), sock), Logging.Info))
                add_include = Meta.parse("include(arg)=Base.include(@__MODULE__,arg)")
                Base.eval(m, add_include)

                if !isempty(args)
                    add_params =  Meta.parse(string("ARGS = [\"", join(args, "\",\""), "\"]"))
                else
                    add_params =  Meta.parse(string("empty!(ARGS)"))
                end

                Base.eval(m, add_params)
                # Following code is not needed, the real problem was global ARGS, not
                add_redirect = Meta.parse("const stdout=IOBuffer(); println(io, x...) = Base.println(io,x...); println(x)=Base.println(stdout, x); println(x...)=Base.println(stdout, x...); println(io, x...)=Base.println(io, x...); print(x...)=Base.print(stdout, x...); stdout")
                out = Base.eval(m, add_redirect)
                add_redirect_err = Meta.parse("const stderr=IOBuffer(); stderr")
                err = Base.eval(m, add_redirect_err)

                task = @async begin
                    while isopen(sock)
                        text = String(take!(out))
                        print(sock, text)
                        sleep(0.5)
                    end
                end
                task2 = @async begin
                    while isopen(sock)
                        text = String(take!(err))
                        print(sock, text)
                        sleep(0.5)
                    end
                end
                # add_logging = Meta.parse("using ExtraLogging; ExtraLogging.global_logger(MinLevelLogger(FormatLogger(create_mylog(fname), stdout), Logging.Info))")
                # Base.eval(m, add_logging)
                run(m)
                # If there is missing message I write it
                text = String(take!(out))

                if !isempty(text)
                    print(sock, text)
                end
                # If there is missing message I write it
                text = String(take!(err))

                if !isempty(text)
                    print(sock, text)
                end
             end
            end
        end
        println(sock, token_ok_end)
    catch e
        if print_stack
            serverReplyError(sock, e, catch_backtrace(), fname)
        else
            serverReplyError(sock, e)
        end
    end

end

"""
serverRunFile(sock, shared)

Run the source code of the filename push through the socket.

# Parameters

- sock: socket in which is going to receive the dir, the filename, and args to run.
- shared: Share the environment between calls. If it is false (default) each run
has its own environment, so the variables/functions are not shared.
"""
function serverRunFile(sock, shared, print_stack)
    fname = ""

    try
        dir = readline(sock)
        fname = readline(sock)
        args_str = readline(sock)
        args = split(args_str, " ")

        if !isempty(args) && isempty(args[1])
            empty!(args)
        end

        # Add it to allow ArgParse and similar packages
        empty!(ARGS)

        for arg in args
            push!(ARGS, arg)
        end

        first_time[] = true

        cd(dir) do
            content = read(fname, String)
            serverRun(sock, shared, print_stack, fname, args) do mod
                include_string(mod, content)
            end
        end
    catch e
        serverReplyError(sock, e)
    end

end

"""
serverRunExpr(sock, shared)

Run the source code of the filename push through the socket.

# Parameters

- sock: socket in which is going to receive the code to run.
- shared: Share the environment between calls. If it is false (default) each run
has its own environment, so the variables/functions are not shared.
"""
function serverRunExpr(sock, shared, print_stack)

    try
        dir = readline(sock)
        expr = readuntil(sock, token_end) # Read until token_end to handle multi-line expressions
        parsedExpr = Meta.parse(strip(expr))

        cd(dir) do
            serverRun(sock, shared, print_stack, "",  String[]) do mod
                Base.eval(mod, parsedExpr)
            end
        end
    catch e
        serverReplyError(sock, e)
    end

end

"""
    runexpr(expr::AbstractString)

Ask the server to run julia code in a string pass as parameters.

# Parameters

- expr: Julia code to run in the server.

# Optionals

- port: Port (default=3000).
- output: stream in which it is shown the output of the run.
"""
function runexpr(expr::AbstractString ; output = stdout, port = PORT)
    try
        sock = Sockets.connect(port)
        println(sock, token_runexpr)
        println(sock, pwd())
        println(sock, expr)
        println(sock, token_end)

        line = readline(sock)
        while (line != token_ok_end && line != token_error_end)
            println(output, line)
            line = readline(sock)
        end

    catch e
        println(stderr, "Error, cannot connect with server. Is it running?")
    end
end

"""
    runfile(fname::AbstractString)

Ask the server to run a specific filename.

# Parameters

- fname: Filename to run.

# Optionals

- args: List of arguments (array of String, default=[]).
- port: Port (default=3000)
- output: stream in which it is shown the output of the run.
"""
function runfile(fname::AbstractString; args=String[], port = PORT, output=stdout)
    result = 0
    dir = dirname(fname)

    if isempty(dir)
        fcompletename = joinpath(dir, fname)
    else
        fcompletename = fname
    end
    try
        sock = Sockets.connect(port)
        println(sock, token_runfile)
        println(sock, pwd())
        println(sock, fcompletename)
        println(sock, join(args, " "))
        line = readline(sock)
        token_size = length(token_ok_end)

        while (length(line) < token_size || !occursin(token_end, line))
            println(output, line)
            line = readline(sock)
        end

        if occursin(token_error_end, line)
            result = 1
        end

        if length(line) > token_size
            end_line = replace(line, token_ok_end => "")
            end_line = replace(end_line, token_error_end => "")
            print(output, end_line)
        end
    catch e
        println(stderr, "Error, cannot connect with server. Is it running?")
        # println(stderr, e)
        exit(1)
    end
    return result
end

"""
    sendExitCode(port)

send the exit code, it closes the server.

# Optionals

- port: port to connect (default=3000).
"""
function sendExitCode(port=PORT)
    sock = Sockets.connect(port)
    println(sock, token_client_end)
    line = readline(sock)
    error = true

    if line == token_server_sync_end
        error = false
    elseif line == token_server_async_end
        error = false

        # For async it is need to create a new connection
        try
            sock = Sockets.connect(port)
        catch e
        end
    end

    return error
end

"""
    runargs(port=PORT)

Ask the server to run all files in ARGS.

# Optionals

- port: Port of the server to ask (default=3000).
"""
function runargs(port=PORT)
    if isempty(ARGS)
        println(stderr, "Error: missing filename")
        exit(1)
    elseif !isfile(ARGS[1])
        println(stderr, "Error: file '$(ARGS[1])' doest not exist")
        exit(1)
    end

    result = runfile(ARGS[1], args=ARGS[2:end], port=port)
    exit(result)
end

end # module


text:
# custom help for byrow()
function Docs.getdoc(x::typeof(byrow), y)
    if y == Union{}
        return _get_doc_byrow("default")
    elseif y == Tuple{typeof(sum)}
        return _get_doc_byrow("sum")
    elseif y == Tuple{typeof(mean)}
        return _get_doc_byrow("mean")
    elseif y == Tuple{typeof(all)}
        return _get_doc_byrow("all")
    elseif y == Tuple{typeof(any)}
        return _get_doc_byrow("any")
    elseif y == Tuple{typeof(count)}
        return _get_doc_byrow("count")
    elseif y == Tuple{typeof(prod)}
        return _get_doc_byrow("prod")
    elseif y == Tuple{typeof(isequal)}
        return _get_doc_byrow("isequal")
    elseif y == Tuple{typeof(isless)}
        return _get_doc_byrow("isless")
    elseif y == Tuple{typeof(in)}
        return _get_doc_byrow("in")
    elseif y == Tuple{typeof(findfirst)}
        return _get_doc_byrow("findfirst")
    elseif y == Tuple{typeof(findlast)}
        return _get_doc_byrow("findlast")
    elseif y == Tuple{typeof(select)}
        return _get_doc_byrow("select")
    elseif y == Tuple{typeof(fill!)}
        return _get_doc_byrow("fill!")
    elseif y == Tuple{typeof(fill)}
        return _get_doc_byrow("fill")
    elseif y == Tuple{typeof(coalesce)}
        return _get_doc_byrow("coalesce")
    elseif y == Tuple{typeof(maximum)}
        return _get_doc_byrow("maximum")
    elseif y == Tuple{typeof(minimum)}
        return _get_doc_byrow("minimum")
    elseif y == Tuple{typeof(argmax)}
        return _get_doc_byrow("argmax")
    elseif y == Tuple{typeof(argmin)}
        return _get_doc_byrow("argmin")
    elseif y == Tuple{typeof(issorted)}
        return _get_doc_byrow("issorted")
    elseif y == Tuple{typeof(join)}
        return _get_doc_byrow("join")
    elseif y == Tuple{typeof(hash)}
        return _get_doc_byrow("hash")
    elseif y == Tuple{typeof(nunique)}
        return _get_doc_byrow("nunique")
    elseif y == Tuple{typeof(mapreduce)}
        return _get_doc_byrow("mapreduce")
    elseif y == Tuple{typeof(var)}
        return _get_doc_byrow("var")
    elseif y == Tuple{typeof(std)}
        return _get_doc_byrow("std")
    elseif y == Tuple{typeof(cumsum!)}
        return _get_doc_byrow("cumsum!")
    elseif y == Tuple{typeof(cumsum)}
        return _get_doc_byrow("cumsum")
    elseif y == Tuple{typeof(cumprod!)}
        return _get_doc_byrow("cumprod!")
    elseif y == Tuple{typeof(cumprod)}
        return _get_doc_byrow("cumprod")
    elseif y == Tuple{typeof(cummax!)}
        return _get_doc_byrow("cummax!")
    elseif y == Tuple{typeof(cummax)}
        return _get_doc_byrow("cummax")
    elseif y == Tuple{typeof(cummin!)}
        return _get_doc_byrow("cummin!")
    elseif y == Tuple{typeof(cummin)}
        return _get_doc_byrow("cummin")
    elseif y == Tuple{typeof(sort!)}
        return _get_doc_byrow("sort!")
    elseif y == Tuple{typeof(sort)}
        return _get_doc_byrow("sort")
    elseif y == Tuple{typeof(stdze!)}
        return _get_doc_byrow("stdze!")
    elseif y == Tuple{typeof(stdze)}
        return _get_doc_byrow("stdze")
    else
        return _get_doc_byrow("generic")
    end
end

function _get_doc_byrow(fun; text = byrow_docs_text)
    split_text = split(text, "@@@@")
    loc = findfirst(==(fun), split_text)
    Markdown.parse(split_text[loc+1])
end

byrow_docs_text = """
@@@@default@@@@
    byrow(ds::AbstractDataset, fun, cols; ...)

Perform a row-wise operation specified by `fun` on selected columns `cols`. Generally,
`fun` can be any function that returns a scalar value for each row.

> User can pass a type as `fun` when `cols` is referring to a single column. In this case, `byrow` simply converts the selected column to vector of type `fun`.

`byrow` is fine tuned for the following operations. To get extra help for each of them search help for `byrow(fun)`, e.g. `?byrow(sum)`;

# Reduction operations

- `all`
- `any`
- `argmax`
- `argmin`
- `coalesce`
- `count`
- `findfirst`
- `findlast`
- `hash`
- `in`
- `isequal`
- `isless`
- `issorted`
- `join`
- `mapreduce`
- `maximum`
- `mean`
- `minimum`
- `nunique`
- `prod`
- `select`
- `std`
- `sum`
- `var`

# Special operations

- `cummax`
- `cummax!`
- `cummin`
- `cummin!`
- `cumprod`
- `cumprod!`
- `cumsum`
- `cumsum!`
- `fill`
- `fill!`
- `sort`
- `sort!`
- `stdze`
- `stdze!`
@@@@sum@@@@
    byrow(ds::AbstractDataset, sum, cols = names(ds, Number); [by = identity, threads])

Sum results of calling function `by` on each element of each row of `ds`. If `cols` is not specified, `byrow`
computes sum for all numeric columns in `ds`.

Passing `threads = false` disables multitrheaded computations.

Missing values are removed from the calculation. When all values in a row are missing, it returns `missing`.

## Example

```jldoctest
julia> ds = Dataset(x = [1,2,3], y = [2.0, 1.5, 4.0])
3×2 Dataset
 Row │ x         y
     │ identity  identity
     │ Int64?    Float64?
─────┼────────────────────
   1 │        1       2.0
   2 │        2       1.5
   3 │        3       4.0

julia> byrow(ds, sum, :)
3-element Vector{Float64}:
 3.0
 3.5
 7.0
```
@@@@mean@@@@
    byrow(ds::AbstractDataset, mean, cols = names(ds, Number); [by = identity, threads])

Compute mean of the results of calling function `by` on each element of each row of `ds`. If `cols` is not specified, `byrow`
computes mean for all numeric columns in `ds`.

Passing `threads = false` disables multitrheaded computations.

Missing values are removed from the calculation. When all values in a row are missing, it returns `missing`.

## Example

```jldoctest
julia> ds = Dataset(x = [1,2,3], y = [2.0, 1.5, 4.0])
3×2 Dataset
 Row │ x         y
     │ identity  identity
     │ Int64?    Float64?
─────┼────────────────────
   1 │        1       2.0
   2 │        2       1.5
   3 │        3       4.0

julia> byrow(ds, mean, :)
3-element Vector{Float64}:
 1.5
 1.75
 3.5
```
@@@@all@@@@
    byrow(ds::AbstractDataset, all, cols = :; [by = isequal(true), threads, mapformats = false])

Test whether all elements in each row in selected columns are `true`, when `by` is passed, determine whether predicate `by` returns `true` for all elements in the row.

By default, `byrow` uses the actual values for test the elements, however, passing `mapformats = true`
change this to the formatted values.

Each columns in `cols` may have its own `by`. This may be achieved by passing a vector of predicates to `by`.

Passing `threads = false` disables multitrheaded computations.

See [`filter`](@ref) or [`filter!`](@ref)

## Example

```jldoctest
julia> ds = Dataset(x = [1,2,3], y = [2.0, 1.5, 4.0])
3×2 Dataset
 Row │ x         y
     │ identity  identity
     │ Int64?    Float64?
─────┼────────────────────
   1 │        1       2.0
   2 │        2       1.5
   3 │        3       4.0

julia> byrow(ds, all, :, by = [==(2), >(1)])
3-element Vector{Bool}:
 0
 1
 0
```
@@@@any@@@@
    byrow(ds::AbstractDataset, any, cols = :; [by = isequal(true), threads, mapformats = false])

Test whether any elements in each row in selected columns is `true`, when `by` is passed, determine whether predicate `by` returns `true` for any elements in the row.

By default, `byrow` uses the actual values for test the elements, however, passing `mapformats = true`
change this to the formatted values.

Each columns in `cols` may have its own `by`. This may be achieved by passing a vector of predicates to `by`.

Passing `threads = false` disables multitrheaded computations.

See [`filter`](@ref) or [`filter!`](@ref)

## Example

```jldoctest
julia> ds = Dataset(x = [1,2,3], y = [2.0, 1.5, 4.0])
3×2 Dataset
 Row │ x         y
     │ identity  identity
     │ Int64?    Float64?
─────┼────────────────────
   1 │        1       2.0
   2 │        2       1.5
   3 │        3       4.0

julia> byrow(ds, any, :, by = [==(2), >(1)])
3-element Vector{Bool}:
 1
 1
 1
```
@@@@count@@@@
    byrow(ds::AbstractDataset, count, cols = :; [by = isequal(true), threads])

Count the number of elements in each row for selected columns which the function `by` returns `true`.

Passing `threads = false` disables multitrheaded computations.

## Example

```jldoctest
julia> julia> ds = Dataset(x = [1,2,3], y = [2, 6, 5])
3×2 Dataset
 Row │ x         y
     │ identity  identity
     │ Int64?    Int64?
─────┼────────────────────
   1 │        1         2
   2 │        2         6
   3 │        3         5

julia> byrow(ds, count, :, by = isodd)
3-element Vector{Int32}:
 1
 0
 2
```
@@@@prod@@@@
    byrow(ds::AbstractDataset, prod, cols = names(ds, Number); [by = identity, threads])

Return the product of the results of calling function `by` on each element of each row of `ds`. If `cols` is not specified, `byrow`
computes product for all numeric columns in `ds`.

Passing `threads = false` disables multitrheaded computations.

Missing values are removed from the calculation. When all values in a row are missing, it returns `missing`.

## Example

```jldoctest
julia>  ds = Dataset(x = [1,2,3], y = [2.0, 1.5, 4.0])
3×2 Dataset
 Row │ x         y
     │ identity  identity
     │ Int64?    Float64?
─────┼────────────────────
   1 │        1       2.0
   2 │        2       1.5
   3 │        3       4.0

julia> byrow(ds, prod, :)
3-element Vector{Union{Missing, Float64}}:
  2.0
  3.0
 12.0
```
@@@@isequal@@@@
    byrow(ds::AbstractDataset, isequal, cols; [with = nothing, threads])

Returns a boolean vector which is `true` if all values in the corresponding row are equal (using `isequal`).
 Optionally, a vector of values can be passed view the `with` keyword argument to compare values in selected
 columns with the passed vector.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(isless)`](@ref), [`byrow(in)`](@ref), [`byrow(issorted)`](@ref)

## Examples

```jldoctest
julia> ds = Dataset(x1 = [1,2,3,1,2,3], x2 = [1,2,1,2,1,2])
6×2 Dataset
 Row │ x1        x2
     │ identity  identity
     │ Int64?    Int64?
─────┼────────────────────
   1 │        1         1
   2 │        2         2
   3 │        3         1
   4 │        1         2
   5 │        2         1
   6 │        3         2

julia> byrow(ds, isequal, [1,2])
6-element Vector{Bool}:
 1
 1
 0
 0
 0
 0

julia> byrow(ds, isequal, [1,2], with = [2,2,2,3,3,3])
6-element Vector{Bool}:
 0
 1
 0
 0
 0
 0
```
@@@@isless@@@@
    byrow(ds::AbstractDataset, isless, cols, [with, threads, rev = false, lt = isless])

Return a boolean vector which is true if all values in corresponding row for selected `cols` are less than value given by the `with` keyword argument. A vector, or a column name can be passed via `with`.

Passing `rev = true` returns true if all values are greater than passed values via `with`.

By default, the comparison is done via `isless` function, however, user may change it by passing a function via the `lt` keyword argument. The function passed to `lt` must accept two arguments where it takes its first argument from `cols` and its second argument from `with`. However, if `rev = true` the function passed as `lt` will take its first argument from `with` and its second argument from `cols`. The function passed as `lt` must return `true` or `false`.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(isequal)`](@ref), [`byrow(in)`](@ref), [`byrow(issorted)`](@ref)

## Examples

```jldoctest
julia> ds = Dataset(x1 = [1,2,3,1,2,3], x2 = [1,2,1,2,1,2], x3 = 6:-1:1)
6×3 Dataset
 Row │ x1        x2        x3
     │ identity  identity  identity
     │ Int64?    Int64?    Int64?
─────┼──────────────────────────────
   1 │        1         1         6
   2 │        2         2         5
   3 │        3         1         4
   4 │        1         2         3
   5 │        2         1         2
   6 │        3         2         1

julia> byrow(ds, isless, [1,2], with = :x3)
6-element Vector{Bool}:
 1
 1
 1
 1
 0
 0

julia> byrow(ds, isless, 1:2, with = :x3, lt = (x,y) -> isless(x^2, y))
6-element Vector{Bool}:
 1
 1
 0
 0
 0
 0

julia> ds = Dataset(x1 = [1,2,3,1,2,3], x2 = [1,2,1,2,1,2],
                    x3 = [(1,2), (2,3), (1,2), (4,5), (1,2), (1,2)])
6×3 Dataset
 Row │ x1        x2        x3
     │ identity  identity  identity
     │ Int64?    Int64?    Tuple…?
─────┼──────────────────────────────
   1 │        1         1  (1, 2)
   2 │        2         2  (2, 3)
   3 │        3         1  (1, 2)
   4 │        1         2  (4, 5)
   5 │        2         1  (1, 2)
   6 │        3         2  (1, 2)

julia> byrow(ds, isless, 1:2, with = :x3, lt = in)
6-element Vector{Bool}:
 1
 1
 0
 0
 1
 0
```
@@@@in@@@@
    byrow(ds::AbstractDataset, in, cols; [item, threads, eq = isequal])

Return a boolean vector which its elements are true if in a row the value of `item` is equal to any values from `cols`. The equality is checked via the function passed as `eq`. User can pass a vector of values or a column name to `item`.

The function passed as `eq` must accept two arguments where it takes its first argument from `item` and its second argument from `cols`. The function passed as `eq` must return `true` or `false`.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(isequal)`](@ref), [`byrow(isless)`](@ref), [`byrow(issorted)`](@ref)

## Examples

```jldoctest
julia> ds = Dataset(x1 = [1,2,3,1,2,3], x2 = [1,2,1,2,1,2], x3 = 6:-1:1)
6×3 Dataset
 Row │ x1        x2        x3
     │ identity  identity  identity
     │ Int64?    Int64?    Int64?
─────┼──────────────────────────────
   1 │        1         1         6
   2 │        2         2         5
   3 │        3         1         4
   4 │        1         2         3
   5 │        2         1         2
   6 │        3         2         1

julia> byrow(ds, in, r"x", item = [1,2,3,4,5,6])
6-element Vector{Bool}:
 1
 1
 1
 0
 0
 0

julia> byrow(ds, in, [2,3], item = :x1, eq = isless)
6-element Vector{Bool}:
 1
 1
 1
 1
 0
 0
julia> byrow(ds, in, r"x", item = [5,4,5,4,5,4], eq = (x,y) -> x+y == 11)
6-element Vector{Bool}:
 1
 0
 0
 0
 0
 0
```
@@@@findfirst@@@@
    byrow(ds::AbstractDataset, findfirst, cols; [by = identity, item = nothing, eq = isequal, threads])

Return the column name of the first `true` value in `cols` or for which `by` returns `true`. If no such value is found, it returns `missing`. User can pass a vector of values or a column name to `item` to find the column name of the first time that the value of `item` is equal to the value of the column. User may use a customised function for checking the equlity of `item` and `columns` by passing it to the `eq` keyword argument. The function passed as `eq` must be a binary function where its first argument is from `item` and its second argument is from `col`.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(findlast)`](@ref), [`byrow(select)`](@ref)

# Examples
```jldoctest
julia> ds = Dataset(g = [1, 1, 1, 2, 2],
                               x1_int = [0, 0, 1, missing, 2],
                               x2_int = [3, 2, 1, 3, -2],
                               x1_float = [1.2, missing, -1.0, 2.3, 10],
                               x2_float = [missing, missing, 3.0, missing, missing],
                               x3_float = [missing, missing, -1.4, 3.0, -100.0])
5×6 Dataset
 Row │ g         x1_int    x2_int    x1_float   x2_float   x3_float
     │ identity  identity  identity  identity   identity   identity
     │ Int64?    Int64?    Int64?    Float64?   Float64?   Float64?
─────┼───────────────────────────────────────────────────────────────
   1 │        1         0         3        1.2  missing    missing
   2 │        1         0         2  missing    missing    missing
   3 │        1         1         1       -1.0        3.0       -1.4
   4 │        2   missing         3        2.3  missing          3.0
   5 │        2         2        -2       10.0  missing       -100.0

julia> byrow(ds, findfirst, :, by = ismissing)
5-element PooledArrays.PooledVector{Union{Missing, Symbol}, UInt32, Vector{UInt32}}:
 :x2_float
 :x1_float
 missing
 :x1_int
 :x2_float

julia> byrow(ds, findfirst, 1:3, item = [1,1,1,1,1])
5-element PooledArrays.PooledVector{Union{Missing, Symbol}, UInt32, Vector{UInt32}}:
 :g
 :g
 :g
 missing
 missing

julia> ds = Dataset(x1 = [1,2,2], x2 = [5,6,7], x3 = [8,9,10])
3×3 Dataset
 Row │ x1        x2        x3
     │ identity  identity  identity
     │ Int64?    Int64?    Int64?
─────┼──────────────────────────────
   1 │        1         5         8
   2 │        2         6         9
   3 │        2         7        10

julia> byrow(ds, select, :, with = byrow(ds, findfirst, :, by = isodd))
3-element Vector{Union{Missing, Int64}}:
 1
 9
 7
```
@@@@findlast@@@@
    byrow(ds::AbstractDataset, findlast, cols; [by = identity, item = nothing, eq = isequal, threads])

Return the column name of the last `true` value in `cols` or for which `by` returns `true`. If no such value is found, it returns `missing`. User can pass a vector of values or a column name to `item` to find the column name of the last time that the value of `item` is equal to the value of the column. User may use a customised function for checking the equlity of `item` and `columns` by passing it to the `eq` keyword argument. The function passed as `eq` must be a binary function where its first argument is from `item` and its second argument is from `col`.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(findfirst)`](@ref), [`byrow(select)`](@ref)

# Examples
```jldoctest
julia> ds = Dataset(g = [1, 1, 1, 2, 2],
                               x1_int = [0, 0, 1, missing, 2],
                               x2_int = [3, 2, 1, 3, -2],
                               x1_float = [1.2, missing, -1.0, 2.3, 10],
                               x2_float = [missing, missing, 3.0, missing, missing],
                               x3_float = [missing, missing, -1.4, 3.0, -100.0])
5×6 Dataset
 Row │ g         x1_int    x2_int    x1_float   x2_float   x3_float
     │ identity  identity  identity  identity   identity   identity
     │ Int64?    Int64?    Int64?    Float64?   Float64?   Float64?
─────┼───────────────────────────────────────────────────────────────
   1 │        1         0         3        1.2  missing    missing
   2 │        1         0         2  missing    missing    missing
   3 │        1         1         1       -1.0        3.0       -1.4
   4 │        2   missing         3        2.3  missing          3.0
   5 │        2         2        -2       10.0  missing       -100.0

julia> byrow(ds, findlast, :, by = ismissing)
5-element PooledArrays.PooledVector{Union{Missing, Symbol}, UInt32, Vector{UInt32}}:
 :x3_float
 :x3_float
 missing
 :x2_float
 :x2_float

julia> byrow(ds, findlast, 1:3, item = [1,1,1,1,1])
5-element PooledArrays.PooledVector{Union{Missing, Symbol}, UInt32, Vector{UInt32}}:
 :g
 :g
 :x2_int
 missing
 missing

julia> ds = Dataset(x1 = [1,2,2], x2 = [5,6,7], x3 = [8,9,10])
3×3 Dataset
 Row │ x1        x2        x3
     │ identity  identity  identity
     │ Int64?    Int64?    Int64?
─────┼──────────────────────────────
   1 │        1         5         8
   2 │        2         6         9
   3 │        2         7        10

julia> byrow(ds, select, :, with = byrow(ds, findlast, :, by = isodd))
3-element Vector{Union{Missing, Int64}}:
 5
 9
 7
```
@@@@select@@@@
    byrow(ds::AbstractDataset, select, cols; [with, threads])

Select value of `with` among `cols`. The `with` must be a vector of column names(`Symbol` or `String`) or column index (relative to column position in `cols`) or a column name which contains this information.

For heterogeneous column types, `byrow` use `promote_type` for the output. If the column select doesn't exist among `cols`, `byrow` returns `missing`.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(findfirst)`](@ref), [`byrow(findlast)`](@ref)

# Examples
```jldoctest
julia> ds = Dataset(x1 = [1,2,3,4],
            x2 = [1.5,6.5,3.4,2.4],
            x3 = [true, false, true, false],
            y1 = ["x2", "x1", missing, "x2"],
            y2 = [:x2, :x1, missing, :x2],
            y3 = [3,1,1,2])
4×6 Dataset
 Row │ x1        x2        x3        y1        y2        y3
     │ identity  identity  identity  identity  identity  identity
     │ Int64?    Float64?  Bool?     String?   Symbol?   Int64?
─────┼────────────────────────────────────────────────────────────
   1 │        1       1.5      true  x2        x2               3
   2 │        2       6.5     false  x1        x1               1
   3 │        3       3.4      true  missing   missing          1
   4 │        4       2.4     false  x2        x2               2

julia> byrow(ds, select, 1:2, with = :y1)
4-element Vector{Union{Missing, Float64}}:
 1.5
 2.0
  missing
 2.4

julia> byrow(ds, select, [2,1,3], with = :y3)
4-element Vector{Union{Missing, Float64}}:
 1.0
 6.5
 3.4
 4.0

julia> byrow(ds, select, [2,1,3], with = [3,1,1,2])
4-element Vector{Union{Missing, Float64}}:
 1.0
 6.5
 3.4
 4.0

julia> ds = Dataset(x1 = [1,2,2], x2 = [5,6,7], x3 = [8,9,10])
3×3 Dataset
 Row │ x1        x2        x3
     │ identity  identity  identity
     │ Int64?    Int64?    Int64?
─────┼──────────────────────────────
   1 │        1         5         8
   2 │        2         6         9
   3 │        2         7        10

julia> byrow(ds, select, :, with = byrow(ds, findfirst, :, by = isodd))
3-element Vector{Union{Missing, Int64}}:
 1
 9
 7
```
@@@@fill!@@@@
    byrow(ds::AbstractDataset, fill!, cols; [with, by = ismissing, rolling = false, threads])

Fill missing (default behaviour) values in `cols` with values from `with`. User can pass a vector of values or a column name to `with`. `byrow` fills the values in-place, so the type of `cols` and `with` must match. By default, `byrow` fills only missing values in `cols`, but, user can pass any function to `by` which `byrow` fills only the values that returns `true` when `by` is called on them.

When `rolling = true`, `byrow` uses `with` to fill the missing values in the first column among `cols` and replace `with` with the updated values in the first column and uses these values to fill the missing values in the second column among `cols` and replace `with` with the updated values in the second column, and continues this process.

Passing `threads = false` disables multitrheaded computations.

`fill!` is a special `byrow` operations, because it changes the input data set rather than producing a vector.

See [`byrow(fill)`](@ref)

# Examples
```jldoctest
julia> ds = Dataset(x = [1,missing,3], y = [missing,2, 3])
3×2 Dataset
 Row │ x         y
     │ identity  identity
     │ Int64?    Int64?
─────┼────────────────────
   1 │        1   missing
   2 │  missing         2
   3 │        3         3

julia> byrow(ds, fill!, :, with = 1:3)
3×2 Dataset
 Row │ x         y
     │ identity  identity
     │ Int64?    Int64?
─────┼────────────────────
   1 │        1         1
   2 │        2         2
   3 │        3         3

julia> ds = Dataset(x = [1,0,3], y = [0,2,3])
3×2 Dataset
 Row │ x         y
     │ identity  identity
     │ Int64?    Int64?
─────┼────────────────────
   1 │        1         0
   2 │        0         2
   3 │        3         3

julia> byrow(ds, fill!, :, with = 1:3, by = isequal(0))
3×2 Dataset
 Row │ x         y
     │ identity  identity
     │ Int64?    Int64?
─────┼────────────────────
   1 │        1         1
   2 │        2         2
   3 │        3         3

julia> ds = Dataset(x = [2,0,3], y = [0,2,3], z = [5,0,0])
3×3 Dataset
 Row │ x         y         z
     │ identity  identity  identity
     │ Int64?    Int64?    Int64?
─────┼──────────────────────────────
   1 │        2         0         5
   2 │        0         2         0
   3 │        3         3         0

julia> byrow(ds, fill!, :, with = [missing, missing, missing], by = isequal(0), rolling = true)
3×3 Dataset
 Row │ x         y         z
     │ identity  identity  identity
     │ Int64?    Int64?    Int64?
─────┼──────────────────────────────
   1 │        2         2         5
   2 │  missing         2         2
   3 │        3         3         3
```
@@@@fill@@@@
    byrow(ds::AbstractDataset, fill, cols; [with, by = ismissing, rolling = false, threads])

Variant of `byrow(fill!)` which passes a copy of `ds` and leaves `ds` unchanged.

See [`byrow(fill!)`](@ref)
@@@@coalesce@@@@
    byrow(ds::AbstractDataset, coalesce, cols; [threads])

Return the first value in each row of `cols` which is not equal to `missing`, if any. Otherwise return `missing`.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(select)`](@ref), [`byrow(findfirst)`](@ref), [`byrow(findlast)`](@ref)

# Examples
```jldoctest
julia> ds = Dataset(x = [1,missing, missing], y = [missing, missing, 3], z = [5, missing, missing])
3×3 Dataset
 Row │ x         y         z
     │ identity  identity  identity
     │ Int64?    Int64?    Int64?
─────┼──────────────────────────────
   1 │        1   missing         5
   2 │  missing   missing   missing
   3 │  missing         3   missing

julia> byrow(ds, coalesce, :)
3-element Vector{Union{Missing, Int64}}:
 1
  missing
 3

julia> byrow(ds, coalesce, [:z, :y, :x])
3-element Vector{Union{Missing, Int64}}:
 5
  missing
 3
```
@@@@maximum@@@@
    byrow(ds::AbstractDataset, maximum, cols; [by = identity, threads])

Return the largest result of calling function `by` on each values in each row (of selected columns). If `cols` is not specified, `byrow`
returns maximum for all numeric columns in `ds`.

Missing values are removed from the calculation. When all values in a row are missing, it returns `missing`.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(minimum)`](@ref), [`byrow(argmax)`](@ref), [`byrow(argmin)`](@ref)

# Examples

```jldoctest
julia> ds = Dataset(x=[1,2,3], y=[1.3,-2.4,5.5], z=[missing, 4,1])
3×3 Dataset
 Row │ x         y         z
     │ identity  identity  identity
     │ Int64?    Float64?  Int64?
─────┼──────────────────────────────
   1 │        1       1.3   missing
   2 │        2      -2.4         4
   3 │        3       5.5         1

julia> byrow(ds, maximum)
3-element Vector{Union{Missing, Float64}}:
 1.3
 4.0
 5.5

julia> byrow(ds, maximum, [1,3])
3-element Vector{Union{Missing, Int64}}:
 1
 4
 3
```
@@@@minimum@@@@
    byrow(ds::AbstractDataset, minimum, cols; [by = identity, threads])

Return the smallest result of calling function `by` on each values in each row (of selected columns). If `cols` is not specified, `byrow`
returns minimum for all numeric columns in `ds`.

Missing values are removed from the calculation. When all values in a row are missing, it returns `missing`.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(maximum)`](@ref), [`byrow(argmax)`](@ref), [`byrow(argmin)`](@ref)

# Examples

```jldoctest
julia> ds = Dataset(x=[1,2,3], y=[1.3,-2.4,5.5], z=[missing, 4,1])
3×3 Dataset
 Row │ x         y         z
     │ identity  identity  identity
     │ Int64?    Float64?  Int64?
─────┼──────────────────────────────
   1 │        1       1.3   missing
   2 │        2      -2.4         4
   3 │        3       5.5         1

julia> byrow(ds, minimum)
3-element Vector{Union{Missing, Float64}}:
  1.0
 -2.4
  1.0

julia> byrow(ds, minimum, [1,3])
3-element Vector{Union{Missing, Int64}}:
 1
 2
 1
```
@@@@argmax@@@@
    byrow(ds::AbstractDataset, argmax, cols; [by = identity, threads])

Return the column name of the maximum result of calling function `by` on each values in each row (of selected columns). If `cols` is not specified, `byrow`
passes all numeric columns in `ds`.

Missing values are removed from the calculation. When all values in a row are missing, it returns `missing`.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(maximum)`](@ref), [`byrow(minimum)`](@ref), [`byrow(argmin)`](@ref)

# Examples
```jldoctest
julia> ds = Dataset(x=[1,2,missing], y=[1.3,-2.4,missing], z=[missing, 4,missing])
3×3 Dataset
 Row │ x         y          z
     │ identity  identity   identity
     │ Int64?    Float64?   Int64?
─────┼───────────────────────────────
   1 │        1        1.3   missing
   2 │        2       -2.4         4
   3 │  missing  missing     missing

julia> byrow(ds, argmax)
3-element PooledArrays.PooledVector{Union{Missing, Symbol}, UInt32, Vector{UInt32}}:
 :y
 :z
 missing

julia> byrow(ds, argmax, 1:2, by = abs)
3-element PooledArrays.PooledVector{Union{Missing, Symbol}, UInt32, Vector{UInt32}}:
 :y
 :y
 missing
```
@@@@argmin@@@@
    byrow(ds::AbstractDataset, argmin, cols; [by = identity, threads])

Return the column name of the minimum result of calling function `by` on each values in each row (of selected columns). If `cols` is not specified, `byrow`
passes all numeric columns in `ds`.

Missing values are removed from the calculation. When all values in a row are missing, it returns `missing`.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(maximum)`](@ref), [`byrow(minimum)`](@ref), [`byrow(argmax)`](@ref)

# Examples
```jldoctest
julia> ds = Dataset(x=[1,2,missing], y=[1.3,-2.4,missing], z=[missing, 4,missing])
3×3 Dataset
 Row │ x         y          z
     │ identity  identity   identity
     │ Int64?    Float64?   Int64?
─────┼───────────────────────────────
   1 │        1        1.3   missing
   2 │        2       -2.4         4
   3 │  missing  missing     missing

julia> byrow(ds, argmin)
3-element PooledArrays.PooledVector{Union{Missing, Symbol}, UInt32, Vector{UInt32}}:
 :x
 :y
 missing

julia> byrow(ds, argmin, 1:2, by = abs)
3-element PooledArrays.PooledVector{Union{Missing, Symbol}, UInt32, Vector{UInt32}}:
 :x
 :x
 missing
```
@@@@issorted@@@@
    byrow(ds::AbstractDataset, issorted, cols; [rev = false, lt = isless, threads])

Test whether the values in rows (in selected `cols`) are in sorted order. Passing `rev = true` test whether the values in rows are in descending order. By default, the order of values is check by the `isless` function, however, user may pass any function to `lt`. The passed function to `lt` must accept two arguments where `byrow` calls `!lt(x2, x1)` when `rev = false` and `!lt(x1, x2)` when `rev = true` on consecutive column values.

Missing values are larger than any other values. User may pass a customised funtion to `lt` to skip missing values.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(isequal)`](@ref), [`byrow(isless)`](@ref), [`byrow(in)`](@ref)

# Examples
```jldoctest
julia> ds = Dataset(x=[1,2,missing], y=[1.3,-2.4,missing], z=[missing, 4,missing])
3×3 Dataset
 Row │ x         y          z
     │ identity  identity   identity
     │ Int64?    Float64?   Int64?
─────┼───────────────────────────────
   1 │        1        1.3   missing
   2 │        2       -2.4         4
   3 │  missing  missing     missing

julia> byrow(ds, issorted, :)
3-element Vector{Bool}:
 1
 0
 1

julia> byrow(ds, issorted, :, lt = (x,y)->isless(abs(x), abs(y)))
3-element Vector{Bool}:
 1
 1
 1

julia> byrow(ds, issorted, :, lt = isequal) # byrow checks !lt(y, x)
3-element Vector{Bool}:
 1
 1
 0

julia> byrow(ds, issorted, :, lt = !isequal)
3-element Vector{Bool}:
 0
 0
 1
```
@@@@join@@@@
    byrow(ds::AbstractDataset, join, cols; [delim = "", last = "", threads])

For each row and selected columns convert values to string and join them into a single string, inserting the given delimiter (if any) between adjacent strings. If `last` is given, it will be used instead of `delim` between the last two strings. Missing values are converted to empty string and `true` and `false` converted to `1` and `0`, respectively.

Passing `threads = false` disables multitrheaded computations.
@@@@hash@@@@
    byrow(ds::AbstractDataset, hash, cols; [by = identity, threads])

Compute an integer hash code of result of calling `by` on each values in each row of selected `cols`. When `cols` is not specified `byrow` compute hash code for all columns in `ds`.

Passing `threads = false` disables multitrheaded computations.
@@@@nunique@@@@
    byrow(ds::AbstractDataset, nunique, cols; [by = identity, count_missing = true])

Return the number of unique values of the result of calling `by` on each values in each row of selected `cols`. When `cols` is not specified, `byrow` returns the number of unique values for all numeric columns. `missing` are counted as distinct value, and passing `count_missing = false` drop missings from the final count.
@@@@mapreduce@@@@
    byrow(ds::AbstractDataset, mapreduce, cols; op = .+, f = identity, init, kwargs...)

Map `f` on each values in each row of selected `cols` and reduce the result by using `op`. Keyword arguments `op` and `init` must be passed.
@@@@var@@@@
    byrow(ds::AbstractDataset, var, cols; [dof = true, by = identity, threads])

Compute the variance of result of calling `by` on each value in each row of `ds` for selected `cols`. When `cols` is not specified `byrow` computes the variance for all numeric columns. By default, degree of freedom is used for denominator, and passing `dof = false` change it to number of values.

Missing values are droped from calculations, and when all values in a row are `missing` it returns `missing`.

Passing `threads = false` disables multitrheaded computations.

# Examples
```jldoctest
julia> ds = Dataset(x1 = [1.0,missing,missing], x2 = [5,6,missing])
3×2 Dataset
 Row │ x1         x2
     │ identity   identity
     │ Float64?   Int64?
─────┼─────────────────────
   1 │       1.0         5
   2 │ missing           6
   3 │ missing     missing

julia> byrow(ds, var, :)
3-element Vector{Union{Missing, Float64}}:
 8.0
  missing
  missing

julia> byrow(ds, var, :, dof = false)
3-element Vector{Union{Missing, Float64}}:
 4.0
 0.0
  missing
```
@@@@std@@@@
    byrow(ds::AbstractDataset, std, cols; [dof = true, by = identity, threads])

Compute the standard deviation of result of calling `by` on each value in each row of `ds` for selected `cols`. When `cols` is not specified `byrow` computes the standard deviation for all numeric columns. By default, degree of freedom is used for denominator, and passing `dof = false` change it to number of values.

Missing values are droped from calculations, and when all values in a row are `missing` it returns `missing`.

Passing `threads = false` disables multitrheaded computations.

# Examples
```jldoctest
julia> ds = Dataset(x1 = [1.0,missing,missing], x2 = [5,6,missing])
3×2 Dataset
 Row │ x1         x2
     │ identity   identity
     │ Float64?   Int64?
─────┼─────────────────────
   1 │       1.0         5
   2 │ missing           6
   3 │ missing     missing

julia> byrow(ds, std, :)
3-element Vector{Union{Missing, Float64}}:
 2.8284271247461903
  missing
  missing

julia> byrow(ds, std, :, dof = false)
3-element Vector{Union{Missing, Float64}}:
 2.0
 0.0
  missing
```
@@@@cumsum!@@@@
    byrow(ds::AbstractDataset, cumsum!, cols; [missings = :ignore, threads])

Replace each value in `cols` by the result of `cumsum` on each row. When `cols` is not specified `byrow` replace every numeric columns. The type of selected column will be promoted to be able to contain the result of computations. By default missing values are filled with the result of preceding calculations, and passing `missings = :skip` leaves `missing` values untouched.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(cumsum)`](@ref), [`byrow(cumprod!)`](@ref), [`byrow(cummax!)`](@ref), [`byrow(cummin!)`](@ref)

# Examples
```jldoctest
julia> ds = Dataset(x1 = [1.0,missing,2.0], x2 = [5,6,missing])
3×2 Dataset
 Row │ x1         x2
     │ identity   identity
     │ Float64?   Int64?
─────┼─────────────────────
   1 │       1.0         5
   2 │ missing           6
   3 │       2.0   missing

julia> byrow(ds, cumsum!)
3×2 Dataset
 Row │ x1         x2
     │ identity   identity
     │ Float64?   Float64?
─────┼─────────────────────
   1 │       1.0       6.0
   2 │ missing         6.0
   3 │       2.0       2.0

julia> ds = Dataset(x1 = [1.0,missing,2.0], x2 = [5,6,missing])
3×2 Dataset
 Row │ x1         x2
     │ identity   identity
     │ Float64?   Int64?
─────┼─────────────────────
   1 │       1.0         5
   2 │ missing           6
   3 │       2.0   missing

julia> byrow(ds, cumsum!, missings = :skip)
3×2 Dataset
 Row │ x1         x2
     │ identity   identity
     │ Float64?   Float64?
─────┼──────────────────────
   1 │       1.0        6.0
   2 │ missing          6.0
   3 │       2.0  missing
```
@@@@cumsum@@@@
    byrow(ds::AbstractDataset, cumsum, cols; [missings = :ignore, threads])

Variant of `byrow(cumsum!)` which pass a copy of `ds` and leave `ds` untouched.
@@@@cumprod!@@@@

    byrow(ds, cumprod!, cols; [missings = :ignore, threads])

Replace each value in `cols` by the result of `cumprod` on each row. When `cols` is not specified `byrow` replace every numeric columns. The type of selected column will be promoted to be able to contain the result of computations. By default missing values are filled with the result of preceding calculations, and passing `missings = :skip` leaves `missing` values untouched.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(cumprod)`](@ref), [`byrow(cumsum!)`](@ref), [`byrow(cummax!)`](@ref), [`byrow(cummin!)`](@ref)
@@@@cumprod@@@@

    byrow(ds, cumprod!, cols; [missings = :ignore, threads])

Variant of `byrow(cumprod!)` which pass a copy of `ds` and leave `ds` untouched.
@@@@cummax!@@@@
    byrow(ds::AbstractDataset, cummax!, cols; [missings = :ignore, threads])

Replace each value in `cols` by the result of cumulative maximum on each row. When `cols` is not specified `byrow` replace every numeric columns. The type of selected column will be promoted to be able to contain the result of computations. By default missing values are filled with the result of preceding calculations, and passing `missings = :skip` leaves `missing` values untouched.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(cummax)`](@ref), [`byrow(cumsum!)`](@ref), [`byrow(cumprod!)`](@ref), [`byrow(cummin!)`](@ref)
@@@@cummax@@@@
    byrow(ds::AbstractDataset, cummax, cols; [missings = :ignore, threads])

Variant of `byrow(cummax!)` which pass a copy of `ds` and leave `ds` untouched.
@@@@cummin!@@@@
    byrow(ds::AbstractDataset, cummin!, cols; [missings = :ignore, threads])

Replace each value in `cols` by the result of cumulative minimum on each row. When `cols` is not specified `byrow` replaces every numeric columns. The type of selected column will be promoted to be able to contain the result of computations. By default missing values are filled with the result of preceding calculations, and passing `missings = :skip` leaves `missing` values untouched.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(cummin)`](@ref), [`byrow(cumsum!)`](@ref), [`byrow(cumprod!)`](@ref), [`byrow(cummax!)`](@ref)
@@@@cummin@@@@
    byrow(ds::AbstractDataset, cummin, cols; [missings = :ignore, threads])

Variant of `byrow(cummin!)` which pass a copy of `ds` and leave `ds` untouched.
@@@@sort!@@@@
    byrow(ds::AbstractDataset, sort!, cols; [threads, kwargs...])

Update `ds` in place with sorted values in each row of selected `cols`. When `cols` is not specified `byrow` uses every numeric columns. User can pass any keyword argument support by Julia `sort` function. Columns in `cols` will be promoted to be able to contain the new sorted values.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(sort)`](@ref), [`sort!`](@ref)

# Examples
```jldoctest
julia> ds = Dataset(x1 = [1.0,missing,2.0], x2 = [5,6,missing])
3×2 Dataset
 Row │ x1         x2
     │ identity   identity
     │ Float64?   Int64?
─────┼─────────────────────
   1 │       1.0         5
   2 │ missing           6
   3 │       2.0   missing

julia> byrow(ds, sort!, :)
3×2 Dataset
 Row │ x1        x2
     │ identity  identity
     │ Float64?  Float64?
─────┼─────────────────────
   1 │      1.0        5.0
   2 │      6.0  missing
   3 │      2.0  missing

julia> ds = Dataset(x1 = [1.0,missing,2.0], x2 = [5,6,missing])
3×2 Dataset
 Row │ x1         x2
     │ identity   identity
     │ Float64?   Int64?
─────┼─────────────────────
   1 │       1.0         5
   2 │ missing           6
   3 │       2.0   missing

julia> byrow(ds, sort!, :, rev = true)
3×2 Dataset
 Row │ x1         x2
     │ identity   identity
     │ Float64?   Float64?
─────┼─────────────────────
   1 │       5.0       1.0
   2 │ missing         6.0
   3 │ missing         2.0
```
@@@@sort@@@@
    byrow(ds::AbstractDataset, sort, cols; [threads, kwargs...])

Variant of `byrow(sort!)` which pass a copy of `ds` and leave `ds` untouched.
@@@@stdze!@@@@
    byrow(ds::AbstractDataset, stdze!, cols; [threads])

Replace each value in each row of `ds` for selected `cols` by its standardised values.

Passing `threads = false` disables multitrheaded computations.

See [`byrow(stdze)`](@ref)
@@@@stdze@@@@
    byrow(ds::AbstractDataset, stdze, cols; [threads])

Variant of `byrow(stdze!)` which pass a copy of `ds` and leave `ds` untouched.
@@@@generic@@@@
    byrow(ds::AbstractDataset, fun, cols; [threads])

Return the result of calling `fun` on each row of `ds` selected by `cols`. The `fun` function must accept one argument which contains the values of each row as a vector of values and return a scalar.

When user passes a type as `fun` and a single column as `cols`,  `byrow` convert the corresponding column to the type specified by `fun`.

For generic functions there are two special cases:

* When `cols` is a single column, `byrow(ds, fun, cols)` acts like `fun.(ds[:, cols])`
* When `cols` is referring to exactly two columns and it is possible to pass two vectors as arguments of `fun`, `byrow` returns `fun.(ds[:, col1], ds[:, col2])` when possible.
"""


text:
include("hw1/prob4.jl")

using .hw1_4

const μ = 3.986004e5 #km^3/s^2
r1 = 8000. #km
r2 = 10r1

# println(hw1_4.oe_to_xyz([8000, 0, 0, 0, 0, 0], μ))
hw1_4.plot_maneuver(r1, r2, π/2, μ, π/2)
# hw1_4.plot_maneuver(r1, r2, π/2, μ, 0.1)


text:
using ArgParse, Random, SparseArrays, HDF5, H5Sparse, SAG

const METADATA_BYTES = 8
const ELEMENT_TYPE = Float64
const CANARY_VALUE = UInt16(2^16 - 1)

include("gradients.jl")

function update_argsettings!(s::ArgParseSettings)
    @add_arg_table s begin
        "--lambda"
            help = "L2 regularizer coefficient"
            arg_type = Float64
            default = 0.0
        "--labeldataset"
            help = "Label dataset name"
            default = "b"
            arg_type = String
        "--iteratedataset"
            help = "Initial iterate dataset name (chosen to be all-zeros if not provided)"
            arg_type = String
        "--nsubpartitions"
            help = "Number of sub-partitions to split the data stored at each worker into"
            arg_type = Int
            default = 1
            range_tester = (x) -> x >= 1
        "--stepsize"
            help = "Gradient descent step size"
            arg_type = Float64
            default = 1.0
            range_tester = (x) -> x > 0
        "--variancereduced"
            help = "Compute a variance-reduced gradient in each iteration"
            action = :store_true
        "--vralgo"
            help = "Specifies which variance-reduced optimizer to use (tree or table)"
            arg_type = String
            default = "tree"
            range_tester = (x) -> x in ["tree", "table"]
        "--nostale"
            help = "If set, do not store stale gradients (to conform with SAG)"
            action = :store_true
        "--nreplicas"
            help = "Number of replicas of each data partition"
            default = 1
            arg_type = Int
            range_tester = (x) -> x >= 1
        "--nwait"
            help = "Number of replicas to wait for in each iteration (defaults to all replicas)"
            arg_type = Int
            range_tester = (x) -> x >= 1
        "--nwaitschedule"
            help = "Factor by which nwait is reduced per iteration"
            arg_type = Float64
            default = 1.0
            range_tester = (x) -> 0 < x <= 1
        "--kickstart"
            help = "Wait for all partitions in the first iteration"
            action = :store_true
    end
end

function update_parsed_args!(s::ArgParseSettings, parsed_args)
    parsed_args[:algorithm] = "logreg.jl"
    nworkers::Int = parsed_args[:nworkers]
    nreplicas::Int = parsed_args[:nreplicas]
    mod(nworkers, nreplicas) == 0 || throw(ArgumentError("nworkers is $nworkers, but must be divisible by nreplicas"))
    npartitions = div(nworkers, nreplicas)
    parsed_args[:nwait] = isnothing(parsed_args[:nwait]) ? npartitions : parsed_args[:nwait]
    0 < parsed_args[:nwait] <= npartitions || throw(ArgumentError("nwait is $(parsed_args[:nwait]), but nworkers/nreplicas is $npartitions"))

    # record the number of rows and columns of the dataset
    nrows, ncolumns = problem_size(parsed_args[:inputfile], parsed_args[:inputdataset])
    parsed_args[:nrows] = nrows
    parsed_args[:ncolumns] = ncolumns

    parsed_args
end

"""

Called inside `asyncmap!` to determine if enough workers have responded. Returns `true` if at
least `nwait` workers have responded and `false` otherwise.
"""
function fwait(epoch, repochs; nworkers, nwait, kickstart, nwaitschedule, kwargs...)
    length(repochs) == nworkers || throw(DomainError(nworkers, "repochs must have length nworkers"))
    0 < nwait <= nworkers || throw(ArgumentError("nwait is $nwait, but must be in [1, npartitions]"))
    nwait = min(max(1, ceil(Int, nwait*nwaitschedule^epoch)), nworkers)
    nrec = 0
    for repoch in repochs
        nrec += repoch == epoch
    end
    if epoch == 1 && kickstart
        return nrec == nwait
    else
        return nrec >= nwait
    end
end

"""

Return the size of the data matrix.
"""
function problem_size(filename::String, dataset::String)
    HDF5.ishdf5(filename) || throw(ArgumentError("$filename isn't an HDF5 file"))
    h5open(filename, "r") do fid
        dataset in keys(fid) || throw(ArgumentError("$dataset is not in $fid"))
        if H5Sparse.h5isvalidcsc(fid, dataset)
            return size(H5SparseMatrixCSC(fid, dataset))
        else
            return size(fid[dataset])
        end
    end
end

function partition_samples(X::H5SparseMatrixCSC, nsubpartitions::Integer)
    nsamples = size(X, 2)
    dividers = round.(Int, range(1, nsamples+1, length=nsubpartitions+1))
    [sparse(X[:, dividers[i]:(dividers[i+1]-1)]) for i in 1:nsubpartitions]
end

function partition_samples(X::Matrix, nsubpartitions::Integer)
    nsamples = size(X, 2)
    dividers = round.(Int, range(1, nsamples+1, length=nsubpartitions+1))
    [view(X, :, dividers[i]:(dividers[i+1]-1)) for i in 1:nsubpartitions]
end

function partition_samples(v::Vector, nsubpartitions::Integer)
    nsamples = length(v)
    dividers = round.(Int, range(1, nsamples+1, length=nsubpartitions+1))
    [view(v, dividers[i]:(dividers[i+1]-1)) for i in 1:nsubpartitions]
end

function read_localdata(i::Integer, nworkers::Integer; inputfile::String, inputdataset::String, labeldataset::String, nreplicas::Integer, kwargs...)
    HDF5.ishdf5(inputfile) || throw(ArgumentError("$inputfile isn't an HDF5 file"))
    0 < nworkers || throw(DomainError(nworkers, "nworkers must be positive"))
    0 < nreplicas || throw(DomainError(nreplicas, "nreplicas must be positive"))
    0 < i <= nworkers || throw(DomainError(i, "i must be in [1, nworkers]"))
    mod(nworkers, nreplicas) == 0 || throw(ArgumentError("nworkers must be divisible by nreplicas"))
    npartitions = div(nworkers, nreplicas)
    partition_index = ceil(Int, i/nreplicas)
    _, nsamples = problem_size(inputfile, inputdataset)
    h5open(inputfile, "r") do fid
        inputdataset in keys(fid) || throw(ArgumentError("$inputdataset is not in $fid"))
        labeldataset in keys(fid) || throw(ArgumentError("$labeldataset is not in $fid"))

        # determine the indices of the samples to be stored locally
        Is = partition(nsamples, npartitions, partition_index)

        # load those samples into memory
        labels = fid[labeldataset][Is]
        if H5Sparse.h5isvalidcsc(fid, inputdataset)
            return sparse(H5SparseMatrixCSC(fid, inputdataset, :, Is))::SparseMatrixCSC, labels
        else
            return fid[inputdataset][:, Is]::Matrix, labels
        end
    end
end

function worker_setup(rank::Integer, nworkers::Integer; kwargs...)
    0 < nworkers || throw(DomainError(nworkers, "nworkers must be positive"))
    features, labels = read_localdata(rank, nworkers; kwargs...)
    dimension = size(features, 1)
    localdata = (features, labels)
    to_worker_metadata_bytes = sizeof(UInt16) * 2 * nworkers
    recvbuf = Vector{UInt8}(undef, sizeof(ELEMENT_TYPE)*(dimension+1) + to_worker_metadata_bytes)
    sendbuf = Vector{UInt8}(undef, sizeof(ELEMENT_TYPE)*(dimension+1) + COMMON_BYTES + METADATA_BYTES)
    localdata, recvbuf, sendbuf
end

function coordinator_setup(nworkers::Integer; inputfile::String, inputdataset::String, iteratedataset, kwargs...)
    0 < nworkers || throw(DomainError(nworkers, "nworkers must be positive"))
    dimension, nsamples = problem_size(inputfile, inputdataset)

    # initial iterate
    V = zeros(Float64, dimension+1)
    if !isnothing(iteratedataset) # given as an argument and loaded from disk
        h5open(inputfile) do fid
            iteratedataset in keys(fid) || throw(ArgumentError("iterate dataset $iteratedataset not found"))
            size(fid[iteratedataset]) == (dimension+1,) || throw(DimensionMismatch("Expected iterate to have dimensions $((dimension+1,)), but it has dimensions $(size(fid[iteratedataset]))"))
            V .= fid[iteratedataset][:]
        end
    end

    # communication buffers
    to_worker_metadata_bytes = sizeof(UInt16) * 2 * nworkers
    sendbuf = Vector{UInt8}(undef, sizeof(ELEMENT_TYPE)*(dimension+1) + to_worker_metadata_bytes)
    recvbuf = Vector{UInt8}(undef, (sizeof(ELEMENT_TYPE)*(dimension+1) + COMMON_BYTES + METADATA_BYTES) * nworkers)
    reinterpret(ELEMENT_TYPE, view(sendbuf, (to_worker_metadata_bytes+1):length(sendbuf))) .= view(V, :)

    # setup state for the gradient tracking method
    gradient_state, gradient_recvf = setup_gradient(V; kwargs...)
    iterate_state = nothing
    state = (gradient_state, iterate_state)

    V, recvbuf, sendbuf, state, gradient_recvf
end

metadata_view(buffer) = view(buffer, COMMON_BYTES+1:(COMMON_BYTES + METADATA_BYTES))
data_view(buffer) = reinterpret(ELEMENT_TYPE, @view buffer[(COMMON_BYTES + METADATA_BYTES+1):end])

function worker_task!(recvbuf, sendbuf, localdata; state=nothing, ncolumns::Integer, nworkers::Integer, kwargs...)
    to_worker_metadata_bytes = sizeof(UInt16) * 2 * nworkers
    features, labels = localdata
    size(features, 2) == length(labels) || throw(DimensionMismatch("features has dimensions $(size(features)), but labels has dimension $(size(labels))"))
    dimension, nlocalsamples = size(features)

    # prepare working memory
    if isnothing(state) # first iteration
        nsubpartitions_prev = 0
        subpartition_index = 1
        w = zeros(Float64, nlocalsamples) # temp. storage
    else # subsequent iterations
        nsubpartitions_prev::Int, subpartition_index::Int, w::Vector{Float64} = state
    end

    # get the number of sub-partitions from the coordinator
    vs = reinterpret(Tuple{UInt16,UInt16}, view(recvbuf, 1:to_worker_metadata_bytes))
    rank <= length(vs) || throw(DimensionMismatch("vs has length $(length(vs)), but rank is $rank"))
    nsubpartitions, _ = vs[rank]

    # increment the sub-partition index
    if !iszero(nsubpartitions_prev)
        subpartition_index = mod(subpartition_index, nsubpartitions_prev) + 1
        subpartition_index = align_partitions(nlocalsamples, nsubpartitions_prev, nsubpartitions, subpartition_index)
    end
    nsubpartitions_prev = nsubpartitions
    0 < nsubpartitions <= nlocalsamples || throw(DimensionMismatch("nsubpartitions is $nsubpartitions, but nlocalsamples is $nlocalsamples"))
    0 < subpartition_index <= nsubpartitions || throw(ArgumentError("subpartition_index is $subpartition_index, but nsubpartitions is $nsubpartitions"))

    # format the recvbuf into a matrix we can operate on
    recvdata = view(recvbuf, (to_worker_metadata_bytes+1):length(recvbuf))
    length(reinterpret(ELEMENT_TYPE, recvdata)) == dimension+1 || throw(DimensionMismatch("recvdata has length $(length(reinterpret(ELEMENT_TYPE, recvdata))), but the data dimension is $dimension"))
    v = reinterpret(ELEMENT_TYPE, recvdata)

    # indices of the local samples to process in this iteration
    cols = partition(nlocalsamples, nsubpartitions, subpartition_index)

    # gradient computation
    tcolsmul!(w, features, view(v, 2:length(v)), cols)
    @views w[cols] .+= v[1] # implicit intercept
    @views w[cols] .*= labels[cols]
    @views w[cols] .= exp.(w[cols])
    @views w[cols] .+= 1
    @views w[cols] .= labels[cols] ./ w[cols]
    @views w[cols] .*= -1
    @views v[1] = sum(w[cols]) # intercept derivative
    colsmul!(view(v, 2:length(v)), features, w, cols)
    v ./= ncolumns # normalize by the total number of samples

    # populate the send buffer
    metadata = reinterpret(UInt16, metadata_view(sendbuf))
    metadata[1] = CANARY_VALUE
    metadata[2] = rank
    metadata[3] = nsubpartitions
    metadata[4] = subpartition_index
    data_view(sendbuf) .= v
    nsubpartitions_prev, subpartition_index, w
end

function setup_gradient(args...; variancereduced::Bool, vralgo::String, kwargs...)
    if variancereduced
        if vralgo == "table"
            return nothing, nothing
        elseif vralgo == "tree"
            return setup_gradient_vrt(args...; kwargs...), recvf_vrt
        else
            throw(ArgumentError("unexpected vralgo $vralgo"))
        end
    else
        return nothing, nothing
    end
end

function update_gradient!(args...; variancereduced::Bool, vralgo::String, kwargs...)
    if variancereduced
        if vralgo == "table"
            return update_gradient_vr!(args...; kwargs...)
        elseif vralgo == "tree"
            return update_gradient_vrt!(args...; kwargs...)
        else
            throw(ArgumentError("unexpected vralgo $vralgo"))
        end
    else
        return update_gradient_sgd!(args...; kwargs...)
    end
end

function update_iterate!(v, ∇; state=nothing, stepsize, lambda, kwargs...)
    size(v) == size(∇) || throw(DimensionMismatch("v has dimensions $(size(v)), but ∇ has dimensions $(size(∇))"))
    v .*= 1 - stepsize * lambda
    v .-= stepsize .* ∇
    return
end

function coordinator_task!(V, ∇, recvbufs, sendbuf, epoch::Integer, repochs::Vector{<:Integer}; nworkers::Integer, state=nothing, kwargs...)
    isnothing(state) || length(state) == 2 || throw(ArgumentError("expected state to be nothing or a tuple of length 2, but got $state"))
    if isnothing(state)
        gradient_state = update_gradient!(∇, recvbufs, epoch, repochs; kwargs...)
        iterate_state = update_iterate!(V, ∇; kwargs...)
    else
        gradient_state, iterate_state = state
        gradient_state = update_gradient!(∇, recvbufs, epoch, repochs; state=gradient_state, kwargs...)
        iterate_state = update_iterate!(V, ∇; state=iterate_state, kwargs...)
    end
    to_worker_metadata_bytes = sizeof(UInt16) * 2 * nworkers
    reinterpret(ELEMENT_TYPE, view(sendbuf, (to_worker_metadata_bytes+1):length(sendbuf))) .= view(V, :)
    gradient_state, iterate_state
end

include("common.jl")


text:
using Pkg; for p in ("Knet",); haskey(Pkg.installed(),p) || Pkg.add(p); end

using Knet

struct CharLM; input; rnn; output; end

CharLM(vocab::Int,input::Int,hidden::Int; o...) =
    CharLM(Embed(vocab,input), RNN(input,hidden; o...), Linear(hidden,vocab))

function (c::CharLM)(x; pdrop=0, hidden=nothing)
    x = c.input(x)                # (B,T)->(X,B,T)
    x = dropout(x, pdrop)
    x = c.rnn(x, hidden=hidden)   # (H,B,T)
    x = dropout(x, pdrop)
    x = reshape(x, size(x,1), :)  # (H,B*T)
    return c.output(x)            # (V,B*T)
end

struct Embed; w; end

Embed(vocab::Int,embed::Int)=Embed(param(embed,vocab))

(e::Embed)(x) = e.w[:,x]

struct Linear; w; b; end

Linear(input::Int, output::Int)=Linear(param(output,input), param0(output))

(l::Linear)(x) = l.w * x .+ l.b

# Sample from trained model

function generate(model,chars,n)
    function sample(y)
        p = Array(exp.(y)); r = rand()*sum(p)
        for j=1:length(p); (r -= p[j]) < 0 && return j; end
    end
    x = 1
    h = []
    for i=1:n
        y = model([x], hidden=h)
        x = sample(y)
        print(chars[x])
    end
    println()
end

@info("Loading Shakespeare data")
include(Knet.dir("data","gutenberg.jl"))
trn,tst,shake_chars1 = shakespeare()
shake_text = String(shake_chars1[vcat(trn,tst)])

@info("Loading Shakespeare model")
isfile("shakespeare.jld2") || download("http://people.csail.mit.edu/deniz/models/tutorial/shakespeare.jld2","shakespeare.jld2")
shake_model, shake_chars = Knet.load("shakespeare.jld2","model","chars")

@info("Reading Julia files")
base = joinpath(Sys.BINDIR, Base.DATAROOTDIR, "julia")
julia_text = ""
for (root,dirs,files) in walkdir(base)
    for f in files
        global julia_text
        f[end-2:end] == ".jl" || continue
        julia_text *= read(joinpath(root,f),String)
    end
    # println((root,length(files),all(f->contains(f,".jl"),files)))
end

@info("Loading Julia model")
isfile("juliacharlm.jld2") || download("http://people.csail.mit.edu/deniz/models/tutorial/juliacharlm.jld2","juliacharlm.jld2")
julia_model, julia_chars = Knet.load("juliacharlm.jld2","model","chars")

nothing

max_stars_repo_path: Courses/World of machine learning/charlm/charlm.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
function lower_varname(var::Variable, idv, order)
    order == 0 && return var
    name = Symbol(var.name, :ˍ, string(idv.name)^order)
    return Variable{vartype(var)}(name)
end

function flatten_differential(O::Operation)
    @assert is_derivative(O) "invalid differential: $O"
    is_derivative(O.args[1]) || return (O.args[1], O.op.x, 1)
    (x, t, order) = flatten_differential(O.args[1])
    isequal(t, O.op.x) || throw(ArgumentError("non-matching differentials on lhs: $t, $(O.op.x)"))
    return (x, t, order + 1)
end

"""
$(TYPEDSIGNATURES)

Takes a Nth order ODESystem and returns a new ODESystem written in first order
form by defining new variables which represent the N-1 derivatives.
"""
function ode_order_lowering(sys::ODESystem)
    eqs_lowered, new_vars = ode_order_lowering(sys.eqs, sys.iv)
    ODESystem(eqs_lowered, sys.iv, vcat(new_vars, states(sys)), sys.ps)
end

function ode_order_lowering(eqs, iv)
    var_order = Dict{Variable,Int}()
    vars = Variable[]
    new_eqs = Equation[]
    new_vars = Variable[]

    for (i, eq) ∈ enumerate(eqs)
        if isequal(eq.lhs, Constant(0))
            push!(new_eqs, eq)
        else
            var, maxorder = var_from_nested_derivative(eq.lhs)
            if maxorder > get(var_order, var, 0)
                var_order[var] = maxorder
                any(isequal(var), vars) || push!(vars, var)
            end
            var′ = lower_varname(var, iv, maxorder - 1)
            rhs′ = rename_lower_order(eq.rhs)
            push!(new_eqs,Differential(iv())(var′(iv())) ~ rhs′)
        end
    end

    for var ∈ vars
        order = var_order[var]
        for o in (order-1):-1:1
            lvar = lower_varname(var, iv, o-1)
            rvar = lower_varname(var, iv, o)
            push!(new_vars, rvar)

            rhs = rvar(iv())
            eq = Differential(iv())(lvar(iv())) ~ rhs
            push!(new_eqs, eq)
        end
    end

    return (new_eqs, new_vars)
end

function rename_lower_order(O::Expression)
    isa(O, Operation) || return O
    if is_derivative(O)
        (x, t, order) = flatten_differential(O)
        return lower_varname(x.op, t.op, order)(x.args...)
    end
    return Operation(O.op, rename_lower_order.(O.args))
end