text:
Britter_McQuaid_correlations = Dict(
    0.010 => ( αs=[-1.0, -0.7, -0.29, -0.2, 1.0],
               βs=[2.25, 2.25, 2.45, 2.45, 1.83]),
    0.005 => ( αs=[-1.0, -0.67, -0.28, -0.15, 1.0],
               βs=[2.4, 2.4, 2.63, 2.63, 2.07]),
    0.020 => ( αs=[-1.0, -0.69, -0.31, -0.16, 1.0],
               βs=[2.08, 2.08, 2.25, 2.25, 1.62]),
    0.002 => ( αs=[-1.0, -0.69, -0.25, -0.13, 1.0],
               βs=[2.6, 2.6, 2.77, 2.77, 2.21]),
    0.100 => ( αs=[-1.0, -0.55, -0.14, 1.0],
               βs=[1.75, 1.75, 1.85, 1.28]),
    0.050 => ( αs=[-1.0, -0.68, -0.29, -0.18, 1.0],
               βs=[1.92, 1.92, 2.06, 2.06, 1.4]),
)

"""
    britter_plume_factory(scenario)

Generates a Britter-McQuaid dispersion model on the given scenario and returns a
function giving the centerline concentration of the form
c(x, [y, z, t])

Currently only implements the max concentration at a downwind distance x, the
other coordinates are ignored.

"""
function britter_plume_factory(scenario)

    required_params = [:mass_emission_rate, :release_height, :jet_density,
                       :release_temperature, :windspeed, :ambient_density,
                       :ambient_temperature, :pasquill_gifford]
    if all(key -> !(ismissing(getproperty(scenario,key))), required_params)
        Q = scenario.mass_emission_rate
        h = scenario.release_height
        ρⱼ = scenario.jet_density
        Tᵣ = scenario.release_temperature

        u = scenario.windspeed
        ρₐ = scenario.ambient_density
        Tₐ = scenario.ambient_temperature
        class = scenario.pasquill_gifford
    else
        missing_params = [ String(i) for i in filter(key -> ismissing(getproperty(scenario,key)), required_params)]
        error_string = "These parameters cannot be missing: " * join(missing_params, ", ")
        e = MissingException(error_string)
        throw(e)
    end

    # Setting up the Britter-McQuaid curves
    britter_interps = [ ]
    concs = sort(collect(keys(Britter_McQuaid_correlations)), rev=true)

    for conc in concs
        αs, βs = Britter_McQuaid_correlations[conc]
        f = LinearInterpolation(αs, βs, extrapolation_bc=Line())
        push!(britter_interps, (c=conc, it=f))
    end

    # determining the windspeed at 10m
    wind = windspeed(u, h, class)
    u₁₀ = wind(10)

    # relative density
    g = 9.80616  # gravity, m/s^2
    gₒ = g * ((ρⱼ - ρₐ)/ ρₐ)

    # critical distance
    Vr = Q/ρⱼ # volumetric release rate
    D = √(Vr/u₁₀)

    # correlation parameter
    α = 0.2*log10( gₒ^2 * Vr * u₁₀^-5 )

    # setting up correlation for constant α
    # calculates the points for the linear interpolation
    concs = [ elem.c for elem in britter_interps ]
    βs = [ elem.it(α) for elem in britter_interps ]

    # linear interpolation between the short distance correlation
    # and the main correlation
    βs = [ log10(30); βs]
    concs = [ 306*30^-2 / (1+ 306*30^-2); concs ]

    # linear interpolation, extrapolates past the end with a straight line
    it_new = LinearInterpolation(βs, concs, extrapolation_bc=Line())

    function britter_plume(x, y=0, z=0)
        T′ = Tₐ/Tᵣ
        x′ = x/D
        if x′ < 30
            # use short distance correlation
            c′ = x′ > 0 ? 306*x′^-2 / (1+ 306*x′^-2) : 1.0
        else
            # use linear interpolation
            β = log10(x′)
            c′ = it_new(β)
        end
        c = ( ρⱼ*c′*T′) / (1 - c′ + c′*T′)
        return c
    end

    return britter_plume
end


text:
#!/usr/bin/env julia

using Distributed

push!(LOAD_PATH, pwd() )

global iterations = 100000
global top_length = 32768

# Create a random chromosome
function random_chromosome(length)

    chromosome = Bool[]

    for i in 1:1:length
        if ( rand(0:1) > 0.5 )
      	    push!(chromosome,true)
        else
	    push!(chromosome,false)
        end
    end
    chromosome
end

# Mutate and compute time
function time_mutations(number,len)
  inicioTiempo = time()

  for i in 0:1:number
      indi = random_chromosome(len)
      rr = @distributed (+) for j in [indi[n:n+3] for n=1:3:length(indi)-3]
          value = reduce(+, j )
          if value == 0 || value == 4
              1
          else
              0
          end
      end
  end
  time()-inicioTiempo

end

function run_royalroad()
    len = 16
    while len <= top_length
        println("julia-RoyalRoad-d4, " * string(len) * ", " * string(time_mutations( iterations,len )))
        len *= 2
    end
end

run_royalroad()

ext: jl


text:
# TODO:
# * add handling of empty ByRow to filter, and select/transform/combine for GroupedDataFrame
# * add handling of multiple column return rules for select/transform/combine for GroupedDataFrame

# normalize_selection function makes sure that whatever input format of idx is it
# will end up in one of four canonical forms
# 1) AbstractVector{Int}
# 2) Pair{Int, <:Pair{<:Base.Callable, <:Union{Symbol, Vector{Symbol}, Type{AsTable}}}}
# 3) Pair{AbstractVector{Int}, <:Pair{<:Base.Callable, <:Union{Symbol, AbstractVector{Symbol}, Type{AsTable}}}}
# 4) Pair{AsTable, <:Pair{<:Base.Callable, <:Union{Symbol, Vector{Symbol}, Type{AsTable}}}}
# 5) Callable

"""
    ByRow

A type used for selection operations to signal that the wrapped function should
be applied to each element (row) of the selection.

Note that `ByRow` always collects values returned by `fun` in a vector.
"""
struct ByRow{T} <: Function
    fun::T
end

(f::ByRow)(cols::AbstractVector...) = f.fun.(cols...)
(f::ByRow)(table::NamedTuple) = f.fun.(Tables.namedtupleiterator(table))

# add a method to funname defined in other/utils.jl
funname(row::ByRow) = funname(row.fun)

normalize_selection(idx::AbstractIndex, sel, renamecols::Bool) =
    try
        idx[sel]
    catch e
        if e isa MethodError && e.f === getindex && e.args === (idx, sel)
            throw(ArgumentError("Unrecognized column selector: $sel"))
        else
            rethrow(e)
        end
    end

normalize_selection(idx::AbstractIndex, sel::Base.Callable, renamecols::Bool) = sel
normalize_selection(idx::AbstractIndex, sel::Colon, renamecols::Bool) = idx[:]

normalize_selection(idx::AbstractIndex, sel::Pair{typeof(nrow), Symbol},
                    renamecols::Bool) =
    length(idx) == 0 ? (Int[] => (() -> 0) => last(sel)) : (1 => length => last(sel))
normalize_selection(idx::AbstractIndex, sel::Pair{typeof(nrow), <:AbstractString},
                    renamecols::Bool) =
    normalize_selection(idx, first(sel) => Symbol(last(sel)), renamecols)
normalize_selection(idx::AbstractIndex, sel::typeof(nrow), renamecols::Bool) =
    normalize_selection(idx, nrow => :nrow, renamecols)

function normalize_selection(idx::AbstractIndex, sel::ColumnIndex, renamecols::Bool)
    c = idx[sel]
    return c => identity => _names(idx)[c]
end

function normalize_selection(idx::AbstractIndex, sel::Pair{<:ColumnIndex, Symbol},
                             renamecols::Bool)
    c = idx[first(sel)]
    return c => identity => last(sel)
end

normalize_selection(idx::AbstractIndex, sel::Pair{<:ColumnIndex, <:AbstractString},
                    renamecols::Bool) =
    normalize_selection(idx, first(sel) => Symbol(last(sel)), renamecols::Bool)

function normalize_selection(idx::AbstractIndex,
                             sel::Pair{<:ColumnIndex,
                                       <:Pair{<:Base.Callable,
                                              <:Union{Symbol, AbstractString}}},
                             renamecols::Bool)
    src, (fun, dst) = sel
    return idx[src] => fun => Symbol(dst)
end

function normalize_selection(idx::AbstractIndex,
                             sel::Pair{<:Any,
                                       <:Pair{<:Base.Callable,
                                              <:Union{Symbol, AbstractString, DataType,
                                                      AbstractVector{Symbol},
                                                      AbstractVector{<:AbstractString}}}},
                             renamecols::Bool)
    lls = last(last(sel))
    if lls isa DataType
        lls === AsTable || throw(ArgumentError("Only DataType supported as target is AsTable"))
    end
    if first(sel) isa AsTable
        rawc = first(sel).cols
        wanttable = true
    else
        rawc = first(sel)
        wanttable = false
    end
    if rawc isa AbstractVector{Int}
        c = rawc
    elseif rawc isa Union{AbstractVector{Symbol}, AbstractVector{<:AbstractString}}
        c = [idx[n] for n in rawc]
    else
        c = try
                idx[rawc]
            catch e
                if e isa MethodError && e.f === getindex && e.args === (idx, rawc)
                    throw(ArgumentError("Unrecognized column selector: $rawc"))
                else
                    rethrow(e)
                end
            end
    end
    if lls isa AbstractString
        r = Symbol(lls)
    elseif lls isa AbstractVector{<:AbstractString}
        r = Symbol.(lls)
    else
        r = lls
    end
    if r isa AbstractVector{Symbol}
        allunique(r) || throw(ArgumentError("target column names must be unique"))
    end
    return (wanttable ? AsTable(c) : c) => first(last(sel)) => r
end

function normalize_selection(idx::AbstractIndex,
                             sel::Pair{<:ColumnIndex,<:Base.Callable}, renamecols::Bool)
    c = idx[first(sel)]
    fun = last(sel)
    if renamecols
        newcol = Symbol(_names(idx)[c], "_", funname(fun))
    else
        newcol = _names(idx)[c]
    end
    return c => fun => newcol
end

function normalize_selection(idx::AbstractIndex,
                             sel::Pair{<:Any, <:Base.Callable}, renamecols::Bool)
    if first(sel) isa AsTable
        rawc = first(sel).cols
        wanttable = true
    else
        rawc = first(sel)
        wanttable = false
    end
    if rawc isa AbstractVector{Int}
        c = rawc
    elseif rawc isa Union{AbstractVector{Symbol}, AbstractVector{<:AbstractString}}
        c = [idx[n] for n in rawc]
    else
        c = try
                idx[rawc]
            catch e
                if e isa MethodError && e.f === getindex && e.args === (idx, rawc)
                    throw(ArgumentError("Unrecognized column selector: $rawc"))
                else
                    rethrow(e)
                end
            end
    end
    fun = last(sel)
    if length(c) > 3
        prefix = join(@views(_names(idx)[c[1:2]]), '_')
        if renamecols
            newcol = Symbol(prefix, "_etc_", funname(fun))
        else
            newcol = Symbol(prefix, "_etc")
        end
    elseif isempty(c)
        renamecols || throw(ArgumentError("when renamecols=false target column name " *
                                          "must be passed if there are no input columns"))
        newcol = Symbol(funname(fun))
    else
        prefix = join(view(_names(idx), c), '_')
        if renamecols
            newcol = Symbol(prefix, '_', funname(fun))
        else
            newcol = Symbol(prefix)
        end
    end
    return (wanttable ? AsTable(c) : c) => fun => newcol
end

_transformation_helper(df::AbstractDataFrame, col_idx::Nothing, fun) = fun(df)
_transformation_helper(df::AbstractDataFrame, col_idx::Int, fun) = fun(df[!, col_idx])

_empty_astable_helper(fun, len) = [fun(NamedTuple()) for _ in 1:len]

function _transformation_helper(df::AbstractDataFrame, col_idx::AsTable, fun)
    tbl = Tables.columntable(select(df, col_idx.cols, copycols=false))
    if isempty(tbl) && fun isa ByRow
        return _empty_astable_helper(fun.fun, nrow(df))
    else
        return fun(tbl)
    end
end

_empty_selector_helper(fun, len) = [fun() for _ in 1:len]

function _transformation_helper(df::AbstractDataFrame, col_idx::AbstractVector{Int}, fun)
    if isempty(col_idx) && fun isa ByRow
        return _empty_selector_helper(fun.fun, nrow(df))
    else
        cdf = eachcol(df)
        return fun(map(c -> cdf[c], col_idx)...)
    end
end

function _gen_colnames(@nospecialize(res), newname::Union{AbstractVector{Symbol},
                                                          Type{AsTable}, Nothing})
    if res isa AbstractMatrix
        colnames = gennames(size(res, 2))
    else
        colnames = propertynames(res)
    end

    if newname !== AsTable && newname !== nothing
        if length(colnames) != length(newname)
            throw(ArgumentError("Number of returned columns does not match the " *
                                "length of requested output"))
        end
        colnames = newname
    end

    # fix the type to avoid unnecessary compilations of methods
    # this should be cheap
    return colnames isa Vector{Symbol} ? colnames : collect(Symbol, colnames)
end

_expand_to_table(res) = Tables.columntable(res)
_expand_to_table(res::Union{AbstractDataFrame, NamedTuple, DataFrameRow, AbstractMatrix}) = res

function _expand_to_table(res::AbstractVector)
    isempty(res) && return Tables.columntable(res)
    kp1 = keys(res[1])
    prepend = all(x -> x isa Integer, kp1)
    if !(prepend || all(x -> x isa Symbol, kp1) || all(x -> x isa AbstractString, kp1))
        throw(ArgumentError("keys of the returned elements must be " *
                            "`Symbol`s, strings or integers"))
    end
    if any(x -> !isequal(keys(x), kp1), res)
        throw(ArgumentError("keys of the returned elements must be identical"))
    end
    newres = DataFrame()
    for n in kp1
        newres[!, prepend ? Symbol("x", n) : Symbol(n)] = [x[n] for x in res]
    end
    return newres
end

function _insert_row_multicolumn(newdf::DataFrame, df::AbstractDataFrame,
                                 allow_resizing_newdf::Ref{Bool}, colnames::AbstractVector{Symbol},
                                 res::Union{NamedTuple, DataFrameRow})
    if ncol(newdf) == 0
        # if allow_resizing_newdf[] is false we know this is select or transform
        rows = allow_resizing_newdf[] ? 1 : nrow(df)
    else
        # allow squashing a scalar to 0 rows
        rows = nrow(newdf)
    end
    @assert length(colnames) == length(res)
    for (newname, v) in zip(colnames, res)
        # note that newdf potentially can contain newname in general
        newdf[!, newname] = fill!(Tables.allocatecolumn(typeof(v), rows), v)
    end
end

function _fix_existing_columns_for_vector(newdf::DataFrame, df::AbstractDataFrame,
                                          allow_resizing_newdf::Ref{Bool}, lr::Int,
                                          @nospecialize(fun))
    # allow shortening to 0 rows
    if allow_resizing_newdf[] && nrow(newdf) == 1
        newdfcols = _columns(newdf)
        for (i, col) in enumerate(newdfcols)
            newdfcols[i] = fill!(similar(col, lr), first(col))
        end
    end
    # !allow_resizing_newdf[] && ncol(newdf) == 0
    # means that we use `select` or `transform` not `combine`
    if !allow_resizing_newdf[] && ncol(newdf) == 0 && lr != nrow(df)
        throw(ArgumentError("length $(lr) of vector returned from " *
                            "function $fun is different from number of rows " *
                            "$(nrow(df)) of the source data frame."))
    end
    allow_resizing_newdf[] = false
end

function _add_col_check_copy(newdf::DataFrame, df::AbstractDataFrame,
                             col_idx::Union{Nothing, Int, AbstractVector{Int}, AsTable},
                             copycols::Bool, @nospecialize(fun),
                             newname::Symbol, v::AbstractVector)
    cdf = eachcol(df)
    vpar = parent(v)
    parent_cols = col_idx isa AsTable ? col_idx.cols : something(col_idx, 1:ncol(df))
    if copycols && !(fun isa ByRow) && (v isa SubArray || any(i -> vpar === parent(cdf[i]), parent_cols))
        newdf[!, newname] = copy(v)
    else
        newdf[!, newname] = v
    end
end

function _add_multicol_res(res::AbstractDataFrame, newdf::DataFrame, df::AbstractDataFrame,
                           colnames::AbstractVector{Symbol},
                           allow_resizing_newdf::Ref{Bool}, @nospecialize(fun),
                           col_idx::Union{Nothing, Int, AbstractVector{Int}, AsTable},
                           copycols::Bool, newname::Union{Nothing, Type{AsTable}, AbstractVector{Symbol}})
    lr = nrow(res)
    _fix_existing_columns_for_vector(newdf, df, allow_resizing_newdf, lr, fun)
    @assert length(colnames) == ncol(res)
    for (newname, v) in zip(colnames, eachcol(res))
        _add_col_check_copy(newdf, df, col_idx, copycols, fun, newname, v)
    end
end

function _add_multicol_res(res::AbstractMatrix, newdf::DataFrame, df::AbstractDataFrame,
                           colnames::AbstractVector{Symbol},
                           allow_resizing_newdf::Ref{Bool}, @nospecialize(fun),
                           col_idx::Union{Nothing, Int, AbstractVector{Int}, AsTable},
                           copycols::Bool, newname::Union{Nothing, Type{AsTable}, AbstractVector{Symbol}})
    lr = size(res, 1)
    _fix_existing_columns_for_vector(newdf, df, allow_resizing_newdf, lr, fun)
    @assert length(colnames) == size(res, 2)
    for (i, newname) in enumerate(colnames)
        newdf[!, newname] = res[:, i]
    end
end

function _add_multicol_res(res::NamedTuple{<:Any, <:Tuple{Vararg{AbstractVector}}},
                           newdf::DataFrame, df::AbstractDataFrame,
                           colnames::AbstractVector{Symbol},
                           allow_resizing_newdf::Ref{Bool}, @nospecialize(fun),
                           col_idx::Union{Nothing, Int, AbstractVector{Int}, AsTable},
                           copycols::Bool, newname::Union{Nothing, Type{AsTable}, AbstractVector{Symbol}})
    lr = length(res[1])
    _fix_existing_columns_for_vector(newdf, df, allow_resizing_newdf, lr, fun)
    @assert length(colnames) == length(res)
    for (newname, v) in zip(colnames, res)
        _add_col_check_copy(newdf, df, col_idx, copycols, fun, newname, v)
    end
end

function _add_multicol_res(res::NamedTuple, newdf::DataFrame, df::AbstractDataFrame,
                           colnames::AbstractVector{Symbol},
                           allow_resizing_newdf::Ref{Bool}, @nospecialize(fun),
                           col_idx::Union{Nothing, Int, AbstractVector{Int}, AsTable},
                           copycols::Bool, newname::Union{Nothing, Type{AsTable}, AbstractVector{Symbol}})
    if any(v -> v isa AbstractVector, res)
        throw(ArgumentError("mixing single values and vectors in a named tuple is not allowed"))
    else
        _insert_row_multicolumn(newdf, df, allow_resizing_newdf, colnames, res)
    end
end

function _add_multicol_res(res::DataFrameRow, newdf::DataFrame, df::AbstractDataFrame,
                           colnames::AbstractVector{Symbol},
                           allow_resizing_newdf::Ref{Bool}, @nospecialize(fun),
                           col_idx::Union{Nothing, Int, AbstractVector{Int}, AsTable},
                           copycols::Bool, newname::Union{Nothing, Type{AsTable}, AbstractVector{Symbol}})
    _insert_row_multicolumn(newdf, df, allow_resizing_newdf, colnames, res)
end

function select_transform!(@nospecialize(nc::Union{Base.Callable, Pair{<:Union{Int, AbstractVector{Int}, AsTable},
                                                                       <:Pair{<:Base.Callable,
                                                                              <:Union{Symbol,
                                                                                AbstractVector{Symbol},
                                                                                DataType}}}}),
                           df::AbstractDataFrame, newdf::DataFrame,
                           transformed_cols::Set{Symbol}, copycols::Bool,
                           allow_resizing_newdf::Ref{Bool})
    if nc isa Base.Callable
        col_idx, fun, newname = nothing, nc, nothing
    else
        col_idx, (fun, newname) = nc
    end
    if newname isa DataType
        newname === AsTable || throw(ArgumentError("Only DataType supported as target is AsTable"))
    end
    # It is allowed to request a tranformation operation into a newname column
    # only once. This is ensured by the logic related to transformed_cols dictionaly
    # in _manipulate, therefore in select_transform! such a duplicate should not happen
    res = _transformation_helper(df, col_idx, fun)

    if newname === AsTable || newname isa AbstractVector{Symbol}
        res = _expand_to_table(res)
    end

    if res isa Union{AbstractDataFrame, NamedTuple, DataFrameRow, AbstractMatrix}
        if newname isa Symbol
            throw(ArgumentError("Table returned but a single output column was expected"))
        end
        colnames = _gen_colnames(res, newname)
        isempty(colnames) && return # nothing to do

        if any(in(transformed_cols), colnames)
            throw(ArgumentError("Duplicate column name(s) returned: :" *
                                "$(join(intersect(colnames, transformed_cols), ", :"))"))
        else
            startlen = length(transformed_cols)
            union!(transformed_cols, colnames)
            @assert startlen + length(colnames) == length(transformed_cols)
        end
        _add_multicol_res(res, newdf, df, colnames, allow_resizing_newdf, fun,
                          col_idx, copycols, newname)
    elseif res isa AbstractVector
        if newname === nothing
            newname = :x1
        end
        if newname in transformed_cols
            throw(ArgumentError("duplicate output column name: :$newname"))
        else
            push!(transformed_cols, newname)
        end
        lr = length(res)
        _fix_existing_columns_for_vector(newdf, df, allow_resizing_newdf, lr, fun)
        _add_col_check_copy(newdf, df, col_idx, copycols, fun, newname, res)
    else
        if newname === nothing
            newname = :x1
        end
        if newname in transformed_cols
            throw(ArgumentError("duplicate output column name: :$newname"))
        else
            push!(transformed_cols, newname)
        end
        res_unwrap = res isa Union{AbstractArray{<:Any, 0}, Ref} ? res[] : res
        if ncol(newdf) == 0
            # if allow_resizing_newdf[] is false we know this is select or transform
            rows = allow_resizing_newdf[] ? 1 : nrow(df)
        else
            # allow squashing a scalar to 0 rows
            rows = nrow(newdf)
        end
        newdf[!, newname] = fill!(Tables.allocatecolumn(typeof(res_unwrap), rows),
                                  res_unwrap)
    end
end

SELECT_ARG_RULES =
    """
    Arguments passed as `args...` can be:

    * Any index that is allowed for column indexing
      ($COLUMNINDEX_STR; $MULTICOLUMNINDEX_STR).
    * A function or a type
    * Column transformation operations using the `Pair` notation that is
      described below and vectors or matrices of such pairs.

    Columns can be renamed using the `old_column => new_column_name` syntax, and
    transformed using the `old_column => fun => new_column_name` syntax.
    `new_column_name` must be a `Symbol` or a string, a vector of `Symbol`s or
    strings, or `AsTable`. `fun` must be a function or a type. If `old_column` is a
    `Symbol`, a string, or an integer then `fun` is applied to the corresponding
    column vector. Otherwise `old_column` can be any column indexing syntax, in
    which case `fun` will be passed the column vectors specified by `old_column`
    as separate arguments. The only exception is when `old_column` is an
    `AsTable` type wrapping a selector, in which case `fun` is passed a
    `NamedTuple` containing the selected columns.

    Column renaming and transformation operations can be passed wrapped in
    vectors or matrices (this is useful when combined with broadcasting).

    # Rules when `new_column_name` is a `Symbol` or a string or is absent

    If `fun` returns a value of type other than `AbstractVector` then it will be
    repeated in a vector matching the target number of rows in the data
    frame, unless its type is one of `AbstractDataFrame`, `NamedTuple`,
    `DataFrameRow`, `AbstractMatrix`, in which case an error is thrown. As a
    particular rule, values wrapped in a `Ref` or a `0`-dimensional
    `AbstractArray` are unwrapped and then repeated.

    To apply `fun` to each row instead of whole columns, it can be wrapped in a
    `ByRow` struct. In this case if `old_column` is a `Symbol`, a string, or an
    integer then `fun` is applied to each element (row) of `old_column` using
    broadcasting. Otherwise `old_column` can be any column indexing syntax, in
    which case `fun` will be passed one argument for each of the columns
    specified by `old_column`. If `ByRow` is used it is allowed for
    `old_column` to select an empty set of columns, in which case `fun`
     is called for each row without any arguments.

    Column transformation can also be specified using the short `old_column =>
    fun` form. In this case, `new_column_name` is automatically generated as
    `\$(old_column)_\$(fun)` if `renamecols=true` and `\$(old_column)` if
    `renamecols=false`. Up to three column names are used for multiple input
    columns and they are joined using `_`; if more than three columns are passed
    then the name consists of the first two names and `etc` suffix then, e.g.
    `[:a,:b,:c,:d] => fun` produces the new column name `:a_b_etc_fun` if
    `renamecols=true` and ``:a_b_etc` if `renamecols=false`.
    It is not allowed to pass `renamecols=false` if `old_column` is empty
    as it would generate an empty column name.

    # Rules when `new_column_name` is a vector of `Symbol`s or strings or is `AsTable`

    In this case it is assumed that `fun` returns multiple columns.

    If `fun` returns one of `AbstractDataFrame`, `NamedTuple`, `DataFrameRow`,
    `AbstractMatrix` then rules described in the section describing the case
    when `args` is a function or a type apply.

    If `fun` returns an `AbstractVector` then each element of this vector must
    support the `keys` function, which must return a collection of `Symbol`s, strings
    or integers; the return value of `keys` must be identical for all elements.
    Then as many columns are created as there are elements in the return value
    of the `keys` function. If `new_column_name` is `AsTable` then their names
    are set to be equal to the key names except if `keys` returns integers, in
    which case they are prefixed by `x` (so the column names are e.g. `x1`,
    `x2`, ...). If `new_column_name` is a vector of `Symbol`s or strings then
    column names produced using the rules above are ignored and replaced by
    `new_column_name` (the number of columns must be the same as the length of
    `new_column_name` in this case).

    If `fun` returns a value of any other type then it is assumed that it is a
    table conforming to the Tables.jl API and the `Tables.columntable` function
    is called on it to get the resulting columns and their names. The names are
    retained when `new_column_name` is `AsTable` and are replaced if
    `new_column_name` is a vector of `Symbol`s or strings.

    # Rules when element of `args` is a function or a type

    In this case the function or type is called with `df` as a single argument.

    If the return value of the transformation is one of `AbstractDataFrame`,
    `NamedTuple`, `DataFrameRow` or `AbstractMatrix` then it is treated as
    containing multiple columns. For `AbstractMatrix` column names are generated
    as `x1`, `x2`, etc. For `AbstractDataFrame`, `NamedTuple` of vectors and
    `AbstractMatrix` the columns are taken as is from the returned value. For
    `DataFrameRow` and` NamedTuple` not containing any vectors the returned
    value is broadcasted to a vector matching the target number of rows in the data
    frame.

    If the return value is an `AbstractVector` then it is used as-is. The resulting
    column gets the name `x1`.

    In all other cases the return value is repeated in a vector matching
    the target number of rows in the data frame. As a particular rule, values
    wrapped in a `Ref` or a `0`-dimensional `AbstractArray` are unwrapped and
    then repeated. The resulting column gets the name `x1`.

    # Special rules

    As a special rule passing `nrow` without specifying `old_column` creates a
    column named `:nrow` containing a number of rows in a source data frame, and
    passing `nrow => new_column_name` stores the number of rows in source data
    frame in `new_column_name` column.

    If a collection of column names is passed to `select!` or `select` then
    requesting duplicate column names in target data frame are accepted (e.g.
    `select!(df, [:a], :, r"a")` is allowed) and only the first occurrence is
    used. In particular a syntax to move column `:col` to the first position in
    the data frame is `select!(df, :col, :)`. On the contrary, output column
    names of renaming, transformation and single column selection operations
    must be unique, so e.g. `select!(df, :a, :a => :a)` or
    `select!(df, :a, :a => ByRow(sin) => :a)` are not allowed.
    """

"""
    select!(df::DataFrame, args...; renamecols::Bool=true)
    select!(args::Callable, df::DataFrame; renamecols::Bool=true)

Mutate `df` in place to retain only columns specified by `args...` and return it.
The result is guaranteed to have the same number of rows as `df`, except when no
columns are selected (in which case the result has zero rows).

$SELECT_ARG_RULES

Note that including the same column several times in the data frame via renaming
or transformations that return the same object without copying will create
column aliases. An example of such a situation is
`select!(df, :a, :a => :b, :a => identity => :c)`.

# Examples
```jldoctest
julia> df = DataFrame(a=1:3, b=4:6)
3×2 DataFrame
│ Row │ a     │ b     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 4     │
│ 2   │ 2     │ 5     │
│ 3   │ 3     │ 6     │

julia> select!(df, 2)
3×1 DataFrame
│ Row │ b     │
│     │ Int64 │
├─────┼───────┤
│ 1   │ 4     │
│ 2   │ 5     │
│ 3   │ 6     │

julia> df = DataFrame(a=1:3, b=4:6);

julia> select!(df, :a => ByRow(sin) => :c, :b)
3×2 DataFrame
│ Row │ c        │ b     │
│     │ Float64  │ Int64 │
├─────┼──────────┼───────┤
│ 1   │ 0.841471 │ 4     │
│ 2   │ 0.909297 │ 5     │
│ 3   │ 0.14112  │ 6     │

julia> select!(df, :, [:c, :b] => (c,b) -> c .+ b .- sum(b)/length(b))
3×3 DataFrame
│ Row │ c        │ b     │ c_b_function │
│     │ Float64  │ Int64 │ Float64      │
├─────┼──────────┼───────┼──────────────┤
│ 1   │ 0.841471 │ 4     │ -0.158529    │
│ 2   │ 0.909297 │ 5     │ 0.909297     │
│ 3   │ 0.14112  │ 6     │ 1.14112      │

julia> df = DataFrame(a=1:3, b=4:6);

julia> select!(df, names(df) .=> [minimum maximum]);

julia> df
3×4 DataFrame
│ Row │ a_minimum │ b_minimum │ a_maximum │ b_maximum │
│     │ Int64     │ Int64     │ Int64     │ Int64     │
├─────┼───────────┼───────────┼───────────┼───────────┤
│ 1   │ 1         │ 4         │ 3         │ 6         │
│ 2   │ 1         │ 4         │ 3         │ 6         │
│ 3   │ 1         │ 4         │ 3         │ 6         │

julia> df = DataFrame(a=1:3, b=4:6);

julia> using Statistics

julia> select!(df, AsTable(:) => ByRow(mean), renamecols=false)
3×1 DataFrame
│ Row │ a_b     │
│     │ Float64 │
├─────┼─────────┤
│ 1   │ 2.5     │
│ 2   │ 3.5     │
│ 3   │ 4.5     │

julia> df = DataFrame(a=1:3, b=4:6);

julia> select!(first, df)
3×2 DataFrame
│ Row │ a     │ b     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 4     │
│ 2   │ 1     │ 4     │
│ 3   │ 1     │ 4     │

julia> df = DataFrame(a=1:3, b=4:6, c=7:9)
3×3 DataFrame
│ Row │ a     │ b     │ c     │
│     │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┤
│ 1   │ 1     │ 4     │ 7     │
│ 2   │ 2     │ 5     │ 8     │
│ 3   │ 3     │ 6     │ 9     │

julia> select!(df, AsTable(:) => ByRow(x -> (mean=mean(x), std=std(x))) => :stats,
               AsTable(:) => ByRow(x -> (mean=mean(x), std=std(x))) => AsTable)
3×3 DataFrame
│ Row │ stats                   │ mean    │ std     │
│     │ NamedTuple…             │ Float64 │ Float64 │
├─────┼─────────────────────────┼─────────┼─────────┤
│ 1   │ (mean = 4.0, std = 3.0) │ 4.0     │ 3.0     │
│ 2   │ (mean = 5.0, std = 3.0) │ 5.0     │ 3.0     │
│ 3   │ (mean = 6.0, std = 3.0) │ 6.0     │ 3.0     │
```

"""
select!(df::DataFrame, @nospecialize(args...); renamecols::Bool=true) =
    _replace_columns!(df, select(df, args..., copycols=false, renamecols=renamecols))

function select!(arg::Base.Callable, df::AbstractDataFrame; renamecols::Bool=true)
    if arg isa Colon
        throw(ArgumentError("First argument must be a transformation if the second argument is a data frame"))
    end
    return select!(df, arg)
end

"""
    transform!(df::DataFrame, args...; renamecols::Bool=true)
    transform!(args::Callable, df::DataFrame; renamecols::Bool=true)

Mutate `df` in place to add columns specified by `args...` and return it.
The result is guaranteed to have the same number of rows as `df`.
Equivalent to `select!(df, :, args...)`.

See [`select!`](@ref) for detailed rules regarding accepted values for `args`.
"""
transform!(df::DataFrame, @nospecialize(args...); renamecols::Bool=true) =
    select!(df, :, args..., renamecols=renamecols)

function transform!(arg::Base.Callable, df::AbstractDataFrame; renamecols::Bool=true)
    if arg isa Colon
        throw(ArgumentError("First argument must be a transformation if the second argument is a data frame"))
    end
    return transform!(df, arg)
end

"""
    select(df::AbstractDataFrame, args...; copycols::Bool=true, renamecols::Bool=true)
    select(args::Callable, df::DataFrame; renamecols::Bool=true)

Create a new data frame that contains columns from `df` specified by `args` and
return it. The result is guaranteed to have the same number of rows as `df`,
except when no columns are selected (in which case the result has zero rows)..

If `df` is a `DataFrame` or `copycols=true` then column renaming and transformations
are supported.

$SELECT_ARG_RULES

If `df` is a `DataFrame` a new `DataFrame` is returned.
If `copycols=false`, then the returned `DataFrame` shares column vectors with `df`
where possible.
If `copycols=true` (the default), then the returned `DataFrame` will not share
columns with `df`.
The only exception for this rule is the `old_column => fun => new_column`
transformation when `fun` returns a vector that is not allocated by `fun` but is
neither a `SubArray` nor one of the input vectors.
In such a case a new `DataFrame` might contain aliases. Such a situation can
only happen with transformations which returns vectors other than their inputs,
e.g. with `select(df, :a => (x -> c) => :c1, :b => (x -> c) => :c2)`  when `c`
is a vector object or with `select(df, :a => (x -> df.c) => :c2)`.

If `df` is a `SubDataFrame` and `copycols=true` then a `DataFrame` is returned
and the same copying rules apply as for a `DataFrame` input:
this means in particular that selected columns will be copied.
If `copycols=false`, a `SubDataFrame` is returned without copying columns.

Note that including the same column several times in the data frame via renaming
or transformations that return the same object when `copycols=false` will create
column aliases. An example of such a situation is
`select(df, :a, :a => :b, :a => identity => :c, copycols=false)`.

# Examples
```jldoctest
julia> df = DataFrame(a=1:3, b=4:6)
3×2 DataFrame
│ Row │ a     │ b     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 4     │
│ 2   │ 2     │ 5     │
│ 3   │ 3     │ 6     │

julia> select(df, 2)
3×1 DataFrame
│ Row │ b     │
│     │ Int64 │
├─────┼───────┤
│ 1   │ 4     │
│ 2   │ 5     │
│ 3   │ 6     │

julia> select(df, :a => ByRow(sin) => :c, :b)
3×2 DataFrame
│ Row │ c        │ b     │
│     │ Float64  │ Int64 │
├─────┼──────────┼───────┤
│ 1   │ 0.841471 │ 4     │
│ 2   │ 0.909297 │ 5     │
│ 3   │ 0.14112  │ 6     │

julia> select(df, :, [:a, :b] => (a,b) -> a .+ b .- sum(b)/length(b))
3×3 DataFrame
│ Row │ a     │ b     │ a_b_function │
│     │ Int64 │ Int64 │ Float64      │
├─────┼───────┼───────┼──────────────┤
│ 1   │ 1     │ 4     │ 0.0          │
│ 2   │ 2     │ 5     │ 2.0          │
│ 3   │ 3     │ 6     │ 4.0          │

julia> select(df, names(df) .=> [minimum maximum])
3×4 DataFrame
│ Row │ a_minimum │ b_minimum │ a_maximum │ b_maximum │
│     │ Int64     │ Int64     │ Int64     │ Int64     │
├─────┼───────────┼───────────┼───────────┼───────────┤
│ 1   │ 1         │ 4         │ 3         │ 6         │
│ 2   │ 1         │ 4         │ 3         │ 6         │
│ 3   │ 1         │ 4         │ 3         │ 6         │

julia> using Statistics

julia> select(df, AsTable(:) => ByRow(mean), renamecols=false)
3×1 DataFrame
│ Row │ a_b     │
│     │ Float64 │
├─────┼─────────┤
│ 1   │ 2.5     │
│ 2   │ 3.5     │
│ 3   │ 4.5     │

julia> select(first, df)
3×2 DataFrame
│ Row │ a     │ b     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 4     │
│ 2   │ 1     │ 4     │
│ 3   │ 1     │ 4     │

julia> df = DataFrame(a=1:3, b=4:6, c=7:9)
3×3 DataFrame
│ Row │ a     │ b     │ c     │
│     │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┤
│ 1   │ 1     │ 4     │ 7     │
│ 2   │ 2     │ 5     │ 8     │
│ 3   │ 3     │ 6     │ 9     │

julia> select(df, AsTable(:) => ByRow(x -> (mean=mean(x), std=std(x))) => :stats,
              AsTable(:) => ByRow(x -> (mean=mean(x), std=std(x))) => AsTable)
3×3 DataFrame
│ Row │ stats                   │ mean    │ std     │
│     │ NamedTuple…             │ Float64 │ Float64 │
├─────┼─────────────────────────┼─────────┼─────────┤
│ 1   │ (mean = 4.0, std = 3.0) │ 4.0     │ 3.0     │
│ 2   │ (mean = 5.0, std = 3.0) │ 5.0     │ 3.0     │
│ 3   │ (mean = 6.0, std = 3.0) │ 6.0     │ 3.0     │
```

"""
select(df::AbstractDataFrame, @nospecialize(args...); copycols::Bool=true, renamecols::Bool=true) =
    manipulate(df, args..., copycols=copycols, keeprows=true, renamecols=renamecols)

function select(arg::Base.Callable, df::AbstractDataFrame; renamecols::Bool=true)
    if arg isa Colon
        throw(ArgumentError("First argument must be a transformation if the second argument is a data frame"))
    end
    return select(df, arg)
end

"""
    transform(df::AbstractDataFrame, args...; copycols::Bool=true, renamecols::Bool=true)
    transform(args::Callable, df::DataFrame; renamecols::Bool=true)

Create a new data frame that contains columns from `df` and adds columns
specified by `args` and return it.
The result is guaranteed to have the same number of rows as `df`.
Equivalent to `select(df, :, args..., copycols=copycols)`.

See [`select`](@ref) for detailed rules regarding accepted values for `args`.
"""
transform(df::AbstractDataFrame, @nospecialize(args...); copycols::Bool=true, renamecols::Bool=true) =
    select(df, :, args..., copycols=copycols, renamecols=renamecols)

function transform(arg::Base.Callable, df::AbstractDataFrame; renamecols::Bool=true)
    if arg isa Colon
        throw(ArgumentError("First argument to must be a transformation if the second argument is a data frame"))
    end
    return transform(df, arg)
end

"""
    combine(df::AbstractDataFrame, args...; renamecols::Bool=true)
    combine(args::Callable, df::AbstractDataFrame; renamecols::Bool=true)

Create a new data frame that contains columns from `df` specified by `args` and
return it. The result can have any number of rows that is determined by the
values returned by passed transformations.

See [`select`](@ref) for detailed rules regarding accepted values for `args` in
`combine(df, args...)` form. For `combine(arg, df)` the same rules as for
`combine` on `GroupedDataFrame` apply except that a `df` with zero rows is
currently not allowed.

# Examples
```jldoctest
julia> df = DataFrame(a=1:3, b=4:6)
3×2 DataFrame
│ Row │ a     │ b     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 4     │
│ 2   │ 2     │ 5     │
│ 3   │ 3     │ 6     │

julia> combine(df, :a => sum, nrow, renamecols=false)
1×2 DataFrame
│ Row │ a     │ nrow  │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 6     │ 3     │

julia> combine(df, :a => ByRow(sin) => :c, :b)
3×2 DataFrame
│ Row │ c        │ b     │
│     │ Float64  │ Int64 │
├─────┼──────────┼───────┤
│ 1   │ 0.841471 │ 4     │
│ 2   │ 0.909297 │ 5     │
│ 3   │ 0.14112  │ 6     │

julia> combine(df, :, [:a, :b] => (a,b) -> a .+ b .- sum(b)/length(b))
3×3 DataFrame
│ Row │ a     │ b     │ a_b_function │
│     │ Int64 │ Int64 │ Float64      │
├─────┼───────┼───────┼──────────────┤
│ 1   │ 1     │ 4     │ 0.0          │
│ 2   │ 2     │ 5     │ 2.0          │
│ 3   │ 3     │ 6     │ 4.0          │

julia> combine(df, names(df) .=> [minimum maximum])
1×4 DataFrame
│ Row │ a_minimum │ b_minimum │ a_maximum │ b_maximum │
│     │ Int64     │ Int64     │ Int64     │ Int64     │
├─────┼───────────┼───────────┼───────────┼───────────┤
│ 1   │ 1         │ 4         │ 3         │ 6         │

julia> using Statistics

julia> combine(df, AsTable(:) => ByRow(mean), renamecols=false)
3×1 DataFrame
│ Row │ a_b     │
│     │ Float64 │
├─────┼─────────┤
│ 1   │ 2.5     │
│ 2   │ 3.5     │
│ 3   │ 4.5     │

julia> combine(first, df)
1×2 DataFrame
│ Row │ a     │ b     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 4     │

julia> df = DataFrame(a=1:3, b=4:6, c=7:9)
3×3 DataFrame
│ Row │ a     │ b     │ c     │
│     │ Int64 │ Int64 │ Int64 │
├─────┼───────┼───────┼───────┤
│ 1   │ 1     │ 4     │ 7     │
│ 2   │ 2     │ 5     │ 8     │
│ 3   │ 3     │ 6     │ 9     │

julia> combine(df, AsTable(:) => ByRow(x -> (mean=mean(x), std=std(x))) => :stats,
               AsTable(:) => ByRow(x -> (mean=mean(x), std=std(x))) => AsTable)
3×3 DataFrame
│ Row │ stats                   │ mean    │ std     │
│     │ NamedTuple…             │ Float64 │ Float64 │
├─────┼─────────────────────────┼─────────┼─────────┤
│ 1   │ (mean = 4.0, std = 3.0) │ 4.0     │ 3.0     │
│ 2   │ (mean = 5.0, std = 3.0) │ 5.0     │ 3.0     │
│ 3   │ (mean = 6.0, std = 3.0) │ 6.0     │ 3.0     │
```
"""
combine(df::AbstractDataFrame, @nospecialize(args...); renamecols::Bool=true) =
    manipulate(df, args..., copycols=true, keeprows=false, renamecols=renamecols)

function combine(arg::Base.Callable, df::AbstractDataFrame; renamecols::Bool=true)
    if arg isa Colon
        throw(ArgumentError("First argument to select! must be a transformation if the second argument is a data frame"))
    end
    return combine(df, arg)
end

manipulate(df::DataFrame, args::AbstractVector{Int}; copycols::Bool, keeprows::Bool,
           renamecols::Bool) =
    DataFrame(_columns(df)[args], Index(_names(df)[args]), copycols=copycols)

function manipulate(df::DataFrame, c::MultiColumnIndex; copycols::Bool, keeprows::Bool,
                    renamecols::Bool)
    if c isa AbstractVector{<:Pair}
        return manipulate(df, c..., copycols=copycols, keeprows=keeprows,
                          renamecols=renamecols)
    else
        return manipulate(df, index(df)[c], copycols=copycols, keeprows=keeprows,
                          renamecols=renamecols)
    end
end

manipulate(df::DataFrame, c::ColumnIndex; copycols::Bool, keeprows::Bool,
           renamecols::Bool) =
    manipulate(df, [c], copycols=copycols, keeprows=keeprows, renamecols=renamecols)

function manipulate(df::DataFrame, @nospecialize(cs...); copycols::Bool, keeprows::Bool, renamecols::Bool)
    cs_vec = []
    for v in cs
        if v isa AbstractVecOrMat{<:Pair}
            append!(cs_vec, v)
        else
            push!(cs_vec, v)
        end
    end
    return _manipulate(df, [normalize_selection(index(df), c, renamecols) for c in cs_vec],
                    copycols, keeprows)
end

function _manipulate(df::AbstractDataFrame, @nospecialize(normalized_cs), copycols::Bool, keeprows::Bool)
    @assert !(df isa SubDataFrame && copycols==false)
    newdf = DataFrame()
    # the role of transformed_cols is the following
    # * make sure that we do not use the same target column name twice in transformations;
    #   note though that it can appear in no-transformation selection like
    #   `select(df, :, :a => ByRow(sin) => :a), where :a is produced both by `:`
    #   and by `:a => ByRow(sin) => :a`
    # * make sure that if some column is produced by transformation like
    #   `:a => ByRow(sin) => :a` and it appears earlier or later in non-transforming
    #   selection like `:` or `:a` then the transformation is computed and inserted
    #   in to the target data frame once and only once the first time the target column
    #   is requested to be produced.
    #
    # For example in:
    #
    # julia> df = DataFrame(a=1:2, b=3:4)
    # 2×2 DataFrame
    # │ Row │ a     │ b     │
    # │     │ Int64 │ Int64 │
    # ├─────┼───────┼───────┤
    # │ 1   │ 1     │ 3     │
    # │ 2   │ 2     │ 4     │
    #
    # julia> select(df, :, :a => ByRow(sin) => :a, :a, 1)
    # 2×2 DataFrame
    # │ Row │ a        │ b     │
    # │     │ Float64  │ Int64 │
    # ├─────┼──────────┼───────┤
    # │ 1   │ 0.841471 │ 3     │
    # │ 2   │ 0.909297 │ 4     │
    #
    # transformed_cols keeps a set of columns that were generated via a transformation
    # up till the point. Note that single column selection and column renaming is
    # considered to be a transformation
    transformed_cols = Set{Symbol}()
    # we allow resizing newdf only if up to some point only scalars were put
    # in it. The moment we put any vector into newdf its number of rows becomes fixed
    # Also if keeprows is true then we make sure to produce nrow(df) rows so resizing
    # is not allowed
    allow_resizing_newdf = Ref(!keeprows)
    for nc in normalized_cs
        if nc isa AbstractVector{Int} # only this case is NOT considered to be a transformation
            allunique(nc) || throw(ArgumentError("duplicate column names selected"))
            for i in nc
                newname = _names(df)[i]
                # as nc is a multiple column selection without transformations
                # we allow duplicate column names with selections applied earlier
                # and ignore them for convinience, to allow for e.g. select(df, :x1, :)
                if !hasproperty(newdf, newname)
                    # allow shortening to 0 rows
                    if allow_resizing_newdf[] && nrow(newdf) == 1
                        newdfcols = _columns(newdf)
                        for (i, col) in enumerate(newdfcols)
                            newdfcols[i] = fill!(similar(col, nrow(df)), first(col))
                        end
                    end
                    # here even if keeprows is true all is OK
                    newdf[!, newname] = copycols ? df[:, i] : df[!, i]
                    allow_resizing_newdf[] = false
                end
            end
        else
            select_transform!(nc, df, newdf, transformed_cols, copycols,
                              allow_resizing_newdf)
        end
    end
    return newdf
end

manipulate(dfv::SubDataFrame, ind::ColumnIndex; copycols::Bool, keeprows::Bool,
           renamecols::Bool) =
    manipulate(dfv, [ind], copycols=copycols, keeprows=keeprows, renamecols=renamecols)

function manipulate(dfv::SubDataFrame, args::MultiColumnIndex;
                 copycols::Bool, keeprows::Bool, renamecols::Bool)
    if args isa AbstractVector{<:Pair}
        return manipulate(dfv, args..., copycols=copycols, keeprows=keeprows,
                          renamecols=renamecols)
    else
        return copycols ? dfv[:, args] : view(dfv, :, args)
    end
end

function manipulate(dfv::SubDataFrame, @nospecialize(args...); copycols::Bool, keeprows::Bool,
                    renamecols::Bool)
    if copycols
        cs_vec = []
        for v in args
            if v isa AbstractVector{<:Pair}
                append!(cs_vec, v)
            else
                push!(cs_vec, v)
            end
        end
        return _manipulate(dfv, [normalize_selection(index(dfv), c, renamecols) for c in cs_vec],
                           true, keeprows)
    else
        # we do not support transformations here
        # newinds contains only indexing; making it Vector{Any} avoids some compilation
        newinds = []
        seen_single_column = Set{Int}()
        for ind in args
            if ind isa ColumnIndex
                ind_idx = index(dfv)[ind]
                if ind_idx in seen_single_column
                    throw(ArgumentError("selecting the same column multiple times " *
                                        "using Symbol, string or integer is not allowed " *
                                        "($ind was passed more than once"))
                else
                    push!(seen_single_column, ind_idx)
                end
            else
                newind = normalize_selection(index(dfv), ind, renamecols)
                if newind isa Pair
                    throw(ArgumentError("transforming and renaming columns of a " *
                                        "SubDataFrame is not allowed when `copycols=false`"))
                end
                push!(newinds, newind)
            end
        end
        return view(dfv, :, isempty(newinds) ? [] : All(newinds...))
    end
end


text:
mutable struct ElementaryPathLightGraphsDijkstraSolver <: ElementaryPathSolver
  graph::SimpleDiGraph
  weights_matrix::Matrix{Float64}
  source::Int
  destination::Int

  function ElementaryPathLightGraphsDijkstraSolver()
    return new(DiGraph(0), zeros(0, 0), -1, -1)
  end
end

function build!(solver::ElementaryPathLightGraphsDijkstraSolver, graph::SimpleDiGraph, source::Int, destination::Int)
  n = nv(graph)
  solver.graph = graph
  solver.weights_matrix = zeros(n, n)
  solver.source = source
  solver.destination = destination
end

has_lp_formulation(::ElementaryPathLightGraphsDijkstraSolver) = false

function solve_linear(solver::ElementaryPathLightGraphsDijkstraSolver, rewards::Dict{Tuple{Int, Int}, Float64})
  # Make up a rewards matrix by copying the input dictionary into the right data structure for the shortest path computations.
  fill!(solver.weights_matrix, 0.0)
  for (k, v) in rewards
    solver.weights_matrix[k[1], k[2]] = v
  end
  solver.weights_matrix = maximum(solver.weights_matrix) .- solver.weights_matrix

  # Compute the shortest path.
  state = dijkstra_shortest_paths(solver.graph, solver.source, solver.weights_matrix)

  # Retrieve the shortest path based on the returned data structure (using the predecessor node).
  # Everything must be done in reverse, due to the way the path is stored for Dijkstra.
  path = Vector{Tuple{Int, Int}}()
  current_node = solver.destination
  while current_node != solver.source
    push!(path, (state.parents[current_node][1], current_node))
    current_node = state.parents[current_node][1]
  end
  return path
end


text:
using Statistics
using NumericalDataManipulation.Common
using NumericalDataManipulation.CommonMath

function merge_grids(G¹::DataGrid, G²::DataGrid;
    ntail::Int=5, ϵʳᵉˡ::Float64=0.4, meanΔ::Function=mean_log10)
    newknots = merge_knots(G¹.knots, G².knots; ntail=ntail, ϵʳᵉˡ=ϵʳᵉˡ, meanΔ=meanΔ)
    return DataGrid(newknots)
end
function merge_knots(V¹::Vector{Float64}, V²::Vector{Float64};
    ntail::Int=5, ϵʳᵉˡ::Float64=0.4, meanΔ::Function=mean_log10)
    @assert(issorted(V¹), "The first vector argument is not sorted in an ascending order.")
    @assert(issorted(V²), "The second vector argument is not sorted in an ascending order.")
    @assert(allunique(V¹), "The first vector argument must contain uniques elements.")
    @assert(allunique(V²), "The second vector argument must contain uniques elements.")
    V⁰ = vcat(V¹, V²)
    if !issorted(V⁰)
        sort!(V⁰)
    end
    if !allunique(V⁰)
        unique!(V⁰)
    end
    return unique_knots(V⁰; ntail=ntail, ϵʳᵉˡ=ϵʳᵉˡ, meanΔ=meanΔ)
end

function unique_knots(knots::Vector{Float64};
    ntail::Int=5, ϵʳᵉˡ::Float64=0.4, meanΔ::Function=mean_log10)
    @assert issorted(knots)
    @assert allunique(knots)
    @assert ntail ≥ 2
    @assert ϵʳᵉˡ > 0
    Δknots = Δknot_values(knots)
    @assert ntail < length(Δknots)

    Δ𝐆ˡᵒᵍ¹⁰ = meanΔ(nonzero(Δknots))
    @debug "Δ𝐆ˡᵒᵍ¹⁰ = $Δ𝐆ˡᵒᵍ¹⁰"
    @assert Δ𝐆ˡᵒᵍ¹⁰ ≠ 0

    dots = Vector{Float64}()
    K = length(knots)
    Δtail = Vector{Float64}()
    Δhead = Vector{Float64}()
    pushfirst!(Δtail, ∞)
    push!(Δhead, Δknots[1])
    for k = 1:K - 1
        # --- tail ---
        if k < ntail
            pushfirst!(Δtail, Δknots[k])
        elseif ntail ≤ k < K - 1
            pop!(Δtail)
            pushfirst!(Δtail, Δknots[k])
        else @assert k == K - 1
            pop!(Δtail)
            pushfirst!(Δtail, ∞)
        end
        @assert length(Δtail) ≥ 2
        # -------------
        # --- head ---
        if 1 ≤ k < K - 2
            pop!(Δhead)
            push!(Δhead, Δknots[k + 1])
        else
            pop!(Δhead)
            push!(Δhead, ∞)
        end
        @assert length(Δhead) == 1
        # ------------
        # --- current knots ---
        knotᵏ = knots[k]
        knotᵏ⁺¹ = knots[k + 1]
        # ---------------------

        @debug "--- BEGIN --- [$knotᵏ, $knotᵏ⁺¹]"
        @debug "-- begin --- ϵΔtail, ϵΔhead"
        ϵΔtail = relfinity(Δtail; D=Δ𝐆ˡᵒᵍ¹⁰, ϵʳᵉˡ=ϵʳᵉˡ)
        ϵΔhead = relfinity(Δhead; D=Δ𝐆ˡᵒᵍ¹⁰, ϵʳᵉˡ=ϵʳᵉˡ)
        @debug "ϵΔtail: $ϵΔtail"
        @debug "ϵΔhead: $ϵΔhead"
        @debug "--- end ---"
        if [1, ∞] == ϵΔtail[1:2]
            @debug "[1, ∞]"
            push!(dots, knotᵏ)
            @debug "[1, ∞] - pushed: $knotᵏ"
        elseif [0, ∞] == ϵΔtail[1:2]
            @debug "[0, ∞]"
            push!(dots, knotᵏ)
            @debug "[0, ∞] - pushed: $knotᵏ"
        elseif [1, 1] == ϵΔtail[1:2]
            @debug "[1, 1]"
            push!(dots, knotᵏ)
            @debug "[1, 1] - pushed: $knotᵏ"
        elseif [0, 1] == ϵΔtail[1:2]
            @debug "[0, 1]"
            if [1] == ϵΔhead
                @debug "[0, 1] - 1"
                if length(ϵΔtail) > 2 && ϵΔhead[1] == 1
                    push!(dots, (knotᵏ + knotᵏ⁺¹) / 2)
                    @debug "[0, 1] - 1 - pushed: $((knotᵏ + knotᵏ⁺¹) / 2)"
                else
                    @debug "[0, 1] - 1 - nothing pushed"
                end
            elseif [0] == ϵΔhead
                @debug "[0, 1] - 0"
                push!(dots, knotᵏ)
                @debug "[0, 1] - 0 - pushed: $knotᵏ"
            elseif [∞] == ϵΔhead
                @debug "[0, 1] - ∞"
                # The cases [0, ∞], [1, ∞] go next.
            end
        elseif [0, 0] == ϵΔtail[1:2]
            @debug "[0, 0]"
        elseif [1, 0] == ϵΔtail[1:2]
            @debug "[1, 0]"
            if length(ϵΔtail) == 2 || (length(ϵΔtail) > 2 && 0 == ϵΔtail[3])
                # TODO: maybe this branch must be the only one to calculate mean knots.
                last_dot = dots[end]
                Σδ, n = sum_until(
                    (δ::Float64, Σδ::Float64) ->
                        δ ∈ [∞, NaN, 1] || 0 ≠ relfinity(Σδ; D=Δ𝐆ˡᵒᵍ¹⁰, ϵʳᵉˡ=ϵʳᵉˡ),
                    Δtail[2:end])
                δ⁰⁵ = Σδ / 2
                if dots[end] > knotᵏ - δ⁰⁵
                    @warn "unsorted next knot: $(knotᵏ - δ⁰⁵), previous: $(dots[end]); tail: $Δtail, ϵΔtail: $ϵΔtail"
                end
                push!(dots, knotᵏ - δ⁰⁵)
                @debug "[1, 0] - pushed: $(knotᵏ - δ⁰⁵) (knotᵏ = $knotᵏ, δ⁰⁵ = $δ⁰⁵)"
            end
            # The cases [0, 1], [0, 0] were handled.
            # The case [1, 0] will be handled.
        elseif [∞, 0] == ϵΔtail[1:2]
            @debug "[∞, 0]"
            push!(dots, knotᵏ⁺¹)
            @debug "[∞, 0] - pushed: $knotᵏ⁺¹"
        elseif [∞, 1] == ϵΔtail[1:2]
            @debug "[∞, 1]"
            push!(dots, knotᵏ⁺¹)
            @debug "[∞, 1] - pushed: $knotᵏ⁺¹"
        end
        @debug "--- END ---"
    end
    return dots
end


text:
n2sStat(x::T where T<:Real) = num2Str(x, 3, Ints=true)

#calculates the test statistic and T value of a data vector
skewnessStat(data::Vector{Float64})::Float64 = skewness(data)/√(6. / length(data))
kurtosisStat(data::Vector{Float64})::Float64 = kurtosis(data)/√(24. / length(data))
JBStat(data::Vector{Float64})::Float64 =
  skewness(data)^2/(6. / length(data)) + kurtosis(data)^2/(24. / length(data))

#simple scheme for doing the ljungBox test
function ljungBox(ρ::Vector{Float64}, N::Int)::Float64
  m::Int = length(ρ)
  Q::Float64 = N*(N+2) * sum((ρ .* ρ) ./ (N .- collect(1:m)))
  return Q
end

#convenience function, defaults to log(T) for the lag
ljungBox(data::Vector{Float64}, lagRange::T where T <:AbstractRange = 1:(Int(round(log(length(data)))))) =
  ljungBox(autoρ(data, lagRange), length(data))

mutable struct UnivariateTest
  statistic::Function
  distribution::Distribution
  H₀::Float64
  testType::Symbol
end

#convenience constructor for univariate test
UnivariateTest(statistic::Function, distribution::Distribution;
  H₀::Float64=0.0, testType::Symbol=:upper)::UnivariateTest =
  UnivariateTest(statistic, distribution, H₀, testType)

function test(U::UnivariateTest, data::Vector{Float64}; verbose::Bool = true)::NTuple{2,Float64}
  θ::Float64 = U.statistic(data .- U.H₀)

  p::Float64 = cdf(U.distribution, θ)

  #adjust the p value for the specific test type
  if U.testType == :upper
    p = 1.0 - p
  elseif  U.testType == :lower
  elseif U.testType == :twoTailed
    p = p > 0.5 ? (1-p)*2 : 2*p
  else
    warn("Test type not found in univariate test")
  end

  return (θ, p)
end

#convenience function for test
test(U::UnivariateTest, data::T where T <: AbstractVector) = Vector{Float64}(dropmissing(data))

#simple function to get the autocorrelation
function autoρ(v::T, lag::Int)::Float64 where T <: AbstractVector
  vAvg::Float64 = mean(skipmissing(v))
  Δv::Vector{Float64} = v .- vAvg
  rng::S where S <:AbstractRange = (lag+1):length(Δv)
  vCov::Float64 = sum(((i::Int)->Δv[i]*Δv[i-lag]).(rng))
  return vCov/sum(Δv.^2.)
end

#convenience function to take in many lags
autoρ(v::T where T <: AbstractVector, lags::T where T<:AbstractRange)::Vector{Float64} = ((l::Int)->autoρ(v,l)).(lags)

#simple function to get the partial autocorrelation
function pAutoρ(v::T, lag::Int)::Float64 where T <: AbstractVector

  #fill the lag matrix
  Y::Vector{Float64} = v[(lag+1):end]
  X::Matrix{Float64} = ones(length(v)-lag,lag+1)
  for i::Int ∈ 1:lag
    X[:,i+1] .= v[(lag-i+1):(end-i)]
  end

  reg::FMLM = FMLM(X,Y)

  return reg.β[end]
end

#convenience function to take in many lags
pAutoρ(v::T where T <: AbstractVector, lags::T where T<:AbstractRange)::Vector{Float64} =
  ((l::Int)->pAutoρ(v,l)).(lags)

#Lag a series in a dataframe
function lagDF!(DF::DataFrame, target::Symbol, lags::Int = 1; name::Symbol = Symbol("$(target)L$lags"))::DataFrame
  DF[:,name] = similar(DF[:,target])
  DF[:,name] .= missing
  DF[(lags+1):end,name] = DF[1:(end-lags),target]

  return DF
end

function lagDF!(DF::DataFrame, targets::Vector{Symbol}, lags::Int = 1;
  names::Vector{Symbol}= ((s::Symbol)->Symbol("$(s)L$lags")).(targets))

  for i::Int ∈ 1:length(targets)
    lagDF!(DF, targets[i], lags, name = names[i])
  end

  return DF
end

mutable struct ARMA
  c::T where T <: Real
  ϕ::Vector{T}  where T <: Real
  ϑ::Vector{T} where T <: Real
  σ²::T where T <: Real
  dist::T where T<: Distribution
  p::Int
  q::Int
end

#constructor for ARMA
function ARMA(;c::Float64 = 0.0, ϕ::Vector{Float64}=ones(1), ϑ::Vector{Float64} = Vector{Float64}(),
    σ²::Float64 = 1.0, dist::UnivariateDistribution = Normal(0,σ²^0.5),
    p::Int = length(ϕ), q::Int = length(ϑ))::ARMA
  return ARMA(c,ϕ,ϑ,σ²,dist, p, q)
end

#this function allocates but does not paramaterize the model
function ARMA(p::Int, q::Int; c::Float64 = 0.0,
    σ²::Float64 = 1.0, dist::UnivariateDistribution = Normal(0,σ²^0.5))::ARMA

    return ARMA(c,zeros(p),zeros(q),σ²,dist, p, q)
end

#these are helper methods to change the type for things like autodifferentiation
function generalizeARMA!(Θ::ARMA)::Nothing
  Θ.c = Real(Θ.c)
  Θ.ϕ = Vector{Real}(Θ.ϕ)
  Θ.ϑ = Vector{Real}(Θ.ϑ)
  Θ.σ² = Real(Θ.σ²)

  return nothing
end

function narrowARMA!(Θ::ARMA)::Nothing
  Θ.c = Float64(Θ.c)
  Θ.ϕ = Vector{Float64}(Θ.ϕ)
  Θ.ϑ = Vector{Float64}(Θ.ϑ)
  Θ.σ² = Float64(Θ.σ²)

  return nothing
end

#fills the predicted Y values using the BJ methedology as described in Hamilton
function ARMAFillYHat!(Θ::ARMA, Y::Vector{Float64}, YHat::Vector{T} = similar(Y),
    ε::Vector{T} = similar(Y)
    #ε::Vector{T} = Vector{T}(length(Y)+max(Θ.q-Θ.p,0))
    )::NTuple{2,Vector{T}} where T<:Real

  ε .= 0.
  τ::Int = length(Y)

  #start with the observations to kick the model into gear
  YHat[1:Θ.p] .= Y[1:Θ.p]

  AR::T = 0.0
  MA::T = 0.0

  qIndexAdj::Int = max(Θ.q-Θ.p,0)

  #fill the predicted residual arrays
  for i::Int ∈ (Θ.p+1):τ
      for j::Int ∈ 1:Θ.p
          AR+=Y[i-j]*Θ.ϕ[j]
      end
      for j::Int ∈ 1:(min(Θ.q,i-1))
          MA -= ε[i-j]*Θ.ϑ[j]
      end
    #AR = Θ.p>0 ? sum(view(Y,(i-1):-1:(i-Θ.p)) .* Θ.ϕ):0.0
    #MA = Θ.q>0 ? sum( view(ε,(i-1+qIndexAdj):-1:(i-Θ.q+qIndexAdj)) .* Θ.ϑ):0.0
    #      println("got here")
    YHat[i] = AR + MA + Θ.c
    ε[i] = Y[i]-YHat[i]
    #ε[i+qIndexAdj] = Y[i] - YHat[i]
  end

  return (YHat, ε)
end

#this is from Hamilton
function ARMALogLike(Θ::ARMA, Y::Vector{Float64},YHat::Vector{T} = similar(Y),
    #ε::Vector{T} = zeros(T, length(Y)+max(Θ.q-Θ.p,0))
    ε::Vector{T} = similar(Y)
    )::T where T<:Real

  ARMAFillYHat!(Θ,Y,YHat,ε)

  #from foc of σ²
  ε²::T = sum((abs2).(ε))
  σ²::T = (ε²)/(length(Y)-Θ.p)

  #return -(length(Y)-Θ.p)/2*(log(2π)+log(Θ.σ²))-1/(2*Θ.σ²)*sum(abs2,ε[(Θ.p+1+max(Θ.q-Θ.p,0)):end])
  return -(length(Y)-Θ.p)/2*(log(2π)+log(σ²)) - ε²/(2*σ²)
end

function ARMAUpdate!(Θ::ARMA, θ::Vector{T}, suppressC::Bool)::ARMA where T<:Real
  if suppressC
    Θ.ϕ .= θ[1:(Θ.p)]
    Θ.ϑ .= θ[(Θ.p+1):end]
  else
    Θ.c = θ[1]
    Θ.ϕ .= θ[2:(Θ.p+1)]
    Θ.ϑ .= θ[(Θ.p+2):end]
  end

  return Θ
end

#holds estimate of an ARMA model
struct ARMAEstimate
  Θ::ARMA
  YHat::Vector{Float64}
  ε::Vector{Float64}
  cError::Float64
  ϕErrors::Vector{Float64}
  ϑErrors::Vector{Float64}
  σ²Error::Float64
end

#Algorithms: :LD_LBFGS , :LN_COBYLA, :LN_BOBYQA
#  :LD_SLSQPx , :LD_TNEWTON (fails), :LD_VAR2/:LD_VAR1 ,
#  :LD_MMA , :GN_ESCH, :GN_DIRECT , GN_DIRECT_L ,
# GN_ISRES, :LN_NELDERMEAD, :LN_SBPLX, :LD_LBFGS* (818) , :GD_STOGO
#GN_CRS2_LM
#estimates an ARMA model from the data
function ARMAEstimate(Θ::ARMA, Y::Vector{Float64};
    alg::Symbol = :LD_SLSQP, maxTime::Float64=MAX_TIME, fTolRel::Float64 = F_TOL_REL,
    fTolAbs::Float64 = F_TOL_ABS, bound::Float64 = OPT_BOUND, boundσ::Float64 = OPT_BOUND/10.,
    verbose::Bool = false, suppressC::Bool = false)::ARMAEstimate

  params::Vector{Float64} = suppressC ? [Θ.ϕ; Θ.ϑ] : [Θ.c; Θ.ϕ; Θ.ϑ]
  K::Int = length(params)

  #heuristics for the bound
  lowerBounds::Vector{Float64} = ones(K) .* -1. .* bound
  upperBounds::Vector{Float64} = ones(K) .* bound

  #safety check on parameters
  params .= max.(lowerBounds, (min.(upperBounds, params)))

  #preallocate by calling the fill YHat function, which will allocate the space
  YHat::Vector{Real}, ε::Vector{Real} = ARMAFillYHat!(Θ, Y)

  #use this for finite differencing
  Δz::Vector{Float64} = similar(params)

  #need to do this for autodifferentation
  generalizeARMA!(Θ)

  ##############objective functions
  function singleObj!(θ::Vector{T})::T where T<: Real
    Θ = ARMAUpdate!(Θ,θ, suppressC)
    return ARMALogLike(Θ,Y, YHat, ε)
  end

  function update∇!(θ::Vector{Float64}, ∇::Vector{Float64})::Nothing
    if length(∇) > 0
      ∇ .= ForwardDiff.gradient(singleObj!, θ) #∇Obj(θ)
    end
    return nothing
  end

  function obj!(θ::Vector{Float64}, ∇::Vector{Float64})::Float64
    objVal::Float64 = singleObj!(θ)
    update∇!(θ, ∇)

    return objVal
  end
  ###################optimization
  #println("testObj: $(obj!(params, similar(params)))")

  if verbose
    println("\nBeginning Optimization for ARMA($(Θ.p),$(Θ.q))")
    println("testObj: $(obj!(params, similar(params)))")
  end

  opt::NLopt.Opt = Opt(alg,K)

  max_objective!(opt, obj!)
  lower_bounds!(opt,lowerBounds)
  upper_bounds!(opt,upperBounds)
  #ftol_rel!(opt, fTolRel)
  #ftol_abs!(opt, fTolAbs)
  maxtime!(opt, maxTime)

  #run the actual optimziation
  (optVal,optθ,retCode) = optimize(opt, params)

  if verbose
    println("Objective: $optVal, retCode: $retCode")
    println("Parameters: $optθ")
  end
  #now do the standard errors
  #hessθ::Function = (θ::Vector{Float64})->ForwardDiff.hessian(singleObj!, θ)
  Σ::Matrix{Float64} = (-1 .* ForwardDiff.hessian(singleObj!, optθ)) \ eye(length(optθ)) #invert the hessian
  σ::Vector{Float64} = diag(Σ)
  if minimum(σ) < 0.
    println("WARNING: Negative Standard Variance Estimate, forced to 0")
  end
  σ .= ((x::Float64)->max(x,0.)^0.5).(σ)
  #println("σ2: $σ")

  if verbose

    Δh::Float64 = 10. ^ -9.
    ΔhMat::Matrix{Float64} = eye(length(optθ)) .* Δh
    function finite∇(θ::Vector{Float64}, ∇::Vector{Float64}=similar(θ), origObj::Float64 = singleObj!(θ))
      for i::Int=1:length(θ)
        ∇[i] = (origObj - singleObj!(θ .+ view(ΔhMat, :, i)))/(Δh)
      end

      return ∇
    end

    println("Analytical Gradient: $(ForwardDiff.gradient(singleObj!, optθ))")
    println("finite diff gradient: $(finite∇(optθ))")
  end

  narrowARMA!(Θ)
  ARMAUpdate!(Θ,optθ, suppressC)
  ARMAFillYHat!(Θ, Y, YHat, ε)
  Θ.σ² = sum((abs2).(Y .- YHat))/(length(Y)-Θ.p)
  Θ.dist = Normal(Θ.c, Θ.σ²^0.5)
  ΘEst::ARMAEstimate = suppressC ? ARMAEstimate(Θ, YHat, ε, 0.0, σ[1:Θ.p], σ[(Θ.p+1):end], 0.0) : (
    ARMAEstimate(Θ, YHat, ε, σ[1], σ[2:(Θ.p+1)], σ[(Θ.p+2):end], 0.0))

  return ΘEst
end

#quickly calculates the RMS of two vectors
RMS(v1::Vector{Float64}, v2::Vector{Float64})::Float64 =
  √(abs2(v1 .- v2)/length(v1))

#Based on StatsBase version, but handles missing in a different way
#Takes the prop percentile value, and replaces all lower values with that value
#Then takes the (1-prop) percentile value, and replaces all higher values with that value
#All calculations skip missing
function winsorize!(target::AbstractVector{T};
    prop::Float64 = 0.0, sorted::Bool = false)::Nothing where T <: Any

  local vcomplete::SubArray{T,1}

  if Missing <: T
    vcomplete =  view(target, (!ismissing).(target)) #don't trust how quantile handles missing
  else
    vcomplete=target
  end

  local minval::Float64 = quantile(vcomplete, min(prop,1-prop), sorted=sorted)
  local maxval::Float64 = quantile(vcomplete, max(prop,1-prop), sorted=sorted)

  target .= (f::MFloat64 -> max(minval, min(f,maxval))).(target)
  return nothing
end

#lags a variable in the time series
#assumes the order field is intended to be in ascending order
function lagwithin!(df::DataFrame,
  targets::Vector{Symbol},
  groups::Vector{Symbol},
  period::Symbol;
  lags::Int=1,
  laggednames::Vector{Symbol} = (lags ≠ 1 ?
    (s::Symbol->Symbol(:L, lags, s)).(targets) : (s::Symbol->Symbol(:L, s)).(targets)),
  sorted::Bool=false)::Nothing

  T::Type = eltype(df[!,period]) #the type of the period column

  if !sorted
    sort!(df, [groups; period])
  end

  Ntargets::Int = length(targets)   #pre-allocate the space
  for t ∈ 1:Ntargets
    df[!,laggednames[t]] =
      Vector{Union{eltype(df[!,targets[t]]), Missing}}(undef, size(df, 1))
  end

  for subdf ∈ groupby(df, groups)
    Nsub::Int = size(subdf, 1)

    #only run these routines if we need to
    if Nsub > lags
      for i::Int ∈ (lags+1):Nsub #iterate for all values
        for t::Int ∈ 1:Ntargets
          cur::T = subdf[i, period]

          #first see if we can lag the easy way
          if subdf[i-lags, period] == cur - lags
            subdf[i,laggednames[t]] = subdf[i-lags, targets[t]]
          elseif lags ≠ 1 #now check the hard way (finding the lagged year), pointless if lags==1
            location::NInt = findfirst(isequal(cur-lags), subdf[1:(i-1),period])

            #record the lagged value if it is available
            if !isnothing(location)
              subdf[i,laggednames[t]] = subdf[location, targets[t]]
            end
          end
        end #targets for loop
      end #periods for loop
    end
  end

  return nothing
end

#helper method to handle the case of a single target and/or group
lagwithin!(df::DataFrame,
  targets::Union{Symbol, Vector{Symbol}},
  groups::Union{Symbol, Vector{Symbol}},
  period::Symbol;
  lags::Int=1,
  laggedname::NSymbol = nothing,
  laggednames::Vector{Symbol} = (lags ≠ 1 ?
    (s::Symbol->Symbol(:L, lags, s)).([targets;]) : (s::Symbol->Symbol(:L, s)).([targets;])),
  sorted::Bool=false)::Nothing = (lagwithin!(df, [targets;], [groups;], period,
      lags=lags, laggednames=[something(laggedname, laggednames);], sorted=sorted))

#creates a differenced column
function differencewithin!(df::DataFrame,
  targets::Vector{Symbol},
  groups::Vector{Symbol},
  period::Symbol;
  differencednames::Vector{Symbol} = (s::Symbol->Symbol(:D, s)).(targets),
  sorted::Bool=false,
  createlag::Bool=true,
  deletelag::Bool=true,
  laggednames::Vector{Symbol} = (deletelag ?
    (s::Symbol->Symbol(:L, s, :_temp)).(targets) : (s->Symbol(:L, s)).(targets))
  )::Nothing

  createlag && lagwithin!(df, targets, groups, period,
    laggednames=laggednames, sorted=sorted)
  for t ∈ 1:length(targets)
    df[!, differencednames[t]] = df[!, targets[t]] .- df[!, laggednames[t]]
    deletelag && select!(df, Not(laggednames[t]))
  end

  return nothing
end

#helper method to handle the case of a single target and/or group
differencewithin!(df::DataFrame,
  targets::Union{Symbol, Vector{Symbol}},
  groups::Union{Symbol, Vector{Symbol}},
  period::Symbol;
  differencedname::NSymbol = nothing,
  differencednames::Vector{Symbol} = (s::Symbol->Symbol(:D, s)).([targets;]),
  sorted::Bool=false,
  createlag::Bool=true,
  deletelag::Bool=true,
  laggedname::NSymbol = nothing,
  laggednames::Vector{Symbol} = (deletelag ?
    (s::Symbol->Symbol(:L, s, :_temp)).([targets;]) : (s->Symbol(:L, s)).([targets;]))
  )::Nothing = differencewithin!(df, [targets;], [groups;], period,
    sorted=sorted,  createlag=createlag, deletelag=deletelag,
    differencednames=[something(differencedname, differencednames);],
    laggednames=[something(laggedname, laggednames);])

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
export sbesselj, shankelh1, diffsbessel, diffbessel
export gaunt_coefficients, kernelN
export associated_legendre_indices, spherical_harmonics_indices, lm_to_spherical_harmonic_index
export spherical_harmonics

"""
    sbesselj(m,x)

Returns the spherical besselj function. The order is 'm' and the argument is 'x'. Note 'x' can be a complex number.
"""
function sbesselj(m::Int,x::Union{T,Complex{T}}) where T<:AbstractFloat
    if (abs(x) > eps(T))
        return sqrt(pi/(2*x)) * besselj(m+1/2,x)
    else
        return (m > 0 ? zero(typeof(x)) : one(typeof(x)))
    end
end

"""
    shankelh1(m,x)

Returns the spherical hankel function of the first kind. The order is 'm' and the argument is 'x'. Note 'x' can be a complex number.
"""
shankelh1(m,x) = sqrt(pi/(2*x)) * hankelh1(m+1/2,x)

"""
    diffsbessel(f::Function,m,x)

Differentiates the spherical bessel function 'f' (for any spherical bessel). The order is 'm' and the argument is 'x'. Note 'x' can be a complex number.
"""
function diffsbessel(f::Function,n,z)
    return f(n-1,z) - (n+1) * f(n,z) / z
end

"""
    diffbessel(f::Function,m,x,n::Int)

Differentiates 'n' times any bessel function 'f' of order 'm' and at the argument 'x'.
"""
function diffbessel(f::Function,m,z,n::Int)
    if n == 0
        return f(m, z)
    elseif n > 0
        n = n - 1
        return 0.5*(diffbessel(f,m-1,z,n) - diffbessel(f,m+1,z,n))
    else
        error("Can not differentiate a negative number of times")
    end
end

"""
    diffhankelh1(m,x,n::Int)

Differentiates 'n' times the hankelh1 function of order 'm' and at the argument 'x'.
"""
function diffhankelh1(m,z,n::Int)
    if n == 0
        return hankelh1(m, z)
    elseif n > 0
        n = n - 1
        return 0.5*(diffhankelh1(m-1,z,n) - diffhankelh1(m+1,z,n))
    else
        error("Can not differentiate a negative number of times")
    end
end

"Derivative of Hankel function of the first kind"
diffhankelh1(n,z) = 0.5*(hankelh1(-1 + n, z) - hankelh1(1 + n, z))

"m-th Derivative of Hankel function of the first kind"
function diffbesselj(n,z,m::Int)
    if m == 0
        return besselj(n, z)
    elseif m > 0
        m = m - 1
        return 0.5*(diffbesselj(n-1,z,m) - diffbesselj(n+1,z,m))
    else
        error("Can not differentiate a negative number of times")
    end
end

"Derivative of Bessel function of first kind"
diffbesselj(n,z) = 0.5*(besselj(-1 + n, z) - besselj(1 + n, z))

function kernelN(n::Int,x::Union{T,Complex{T}},y::Union{T,Complex{T}}; dim = 2) where T<:AbstractFloat
    if dim == 2
        h = hankelh1(n,x); dh = diffhankelh1(n,x)
        j = besselj(n,y);  dj = diffbesselj(n,y)
    else
        h = shankelh1(n,x); dh = diffsbessel(shankelh1,n,x)
        j = sbesselj(n,y);  dj = diffsbessel(sbesselj,n,y)
    end

    return x * dh * j - y * h * dj
end

"""
    gaunt_coefficients(l1,m1,l2,m2,l3,m3)

A version of the Gaunt coefficients which are used to write the product of two spherical harmonics. If Y_{l,m} is a complex spherical harmonic, with the typical phase conventions from quantum mechanics, then:

    gaunt_coefficients(l1,m1,l2,m2,l3,m3) = 4*π*im^{l2+l3-l1} Integral[Y_{l1,m1}*conj(Y_{l2,m2})*conj(Y_{l3,m3})]

where the integral is over the solid angle.

The most standard gaunt coefficients `G(l1,m1;l2,m2;l3)` are related through the identity:

    4pi * G(l1,m1;l2,m2;l3) = im^(l1-l2-l3) * (-1)^m2 * gaunt_coefficients(l1,m1,l2,-m2,l3,m1+m2)

"""
function gaunt_coefficients(T::Type{<:AbstractFloat},l1::Int,m1::Int,l2::Int,m2::Int,l3::Int,m3::Int)
    # note the wigner3j has only one convention, and is highly symmetric.
    return (one(T)*im)^(l2+l3-l1) * (-T(1))^m1 * sqrt(4pi*(2*l1+1)*(2*l2+1)*(2*l3+1)) *
        wigner3j(T,l1,l2,l3,0,0,0) * wigner3j(T,l1,l2,l3,m1,-m2,-m3)
end
gaunt_coefficients(l1::Int,m1::Int,l2::Int,m2::Int,l3::Int,m3::Int) = gaunt_coefficients(Float64,l1,m1,l2,m2,l3,m3)

lm_to_spherical_harmonic_index(l::Int,m::Int)::Int = l^2 + m + l + 1

function spherical_harmonics_indices(l_max::Int)
    ls = [l for l in 0:l_max for m in -l:l]
    ms = [m for l in 0:l_max for m in -l:l]

    return ls, ms
end

function associated_legendre_indices(l_max::Int)
    ls = [l for l in 0:l_max for m in 0:l]
    ms = [m for l in 0:l_max for m in 0:l]

    return ls, ms
end

"""
`spherical_harmonics(l_max::Int, θ::T, φ::T)`

returns a vector of all spherical harmonics with degree `l <= l_max`. The degree and order (indices) of the elements of the vector are given by `spherical_harmonics_indices(l_max::Int)`.

The associated legendre polynomials are taken from the package GSL.jl.
"""
function spherical_harmonics(l_max::Int, θ::T, φ::T) where T <: AbstractFloat

    ls, ms = associated_legendre_indices(l_max)
    Plm_arr = sf_legendre_array(GSL_SF_LEGENDRE_SPHARM, l_max, cos(θ))[1:length(ls)]

    Ylm_vec = Vector{Complex{T}}(undef, (l_max+1)^2)
    Ylm_vec[1] = Plm_arr[1]

    ind1 = 1
    ind2 = 1
    for i = 1:l_max
        inds1 = (ind1+i):(ind1+2i)
        Ylm_vec[(ind2+i):(ind2+3i)] = [reverse((-1).^ms[inds1[2:end]] .* conj(Plm_arr[inds1[2:end]])); Plm_arr[inds1]]
        ind1 += i
        ind2 += 2i
    end

    ls, ms = spherical_harmonics_indices(l_max)
    Ylm_vec = exp.(ms .* (im*φ)) .* Ylm_vec

    return Ylm_vec
end


text:
# Stricker Klaus 2020 - License is MIT: http://julialang.org/license
# ReadWriteDlm2 - rwd2_readutils.jl- https://github.com/strickek/ReadWriteDlm2.jl

"""

    parseothers(y::AbstractString, doparsetime::Bool, doparsecomplex::Bool, doparserational::Bool)

Parse string `y` for `Time`, `Complex` and `Rational` format and if match return the value.
Otherwise return the input string `y`.
"""
function parseothers(y, doparsetime, doparsecomplex, doparserational)

    if doparsetime # parse Time
        mt = match(r"^ *(0?\d|1\d|2[0-3])[:Hh]([0-5]?\d)(:([0-5]?\d)([\.,](\d{1,3})(\d{1,3})?(\d{1,3})?)?)? *$", y)
        if mt != nothing
            h = parse(Int, mt[1])
            mi = parse(Int, mt[2])
            (mt[4] == nothing) ? s = ms = us = ns = 0 :
            s  = parse(Int, lpad(string(mt[4]), 2, string(0))); (mt[6] == nothing) ? ms = us = ns = 0 :
            ms = parse(Int, rpad(string(mt[6]), 3, string(0))); (mt[7] == nothing) ? us = ns = 0 :
            us = parse(Int, rpad(string(mt[7]), 3, string(0))); (mt[8] == nothing) ? ns = 0 :
            ns = parse(Int, rpad(string(mt[8]), 3, string(0)))
            return Dates.Time(h, mi, s, ms, us, ns)
        end
    end

    if doparsecomplex # parse Complex
        mc = match(r"^ *(-?\d+(\.\d+)?([eE]-?\d+)?|(-?\d+)//(\d+)) ?([\+-]) ?(\d+(\.\d+)?([eE]-?\d+)?|(\d+)//(\d+))(\*im|\*i|\*j|im|i|j) *$", y)
        if mc != nothing
            real =
                ((mc[4] != nothing) && (mc[5] != nothing)) ? //(parse(Int, mc[4]), parse(Int, mc[5])) :
                ((mc[2] == nothing) && (mc[3] == nothing)) ? parse(Int, mc[1]) : parse(Float64, mc[1])
            imag =
                ((mc[10] != nothing) && (mc[11] != nothing)) ? //(parse(Int, mc[6]*mc[10]), parse(Int, mc[11])) :
                ((mc[8] == nothing) && (mc[9] == nothing)) ? parse(Int, mc[6]*mc[7]) : parse(Float64, mc[6]*mc[7])
            return complex(real, imag)
        end
    end

    if doparserational # parse Rational
        mr = match(r"^ *(-?\d+)//(-?\d+) *$", y)
        if mr != nothing
            nu = parse(Int, mr[1])
            de = parse(Int, mr[2])
            return //(nu, de)
        end
    end

    return y
end


text:
recurs_sym_type(ex::Any) =
  (ex==None || typeof(ex)==Symbol || length(ex.args)==1) ? eval(ex) : Expr(ex.head, ex.args[1], recurs_sym_type(ex.args[2]))
macro c(ret_type, func, arg_types, lib)
  local _arg_types = Expr(:tuple, [recurs_sym_type(a) for a in arg_types.args]...)
  local _ret_type = recurs_sym_type(ret_type)
  local _args_in = Any[ symbol(string('a',x)) for x in 1:length(_arg_types.args) ]
  local _lib = eval(lib)
  quote
    $(esc(func))($(_args_in...)) = ccall( ($(string(func)), $(Expr(:quote, _lib)) ), $_ret_type, $_arg_types, $(_args_in...) )
  end
end

macro ctypedef(fake_t,real_t)
  real_t = recurs_sym_type(real_t)
  quote
    typealias $fake_t $real_t
  end
end

# header: /usr/local/include/nvector/nvector_parallel.h
@ctypedef realtype Float64
@ctypedef N_Vector_Ops Ptr{:Void}
@ctypedef N_Vector Ptr{:Void}
@ctypedef N_Vector_S Ptr{:N_Vector}
@c N_Vector N_VClone (:N_Vector,) shlib
@c N_Vector N_VCloneEmpty (:N_Vector,) shlib
@c None N_VDestroy (:N_Vector,) shlib
@c None N_VSpace (:N_Vector,Ptr{:Int64},Ptr{:Int64}) shlib
@c Ptr{:realtype} N_VGetArrayPointer (:N_Vector,) shlib
@c None N_VSetArrayPointer (Ptr{:realtype},:N_Vector) shlib
@c None N_VLinearSum (:realtype,:N_Vector,:realtype,:N_Vector,:N_Vector) shlib
@c None N_VConst (:realtype,:N_Vector) shlib
@c None N_VProd (:N_Vector,:N_Vector,:N_Vector) shlib
@c None N_VDiv (:N_Vector,:N_Vector,:N_Vector) shlib
@c None N_VScale (:realtype,:N_Vector,:N_Vector) shlib
@c None N_VAbs (:N_Vector,:N_Vector) shlib
@c None N_VInv (:N_Vector,:N_Vector) shlib
@c None N_VAddConst (:N_Vector,:realtype,:N_Vector) shlib
@c realtype N_VDotProd (:N_Vector,:N_Vector) shlib
@c realtype N_VMaxNorm (:N_Vector,) shlib
@c realtype N_VWrmsNorm (:N_Vector,:N_Vector) shlib
@c realtype N_VWrmsNormMask (:N_Vector,:N_Vector,:N_Vector) shlib
@c realtype N_VMin (:N_Vector,) shlib
@c realtype N_VWL2Norm (:N_Vector,:N_Vector) shlib
@c realtype N_VL1Norm (:N_Vector,) shlib
@c None N_VCompare (:realtype,:N_Vector,:N_Vector) shlib
@c Int32 N_VInvTest (:N_Vector,:N_Vector) shlib
@c Int32 N_VConstrMask (:N_Vector,:N_Vector,:N_Vector) shlib
@c realtype N_VMinQuotient (:N_Vector,:N_Vector) shlib
@c Ptr{:N_Vector} N_VCloneEmptyVectorArray (:Int32,:N_Vector) shlib
@c Ptr{:N_Vector} N_VCloneVectorArray (:Int32,:N_Vector) shlib
@c None N_VDestroyVectorArray (Ptr{:N_Vector},:Int32) shlib
@ctypedef N_VectorContent_Parallel Ptr{:Void}
@c N_Vector N_VNew_Parallel (:Int32,:Int64,:Int64) shlib
@c N_Vector N_VNewEmpty_Parallel (:Int32,:Int64,:Int64) shlib
@c N_Vector N_VMake_Parallel (:Int32,:Int64,:Int64,Ptr{:realtype}) shlib
@c Ptr{:N_Vector} N_VCloneVectorArray_Parallel (:Int32,:N_Vector) shlib
@c Ptr{:N_Vector} N_VCloneVectorArrayEmpty_Parallel (:Int32,:N_Vector) shlib
@c None N_VDestroyVectorArray_Parallel (Ptr{:N_Vector},:Int32) shlib
@c None N_VPrint_Parallel (:N_Vector,) shlib
@c N_Vector N_VCloneEmpty_Parallel (:N_Vector,) shlib
@c N_Vector N_VClone_Parallel (:N_Vector,) shlib
@c None N_VDestroy_Parallel (:N_Vector,) shlib
@c None N_VSpace_Parallel (:N_Vector,Ptr{:Int64},Ptr{:Int64}) shlib
@c Ptr{:realtype} N_VGetArrayPointer_Parallel (:N_Vector,) shlib
@c None N_VSetArrayPointer_Parallel (Ptr{:realtype},:N_Vector) shlib
@c None N_VLinearSum_Parallel (:realtype,:N_Vector,:realtype,:N_Vector,:N_Vector) shlib
@c None N_VConst_Parallel (:realtype,:N_Vector) shlib
@c None N_VProd_Parallel (:N_Vector,:N_Vector,:N_Vector) shlib
@c None N_VDiv_Parallel (:N_Vector,:N_Vector,:N_Vector) shlib
@c None N_VScale_Parallel (:realtype,:N_Vector,:N_Vector) shlib
@c None N_VAbs_Parallel (:N_Vector,:N_Vector) shlib
@c None N_VInv_Parallel (:N_Vector,:N_Vector) shlib
@c None N_VAddConst_Parallel (:N_Vector,:realtype,:N_Vector) shlib
@c realtype N_VDotProd_Parallel (:N_Vector,:N_Vector) shlib
@c realtype N_VMaxNorm_Parallel (:N_Vector,) shlib
@c realtype N_VWrmsNorm_Parallel (:N_Vector,:N_Vector) shlib
@c realtype N_VWrmsNormMask_Parallel (:N_Vector,:N_Vector,:N_Vector) shlib
@c realtype N_VMin_Parallel (:N_Vector,) shlib
@c realtype N_VWL2Norm_Parallel (:N_Vector,:N_Vector) shlib
@c realtype N_VL1Norm_Parallel (:N_Vector,) shlib
@c None N_VCompare_Parallel (:realtype,:N_Vector,:N_Vector) shlib
@c Int32 N_VInvTest_Parallel (:N_Vector,:N_Vector) shlib
@c Int32 N_VConstrMask_Parallel (:N_Vector,:N_Vector,:N_Vector) shlib
@c realtype N_VMinQuotient_Parallel (:N_Vector,:N_Vector) shlib

# header: /usr/local/include/nvector/nvector_serial.h
@ctypedef N_VectorContent_Serial Ptr{:Void}
@c N_Vector N_VNew_Serial (:Int64,) shlib
@c N_Vector N_VNewEmpty_Serial (:Int64,) shlib
@c N_Vector N_VMake_Serial (:Int64,Ptr{:realtype}) shlib
@c Ptr{:N_Vector} N_VCloneVectorArray_Serial (:Int32,:N_Vector) shlib
@c Ptr{:N_Vector} N_VCloneVectorArrayEmpty_Serial (:Int32,:N_Vector) shlib
@c None N_VDestroyVectorArray_Serial (Ptr{:N_Vector},:Int32) shlib
@c None N_VPrint_Serial (:N_Vector,) shlib
@c N_Vector N_VCloneEmpty_Serial (:N_Vector,) shlib
@c N_Vector N_VClone_Serial (:N_Vector,) shlib
@c None N_VDestroy_Serial (:N_Vector,) shlib
@c None N_VSpace_Serial (:N_Vector,Ptr{:Int64},Ptr{:Int64}) shlib
@c Ptr{:realtype} N_VGetArrayPointer_Serial (:N_Vector,) shlib
@c None N_VSetArrayPointer_Serial (Ptr{:realtype},:N_Vector) shlib
@c None N_VLinearSum_Serial (:realtype,:N_Vector,:realtype,:N_Vector,:N_Vector) shlib
@c None N_VConst_Serial (:realtype,:N_Vector) shlib
@c None N_VProd_Serial (:N_Vector,:N_Vector,:N_Vector) shlib
@c None N_VDiv_Serial (:N_Vector,:N_Vector,:N_Vector) shlib
@c None N_VScale_Serial (:realtype,:N_Vector,:N_Vector) shlib
@c None N_VAbs_Serial (:N_Vector,:N_Vector) shlib
@c None N_VInv_Serial (:N_Vector,:N_Vector) shlib
@c None N_VAddConst_Serial (:N_Vector,:realtype,:N_Vector) shlib
@c realtype N_VDotProd_Serial (:N_Vector,:N_Vector) shlib
@c realtype N_VMaxNorm_Serial (:N_Vector,) shlib
@c realtype N_VWrmsNorm_Serial (:N_Vector,:N_Vector) shlib
@c realtype N_VWrmsNormMask_Serial (:N_Vector,:N_Vector,:N_Vector) shlib
@c realtype N_VMin_Serial (:N_Vector,) shlib
@c realtype N_VWL2Norm_Serial (:N_Vector,:N_Vector) shlib
@c realtype N_VL1Norm_Serial (:N_Vector,) shlib
@c None N_VCompare_Serial (:realtype,:N_Vector,:N_Vector) shlib
@c Int32 N_VInvTest_Serial (:N_Vector,:N_Vector) shlib
@c Int32 N_VConstrMask_Serial (:N_Vector,:N_Vector,:N_Vector) shlib
@c realtype N_VMinQuotient_Serial (:N_Vector,:N_Vector) shlib


text:
using ImageClipboard
using Images
using Test

@testset "ImageClipboard.jl" begin

    @testset "RGB image copy & paste" begin
        img = rand(RGB{N0f8}, 100, 200)
        clipboard_img(img)
        img2 = clipboard_img()

        # FIXME, windows..
        if Sys.iswindows()
            @test RGBA.(img) == img2
        else
            @test img == img2
        end
    end

    @testset "RGBA image copy & paste" begin
        img = rand(RGBA{N0f8}, 100, 200)
        clipboard_img(img)
        img2 = clipboard_img()

        # FIXME, windows..
        if !Sys.iswindows()
            @test img == img2
        end
    end

end


text:
#
# The functions in this file manage Jacobian evaluations and
# factorizations and function evaluations. The function evaluation bits
# are used in the Newton-Krylov solvers too.
#

"""
For nsoli I use
PrepareJac!(FPS, FS, x, ItRules)
and for ptcsoli
PrepareJac!(FPS, FS, x, ItRules, dt)

Compute the Jacobian and perform the factorization. If know something
about the Jacobian, you can tell me what factorization to use.

For example, if your Jacobian is spd, fact=cholesky! would work well.

"""
function PrepareJac!(FPS, FS, x, ItRules)
    F! = ItRules.f
    J! = ItRules.fp
    dx = ItRules.dx
    fact = ItRules.fact
    pdata = ItRules.pdata
    EvalJ!(FPS, FS, x, F!, J!, dx, pdata)
    TF = fact(FPS)
    return TF
end

function PrepareJac!(FPS, FS, x, ItRules, dt)
    dt > 0 || error("dt must be > 0 in PTC")
    F! = ItRules.f
    J! = ItRules.fp
    dx = ItRules.dx
    fact = ItRules.fact
    jknowsdt = ItRules.jknowsdt
    pdata = ItRules.pdata
    EvalJ!(FPS, FS, x, F!, J!, dx, dt, pdata, jknowsdt)
    TF = fact(FPS)
    return TF
end

"""
PrepareJac!(fc, fm::Real, x, xm, ItRules, dt=0)
Scalar equations
"""
function PrepareJac!(fps::Real, fc, x, ItRules, dt = 0)
    newjac = 0
    newfun = 0
    fp = ItRules.fp
    f = ItRules.f
    dx = ItRules.dx
    pdata = ItRules.pdata
    solver = ItRules.solver
    df = fpeval_newton(x, f, fc, fp, dx, pdata)
    dt == 0 || (df += 1.0 / dt)
    newjac = newjac + 1
    return df
end

"""
klfact(A)

Returns the default choice for the factorization unless you tell
me to do something else. QR is the default choice for banded because
that works properly with Float32.

"""
function klfact(A::Array{T,2}) where {T<:Real}
    TF = lu!(A)
end

# The default for sparse is lu. lu! for sparse matrices is
# too complicated to put in here. You can use lu! if you
# set fact = nofact and manage the factorization in your Jacobian
# evaluation code. You'll also get to manage the storage. There's
# a project in chapter 2 about that.
#
function klfact(A::SparseMatrixCSC{Float64,Int64})
    TF = lu(A)
end

# The default for banded matrices is qr, because I do not trust
# you to allocate the extra two upper bands so I can use qr!.
# I'm using qr! in the example in Chapter 2. Look at the source
# to see how I did that.
#
function klfact(A::BandedMatrix)
    TF = qr(A)
end

# Default: do nothing.
function klfact(A)
    TF = nofact(A)
end

function nofact(A)
    TF = A
end

"""
EvalF!(F!, FS, x, pdata)

This is a wrapper for the function evaluation that figures out if
you are using precomputed data or not. No reason to get excited
about this.
"""
function EvalF!(F!, FS, x, q::Nothing)
    FS=F!(FS, x)
    return FS
end

function EvalF!(F!, FS, x, pdata)
    FS=F!(FS, x, pdata)
    return FS
end

function EvalF!(F!, FS::Real, x::Real, q::Nothing)
    FS = F!(x)
    return FS
end

function EvalF!(F!, FS::Real, x::Real, pdata)
    FS = F!(x, pdata)
    return FS
end

"""
If you let me handle dt in PTC
JV!(FPS, FS, x, J!, pdata)

If you put the (1/dt) * I in the Jacobian yourself
JV!(FPS, FS, x, dt, J!, pdata)

This is a wrapper for the Jacobian evaluation that figures out if
you are using precomputed data or not. No reason to get excited
about this.
"""
function JV!(FPS, FS, x, J!, pdata)
    J!(FPS, FS, x, pdata)
end

function JV!(FPS, FS, x, dt, J!, pdata)
    J!(FPS, FS, x, dt, pdata)
end

function JV!(FPS, FS, x, dt, J!, q::Nothing)
    J!(FPS, FS, x, dt)
end

function JV!(FPS, FS, x, J!, q::Nothing)
    J!(FPS, FS, x)
end

"""
for Newton
EvalJ!(FPS, FS, x, F!, J!, dx, pdata)

for PTC
EvalJ!(FPS, FS, x, F!, J!, dx, dt, pdata)

evaluates the Jacobian before the factorization in PrepareJac!

"""

function EvalJ!(FPS, FS, x, F!, J!, dx, dt, pdata, jknowsdt)
    #    if J! != diffjac!
    #        JV!(FPS, FS, x, J!, pdata)
    #    else
    #        diffjac!(FPS, FS, F!, x, dx, pdata)
    #    end
    if jknowsdt
        FPS = JV!(FPS, FS, x, dt, J!, pdata)
    else
        EvalJ!(FPS, FS, x, F!, J!, dx, pdata)
        FPS .= FPS + (1.0 / dt) * I
    end
    return FPS
end

function EvalJ!(FPS, FS, x, F!, J!, dx, pdata)
    if J! != diffjac!
        JV!(FPS, FS, x, J!, pdata)
    else
        diffjac!(FPS, FS, F!, x, dx, pdata)
    end
    return FPS
end

"""
   diffjac!(FPS::Array{T,2}, FS, F!, x, dx, pdata) where T <: Real

Computes a finite-difference dense and unstructured Jacobian.
This is not something an user wants to mess with. Look at the
docstrings to nsold to see more details.

Nothing much to see here. Move along.
"""
function diffjac!(FPS::Array{T,2}, FS, F!, x, dx, pdata) where {T<:Real}
    h = dx * norm(x, Inf) + 1.e-8
    n = length(x)
    y = ones(size(x))
    FY = ones(size(x))
    for ic = 1:n
        y .= x
        y[ic] = y[ic] + h
        EvalF!(F!, FY, y, pdata)
        for ir = 1:n
            FPS[ir, ic] = (FY[ir] - FS[ir]) / h
        end
    end
    return FPS
end

"""
UpdateIteration

Take a trial step. Evaluate the function and the residual norm.
"""
function UpdateIteration(xt::Array{T}, x, FS, lambda, step, ItRules) where {T<:Real}
    F! = ItRules.f
    pdata = ItRules.pdata
    copy!(xt, x)
    BLAS.axpy!(lambda, step, xt)
    EvalF!(F!, FS, xt, pdata)
    resnorm = norm(FS)
    return (xt, FS, resnorm)
end

function UpdateIteration(xt::T, xm, ft, lambda, d, ItRules) where {T<:Real}
    f = ItRules.f
    pdata = ItRules.pdata
    xt = xm + lambda * d
    fc = 0.0
    fc = EvalF!(f, fc, xt, pdata)
    #fc = f(xt)
    residc = norm(fc)
    return (xt, fc, residc)
end

"""
fpeval_newton

Evaluates f' by differences or the user's code.

"""
function fpeval_newton(x, f, fc, fp, h, pdata)
    fps = string(fp)
    df = 0.0
    dps = string(difffp)
    if fps == dps
        df = difffp(x, f, fc, h, pdata)
    else
        df = EvalF!(fp, df, x, pdata)
    end
    return df
end

"""
difffp

forward differencing for scalar equations
"""
function difffp(x, f, fc, h, pdata)
    fph = 0.0
    fph = EvalF!(f, fph, x + h, pdata)
    df = (fph - fc) / h
    #    df = (f(x + h) - fc) / h
    return df
end

"""
test_evaljac(ItRules, itc, newiarm, residratio)

Figures out if it's time to reevaluate and refacto the Jacbian in
Newton's method.
"""
function test_evaljac(ItRules, itc, newiarm, residratio)
    solver = ItRules.solver
    sham = ItRules.sham
    resdec = ItRules.resdec
    evaljacit = (itc % sham == 0 || newiarm > 0 || residratio > resdec)
    chordinit = (solver == "chord") && itc == 0
    evaljac = (evaljacit && solver == "newton") || chordinit || solver == "secant"
end