text:
# TODO:
# * add combine(fun, df) for DataFrame with 0 rows

# 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, Symbol}}
# 3) Pair{AbstractVector{Int}, <:Pair{<:Base.Callable, Symbol}}
# 4) Pair{AsTable, <:Pair{<:Base.Callable, Symbol}}

"""
    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. Therefore,
to allow for future extensions, returning `NamedTuple` or `DataFrameRow`
from `fun` is currently disallowed.
"""
struct ByRow{T} <: Function
    fun::T
end

_by_row_helper(x::Any) = x
_by_row_helper(x::Union{NamedTuple, DataFrameRow}) =
    throw(ArgumentError("return value of type $(typeof(x)) " *
                        "is currently not allowed with ByRow."))

(f::ByRow)(cols::AbstractVector...) = _by_row_helper.(f.fun.(cols...))
(f::ByRow)(table::NamedTuple) =
    _by_row_helper.(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) =
    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::Pair{typeof(nrow), Symbol}) =
    length(idx) == 0 ? (Int[] => (() -> 0) => last(sel)) : (1 => length => last(sel))
normalize_selection(idx::AbstractIndex, sel::Pair{typeof(nrow), <:AbstractString}) =
    normalize_selection(idx, first(sel) => Symbol(last(sel)))
normalize_selection(idx::AbstractIndex, sel::typeof(nrow)) =
    normalize_selection(idx, nrow => :nrow)

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

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

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

function normalize_selection(idx::AbstractIndex,
                             sel::Pair{<:Any,<:Pair{<:Base.Callable, Symbol}})
    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 length(c) == 0 && first(last(sel)) isa ByRow
        throw(ArgumentError("at least one column must be passed to a " *
                            "`ByRow` transformation function"))
    end
    return (wanttable ? AsTable(c) : c) => last(sel)
end

normalize_selection(idx::AbstractIndex,
                    sel::Pair{<:Any,<:Pair{<:Base.Callable,<:AbstractString}}) =
    normalize_selection(idx, first(sel) => first(last(sel)) => Symbol(last(last(sel))))

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

function normalize_selection(idx::AbstractIndex,
                             sel::Pair{<:Any, <:Base.Callable})
    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 length(c) == 0 && last(sel) isa ByRow
        throw(ArgumentError("at least one column must be passed to a " *
                            "`ByRow` transformation function"))
    end
    fun = last(sel)
    if length(c) > 3
        newcol = Symbol(join(@views(_names(idx)[c[1:2]]), '_'), "_etc_", funname(fun))
    elseif isempty(c)
        newcol = Symbol(funname(fun))
    else
        newcol = Symbol(join(view(_names(idx), c), '_'), '_', funname(fun))
    end
    return (wanttable ? AsTable(c) : c) => fun => newcol
end

function select_transform!(nc::Pair{<:Union{Int, AbstractVector{Int}, AsTable},
                                    <:Pair{<:Base.Callable, Symbol}},
                           df::AbstractDataFrame, newdf::DataFrame,
                           transformed_cols::Dict{Symbol, Any}, copycols::Bool,
                           allow_resizing_newdf::Ref{Bool})
    col_idx, (fun, newname) = nc
    # 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
    @assert !hasproperty(newdf, newname)
    cdf = eachcol(df)
    if col_idx isa Int
        res = fun(df[!, col_idx])
    elseif col_idx isa AsTable
        res = fun(Tables.columntable(select(df, col_idx.cols, copycols=false)))
    else
        # it should be fast enough here as we do not expect to do it millions of times
        @assert col_idx isa AbstractVector{Int}
        res = fun(map(c -> cdf[c], col_idx)...)
    end
    if res isa Union{AbstractDataFrame, NamedTuple, DataFrameRow, AbstractMatrix}
        throw(ArgumentError("return value from function $fun " *
                            "of type $(typeof(res)) is currently not allowed."))
    end
    if res isa AbstractVector
        # 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, length(res)), 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 && length(res) != nrow(df)
            throw(ArgumentError("length $(length(res)) of vector returned from " *
                                "function $fun is different from number of rows " *
                                "$(nrow(df)) of the source data frame."))
        end
        allow_resizing_newdf[] = false
        respar = parent(res)
        parent_cols = col_idx isa AsTable ? col_idx.cols : col_idx
        if copycols && !(fun isa ByRow) &&
            (res isa SubArray || any(i -> respar === parent(cdf[i]), parent_cols))
            newdf[!, newname] = copy(res)
        else
            newdf[!, newname] = res
        end
    else
        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
    # mark that column transformation was applied
    # nothing is not possible otherwise as a value in this dict
    transformed_cols[newname] = nothing
end

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

    * Any index that is allowed for column indexing
      ($COLUMNINDEX_STR; $MULTICOLUMNINDEX_STR).
    * Column transformation operations using the `Pair` notation that is
      described below and vectors 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, and `fun` 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.

    If `fun` returns a value of type other than `AbstractVector` then it will be
    broadcasted into 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
    currently these return types are not allowed. As a particular rule, values
    wrapped in a `Ref` or a `0`-dimensional `AbstractArray` are unwrapped and
    then broadcasted.

    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 not allowed for
    `old_column` to select an empty set of columns nor for `fun` to return a
    `NamedTuple` or a `DataFrameRow`.

    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)`. 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`.

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

    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...)

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) .=> sum);

julia> df
3×2 DataFrame
│ Row │ a_sum │ b_sum │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 6     │ 15    │
│ 2   │ 6     │ 15    │
│ 3   │ 6     │ 15    │

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

julia> using Statistics

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

"""
select!(df::DataFrame, args...) =
    _replace_columns!(df, select(df, args..., copycols=false))

"""
    transform!(df::DataFrame, args...)

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, args...) = select!(df, :, args...)

"""
    select(df::AbstractDataFrame, args...; copycols::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, :b)
3×1 DataFrame
│ Row │ b     │
│     │ Int64 │
├─────┼───────┤
│ 1   │ 4     │
│ 2   │ 5     │
│ 3   │ 6     │

julia> select(df, Not(:b)) # drop column :b from df
3×1 DataFrame
│ Row │ a     │
│     │ Int64 │
├─────┼───────┤
│ 1   │ 1     │
│ 2   │ 2     │
│ 3   │ 3     │

julia> select(df, :a => :c, :b)
3×2 DataFrame
│ Row │ c     │ b     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 4     │
│ 2   │ 2     │ 5     │
│ 3   │ 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) .=> sum)
3×2 DataFrame
│ Row │ a_sum │ b_sum │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 6     │ 15    │
│ 2   │ 6     │ 15    │
│ 3   │ 6     │ 15    │

julia> select(df, names(df) .=> sum .=> [:A, :B])
3×2 DataFrame
│ Row │ A     │ B     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 6     │ 15    │
│ 2   │ 6     │ 15    │
│ 3   │ 6     │ 15    │

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

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

"""
    transform(df::AbstractDataFrame, args...; copycols::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, args...; copycols::Bool=true) =
    select(df, :, args..., copycols=copycols)

"""
    combine(df::AbstractDataFrame, args...)
    combine(arg, df::AbstractDataFrame)

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)
1×2 DataFrame
│ Row │ a_sum │ nrow  │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 6     │ 3     │
```
"""
combine(df::AbstractDataFrame, args...) =
    manipulate(df, args..., copycols=true, keeprows=false)

function combine(arg, df::AbstractDataFrame)
    if nrow(df) == 0
        throw(ArgumentError("calling combine on a data frame with zero rows" *
                            " with transformation as a first argument is " *
                            "currently not supported"))
    end
    return combine(arg, groupby(df, Symbol[]))
end

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

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

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

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

function _manipulate(df::AbstractDataFrame, 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     │
    #
    # we compute column :a immediately when we process `:` although it is specified
    # later by `:a=>sin=>:a` because we know from `transformed_cols` variable that
    # it will be computed later via a transformation
    transformed_cols = Dict{Symbol, Any}()
    for nc in normalized_cs
        if nc isa Pair
            newname = last(last(nc))
            @assert newname isa Symbol
            if haskey(transformed_cols, newname)
                throw(ArgumentError("duplicate target column name $newname passed"))
            end
            transformed_cols[newname] = nc
        end
    end
    # 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}
            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)
                    if haskey(transformed_cols, newname)
                        # if newdf does not have a column newname
                        # but a column transformation was requested for this column
                        # then apply the transformation immediately
                        # in such a case nct may not be nothing, as if it were
                        # nothing then newname should be preasent in newdf already
                        nct = transformed_cols[newname]
                        @assert nct !== nothing
                        select_transform!(nct, df, newdf, transformed_cols, copycols,
                                          allow_resizing_newdf)
                    else
                        # 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
            end
        else
            # nc is normalized so it has a form src_cols => fun => Symbol
            newname = last(last(nc))
            if hasproperty(newdf, newname)
                # it is possible that the transformation has already been applied
                # via multiple column selection, like in select(df, :, :x1 => :y1)
                # but then transformed_cols[newname] must be nothing
                @assert transformed_cols[newname] === nothing
            else
                select_transform!(nc, df, newdf, transformed_cols, copycols,
                                  allow_resizing_newdf)
            end
        end
    end
    return newdf
end

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

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

function manipulate(dfv::SubDataFrame, args...; copycols::Bool, keeprows::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) 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)
                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:
function validate_buoyancy(buoyancy, tracers)
    req_tracers = required_tracers(buoyancy)

    all(tracer ∈ tracers for tracer in req_tracers) ||
        error("$(req_tracers) must be among the list of tracers to use $(typeof(buoyancy).name.wrapper)")

    return nothing
end

validate_buoyancy(::Nothing, tracers) = nothing

ext: jl


text:
module HexGeometry

# examples of how to use polygons defined in GeoJSON filled with hexes for indexing and searching
using Plots
using GeoInterface
using LibGEOS
using H3
using H3.Lib
using H3.Lib: H3Index, GeoCoord, Geofence, GeoPolygon
import H3.API.geoToH3
using Geodesy
using Statistics

import Plots.plot, Plots.plot!

struct Hex
    center::Vector{Float64}
    epsilon::Float64
    shape::LibGEOS.Polygon
end

"""
    Hex(center, size; angle=0, adj=(0, 0))

    Creates a regular hexagon as a LibGEOS polygon. In the simplest case
    center is at the center of the hexagon and size gives the length
    of any side. You can rotate the hexagon with the optional keyword
    argument angle or you can adjust what point in the hexagon that the
    center actually represents. Use adj=(-1,0) to specify the left side,
    adj=(0,-1) to specify the center of the bottom side and so on.

    # Examples
    Honey comb
    ```julia
    julia> h1 = Hex((0,0), 1)
    julia> h2 = Hex((0,sqrt(3)), 1)
    julia> h3 = Hex((0,-sqrt(3)), 1)
    julia> h4 = Hex((0.5,sqrt(3)/2), 1, adj=(-1,0))
    julia> h5 = Hex((-0.5,sqrt(3)/2), 1, adj=(-1,0))
    julia> h6 = Hex((0.5,-sqrt(3)/2), 1, adj=(-1,0))
    julia> h7 = Hex((-0.5,-sqrt(3)/2), 1, adj=(-1,0))
    ```
"""
function Hex(center, epsilon; angle=0, adj=[0,0])
    base = [0 0;0.5 sqrt(3)/2; 1.5 sqrt(3)/2; 2 0; 1.5 -sqrt(3)/2; 0.5 -sqrt(3)/2; 0 0] .- [1 0]
    center = center * epsilon
    rotation = [cos(angle) sin(angle); -sin(angle) cos(angle)]
    v = (base .+ + ([2, sqrt(3)] .* adj)') * rotation .+ center'
    z = [ [v[i,1],v[i,2]] for i in 1:size(v)[1] ]
    return Hex(center, epsilon, LibGEOS.Polygon([z]))
end

function kring()
    return reduce(LibGEOS.union, [
        Hex([0,0], 1).shape,
        Hex([0,0],1,adj=[0.75,0.5]).shape,
        Hex([0,0],1,adj=[0.75,-0.5]).shape,
        Hex([0,0],1,adj=[-0.75,-0.5]).shape,
        Hex([0,0],1,adj=[-0.75,0.5]).shape,
        Hex([0,0],1,adj=[-0,1]).shape,
        Hex([0,0],1,adj=[-0,-1]).shape
        ])
end

function hexSafe(center, radius, ring)
    circle = buffer(LibGEOS.Point(center...), radius)
    r = difference(circle, ring)
    return area(r) / area(circle)
end

function hexSafe(radius)
    core = Hex([0,1], 1).shape
    ring = kring()
    r = []
    centers = (rand(10000, 2) .- [0.5 0.5]) .* [2 sqrt(3)]
    centers = centers[ [contains(core, LibGEOS.Point(centers[i,:])) for i in 1:10000], : ]
    return mean([hexSafe(centers[i, :], radius, ring) for i in 1:size(centers)[2]])
end

function plot(h::Hex)
    plot(h.shape)
end

function plot!(h::Hex)
    plot!(h.shape)
end

"""
    geoToH3(p, resolution)

Converts a point p defined in terms of latitude and longitude in degrees into an H3Index.

# Example

geoToH3(LibGEOS.Point(42, -110), 6) => 0x08626b3cafffffff
"""
function geoToH3(p::AbstractPoint, resolution::Int) :: H3Index
    geoToH3(coordinates(p)..., resolution)
end

"""
# geoToH3(latitude, longitude, resolution)

Converts a point defined in terms of latitude and longitude in degrees into an H3Index.

# Example

julia> geoToH3(42, -110, 6)
0x08626b3cafffffff

julia> poly = LibGEOS.Polygon([[[0.0,0.0],[1,1],[2,0],[0,0]]])
julia> (x->geoToH3(x..., 5)).(coordinates(boundary(poly)))
4-element Vector{UInt64}:
 0x085754e67fffffff
 0x0857541affffffff
 0x0857542b7fffffff
 0x085754e67fffffff
"""
function geoToH3(latitude::Number, longitude::Number, resolution::Int) :: H3Index
    geoToH3(H3.API.GeoCoord(deg2rad(latitude), deg2rad(longitude)), resolution)
end

"""
# geoToH3([latitude, longitude], resolution)

Converts a point defined in terms of latitude and longitude in degrees into an H3Index.

# Example

julia> geoToH3(42, -110, 6)
0x08626b3cafffffff

julia> poly = LibGEOS.Polygon([[[0.0,0.0],[1,1],[2,0],[0,0]]])
julia> (x->geoToH3(x..., 5)).(coordinates(boundary(poly)))
4-element Vector{UInt64}:
 0x085754e67fffffff
 0x0857541affffffff
 0x0857542b7fffffff
 0x085754e67fffffff
"""
function geoToH3(point::Vector{Float64}, resolution::Int) :: H3Index
    geoToH3(point..., resolution)
end

"""
# h3ToPolygon(id)

Converts a single H3Index into the corresponding GeoPolygon
"""
function h3ToPolygon(id::H3Index)
    px = map(p -> [rad2deg(p.lat), rad2deg(p.lon)], H3.API.h3ToGeoBoundary(id))
    return closedPolygon(px)
end

"""
# closedPolygon(coordinates)

Add the first point in a vector of coordinates to the end and then creates
a polygon fromt the result.
"""
function closedPolygon(coordinates::Vector{Vector{Float64}}) :: AbstractPolygon
    push!(coordinates, coordinates[1])
    return LibGEOS.Polygon([coordinates])
end

"""
# h3ToPolygon(Vector{H3Index})

Converts a vector of hexagons (represented as H3Index values) into a polygon which
is the union of all of the hexagons.
"""
function h3ToPolygon(ids::Vector{H3Index})
    return reduce(LibGEOS.union, h3ToPolygon.(ids))
end

"""
# polyfill(p::AbstractPolygon, resolution::Int; cover::Float64)

Covers the polygon p (which is expressed in lat/long) with hexagons and returns
the list of H3Indexes for each. Note that this is only an approximate covering.
Some hexagons will extend outside the polygon being covered and not all of the
polygon will be covered.

To avoid undercoverage this function expands poly by an amount determined by the
cover parameter. The default value will result in a cover that rarely (if ever)
has any under-coverage. Increasing the setting to 2 will prevent any under-coverage
but will result in lots of overspray. Setting cover=0 will avoid all expansion, but will
almost always result in significant under-coverage as well. At about cover=1.1, the
undercoverage should be very small (<< 0.1%), but will still happen on occasion.

# Example
julia> p = LibGEOS.Polygon([[[40.0, -110.0], [41,-110], [41.4,-109], [40,-110]]])
julia> plot(h3ToPolygon(polyfill(p, 2)))
julia> plot!(h3ToPolygon(polyfill(p, 3)))
julia> plot!(h3ToPolygon(polyfill(p, 4)))
julia> plot!(h3ToPolygon(polyfill(p, 5)))
julia> plot!(p)
## Highlight under-coverage
julia> px = h3ToPolygon(polyfill(p, 4, cover=0)
julia> plot(p)
julia> plot(difference(p,px))
julia> plot(difference(px,p))
"""
function polyfill(poly::AbstractPolygon, resolution::Int; cover=1.2)::Vector{H3Index}
    cover = max(0.0, cover)
    if cover > 0.0
        # we expand by an amount related to the size of a hex edge
        fudge = cover * Lib.edgeLengthM(resolution)

        # this projection puts the first vertex of the polygon at the origin
        # and is expressed in meters [east, north, up] against the wgs84 ellipsoid
        origin = coordinates(boundary(poly))[1]
        tx = ENUfromLLA(LLA(origin...,0), wgs84)

        # expand our polygon projected form. We also simplify to avoid large number of vertices
        projected = map(x->tx(LLA(x..., 0))[1:2], coordinates(boundary(poly)))
        expanded = simplify(buffer(LibGEOS.Polygon([projected]), fudge), fudge/5)

        # project back
        coords = map(lla-> [lla.lat, lla.lon], map(x->inv(tx)(ENU(x..., 0)), coordinates(boundary(expanded))))
        poly = LibGEOS.Polygon([coords])
    end
    vertices = map(x->GeoCoord(map(deg2rad, x)...), coordinates(boundary(poly)))
    fence = Geofence(length(vertices), pointer(vertices))
    h3poly = GeoPolygon(fence, 0, C_NULL)
    hexagons = Lib.polyfill(Ref(h3poly), resolution)
    return filter(!iszero, hexagons)
end

end # module

max_stars_repo_name: billmetangmo/weather-server

max_issues_repo_name: billmetangmo/weather-server

max_forks_repo_name: billmetangmo/weather-server


text:
"""
    module TightlyBound

Main Module
"""
module TightlyBound

include("GlobalUnits.jl")

include("SetDir.jl")

include("Utility.jl")
include("BandTools.jl")
include("Atomic.jl")
include("Atomdata.jl")
include("Crystal.jl")

using Printf
using .CrystalMod:crystal
using .CrystalMod:makecrys

export makecrys

include("Ewald.jl")
include("DFToutMod.jl")
using .DFToutMod:dftout
using .DFToutMod:makedftout
export dftout
export makedftout

include("TB.jl")

using .TB:tb_crys
using .TB:tb_crys_kspace

using .TB:Hk
using .TB:calc_bands
using .TB:read_tb_crys
using .TB:write_tb_crys

include("BandStruct.jl")

using .BandStruct:plot_compare_tb
using .BandStruct:plot_bandstr
using .BandStruct:plot_compare_dft

include("DOS.jl")
using .DOS:dos

export Hk
export calc_bands
export plot_compare_tb
export plot_bandstr
export plot_compare_dft
export read_tb_crys
export dos
export write_tb_crys

include("RunDFT.jl")

export set_units

#include("RunWannier90.jl")
include("AtomicProj.jl")
include("CalcTB_laguerre.jl")
using .CalcTB:calc_tb_fast
export calc_tb_fast

include("SCF.jl")

include("FitTB_laguerre.jl")
include("Force_Stress.jl")

include("ManageDatabase.jl")

include("MyOptim.jl")

include("Relax.jl")
using .CrystalMod:print_with_force_stress

export scf_energy
export scf_energy_force_stress
export relax_structure

"""
    function set_units(;energy=missing, length=missing, both=missing)

Set global units for energy/length. Run with no arguments to check/return current units.

- Default units are `"eV"` and `"Angstrom"` (or `"Å"` or `"Ang"` ).
- Choose atomic units by setting `energy="Ryd"` and `length="Bohr"`.
- Set both at the same time with `both="atomic"` or `both="eVAng"`
- Internally, all units are atomic. Only main public facing functions actually change units.
"""
function set_units(;energy=missing, length=missing, both=missing)

    if !ismissing(both)
        both = String(both)
        if both == "atomic" || both == "Atomic" || both == "au" || both == "a.u."
            energy="Ryd."
            length="Bohr"
        elseif both == "eVang" || both == "AngeV" || both == "eVAng" || both == "evang"
            energy="eV"
            length="Å"
        else
            println("I don't understand both : ", both)
        end

    end

    if !ismissing(length)
        length = String(length)
        if length == "A" || length == "Ang" || length == "Ang." || length == "Angstrom" || length == "a" || length == "ang" || length == "ang." || length == "angstrom" || length == "Å"

            global global_length_units="Å"
        elseif length == "Bohr" || length == "bohr" || length == "au" || length == "a.u." || length == "atomic"
            global global_length_units="Bohr"
        else
            println("I don't understand length : ", length)
        end
    end

    if !ismissing(energy)
        energy = String(energy)
        if energy == "eV" || energy == "ev" || energy == "EV" || energy == "electronvolts"
            global global_energy_units="eV"
        elseif energy == "Rydberg" || energy == "Ryd." || energy == "Ryd" || energy == "atomic" || energy == "au" || energy == "a.u." || energy == "rydberg" || energy == "ryd" || energy == "ryd."
            global global_energy_units="Ryd."
        else
            println("I don't understand energy : ", energy)
        end

    end

    println("Units are now $global_energy_units and $global_length_units ")

    return deepcopy(global_energy_units), deepcopy(global_length_units)

end

"""
    relax_structure(c::crystal; mode="vc-relax")

Find the lowest energy atomic configuration of crystal `c`.

# Arguments
- `c::crystal`: the structure to relax, only required argument
- `mode="vc-relax"`: Default (variable-cell relax) will relax structure and cell, anything else will relax structure only.
- `database=missing`: coefficent database, default is to use the pre-fit pbesol database
- `smearing=0.01`: smearing temperature (ryd), default = 0.01
- `grid=missing`: k-point grid, e.g. [10,10,10], default chosen automatically
- `nsteps=100`: maximum iterations
- `update_grid=true`: update automatic k-point grid during relaxation
- `conv_thr = 2e-4 `: Convergence threshold for combo of forces, positions, energy
"""
function relax_structure(c::crystal; database=missing, smearing = 0.01, grid = missing, mode="vc-relax", nsteps=50, update_grid=true, conv_thr = 2e-4)

    if ismissing(database)
        ManageDatabase.prepare_database(c)
        database = ManageDatabase.database_cached
    end
    if !ismissing(grid)
        update_grid = false
    end

    cfinal, tbc, energy, force, stress = Relax.relax_structure(c, database, smearing=smearing, grid=grid, mode=mode, nsteps=nsteps, update_grid=update_grid, conv_thr=conv_thr)

    println("Relax done")

#    println("Calculate final energy")
#
#    energy_tot, tbc, conv_flag = scf_energy(cfinal; database=database, smearing=0.01, grid = missing)
#    energy_tot, f_cart, stress = scf_energy_force_stress(tbc, database=database, smearing=smearing, grid=grid)
#
#    println("done with all relax")
#    println()

    energy = convert_energy(energy)
    println("Final Energy $energy ")

    print_with_force_stress(cfinal, force, stress)
    println()

    energy = convert_energy(energy)
    force = convert_force(force)
    stress = convert_stress(stress)

    println()

    return cfinal, tbc, energy, force, stress

end

"""
    scf_energy_force_stress(c::crystal; database = missing, smearing = 0.01, grid = missing)

Calculate energy, force, and stress for a crystal.

`return energy, force, stress, tight_binding_crystal_object`

# Arguments
- `c::crystal`: the structure to calculate on. Only required argument.
- `database=missing`: Source of coeficients. Will be loaded from pre-fit coefficients if missing.
- `smearing=0.01`: Gaussian smearing temperature, in Ryd. Usually can leave as default.
- `grid=missing`: k-point grid, e.g. [10,10,10], default chosen automatically
"""
function scf_energy_force_stress(c::crystal; database = missing, smearing = 0.01, grid = missing)

    energy_tot, tbc, conv_flag = scf_energy(c; database=database, smearing=smearing, grid = grid)

    if ismissing(database)
        database = ManageDatabase.database_cached
    end

    println()
    println("Calculate Force, Stress")

    energy_tot, f_cart, stress = Force_Stress.get_energy_force_stress_fft(tbc, database, do_scf=false, smearing=smearing, grid=grid)

    println("done")
    println("----")

    energy_tot = convert_energy(energy_tot)
    f_cart = convert_force(f_cart)
    stress = convert_stress(stress)

    return energy_tot, f_cart, stress, tbc

end

"""
    scf_energy_force_stress(tbc::tb_crys; database = missing, smearing = 0.01, grid = missing)

Calculate energy, force, and stress for a tight binding crystal
object. This allows the calculation to run without re-doing the SCF
calculation. Assumes SCF already done!

returns energy, force, stress, tight_binding_crystal_object

"""
function scf_energy_force_stress(tbc::tb_crys; database = missing, smearing = 0.01, grid = missing, do_scf=false)

    if ismissing(database)
        ManageDatabase.prepare_database(tbc.crys)
        database = ManageDatabase.database_cached
    end

    println()
    println("Calculate Force, Stress (no scf)")

    energy_tot, f_cart, stress = Force_Stress.get_energy_force_stress_fft(tbc, database, do_scf=false, smearing=smearing, grid=grid)

    println("done")
    println("----")

    energy_tot = convert_energy(energy_tot)
    f_cart = convert_force(f_cart)
    stress = convert_stress(stress)

    return energy_tot, f_cart, stress, tbc

end

"""
    scf_energy(c::crystal)

Calculate energy, force, and stress for a crystal. Does self-consistent-field (SCF) calculation if using self-consistent electrostatics.

returns energy, tight-binding-crystal-object, error-flag

# Arguments
- `c::crystal`: the structure to calculate on. Only required argument.
- `database=missing`: Source of coeficients. Will be loaded from pre-fit coefficients if missing.
- `smearing=0.01`: Gaussian smearing temperature, in Ryd. Usually can leave as default.
- `grid=missing`: k-point grid, e.g. [10,10,10], default chosen automatically
- `conv_thr = 1e-5`: SCF convergence threshold (Ryd).
- `iter = 75`: number of iterations before switch to more conservative settings.
- `mix = -1.0`: initial mixing. -1.0 means use default mixing. Will automagically adjust mixing if SCF is failing to converge.
- `mixing_mode =:pulay`: default is Pulay mixing (DIIS). Other option is :simple, for simple linear mixing of old and new electron-density. Will automatically switch to simple if Pulay fails.
"""
function scf_energy(c::crystal; database = missing, smearing=0.01, grid = missing, conv_thr = 1e-5, iters = 75, mix = -1.0, mixing_mode=:pulay)
    println()
    println("Begin scf_energy-------------")
    println()
    if ismissing(database)
        println("Load TB parameters from file")
        ManageDatabase.prepare_database(c)
        database = ManageDatabase.database_cached
        println()
    end

    energy_tot, efermi, e_den, dq, VECTS, VALS, error_flag, tbc = SCF.scf_energy(c, database, smearing=smearing, grid = grid, conv_thr = conv_thr, iters = iters, mix = mix,  mixing_mode=mixing_mode)

    conv_flag = !error_flag
    if tbc.within_fit == false
        conv_flag = false
    end
    if conv_flag == true
        println("scf_energy success, done")
    else
        println("scf_energy error detected somewhere!!!!!!!!!!, done")
    end
    println()

    energy_tot = convert_energy(energy_tot)

    return energy_tot, tbc, conv_flag

end

"""
    scf_energy(d::dftout)

    SCF energy using crystal structure from DFT object.
"""
function scf_energy(d::dftout; database = Dict(), smearing=0.01, grid = missing, conv_thr = 1e-5, iters = 75, mix = -1.0, mixing_mode=:pulay)

    return scf_energy(d.crys, smearing=smearing, grid = grid, conv_thr = conv_thr, iters = iters, mix = mix, mixing_mode=mixing_mode)

end

"""
    scf_energy(tbc::tbc_crys)

    SCF energy using crystal structure from TBC object.
"""
function scf_energy(tbc::tb_crys; smearing=0.01, grid = missing, e_den0 = missing, conv_thr = 1e-5, iters = 75, mix = -1.0, mixing_mode=:pulay)

    energy_tot, efermi, e_den, dq, VECTS, VALS, error_flag, tbc = SCF.scf_energy(tbc; smearing=smearing, grid = grid, e_den0 = e_den0, conv_thr = conv_thr, iters = iters, mix = mix, mixing_mode=mixing_mode)

    conv_flag = !error_flag
    if tbc.within_fit == false
        conv_flag = false
    end
    if conv_flag == true
        println("success, done")
    else
        println("error detected!!!!!!!!!!, done")
    end
    println()

    energy_tot = convert_energy(energy_tot)

    return energy_tot, tbc, conv_flag

end

end #end module


text:
using JuMP
using SCS

w1 = [0.4, 0.2, 0.2, 0.1, 0.1]
w2 = [0.6, 0.2, 0.2]

C = [ 3 5 2;
      0 2 5;
      1 1 3;
      8 4 3;
      7 6 5]

# The EMD is the sum(T .* C) where T[i, j] from each word/elem i in w1 to
# each word/elem j in w2 where T is the matrix that minimizes sum(T .* C)
# and sum(T, 1) = w2 and sum(T, 2) = w1.

# T = flow = emd(w1, w2, C)

m = Model()
N = length(w1)
M = length(w2)
@variable(m, 0 <= T[1:N,1:M] <= 1)
for j in 1:M
    @constraint(m, sum(T[i,j] for i=1:N) == w2[j])
end
for i in 1:N
    @constraint(m, sum(T[i,j] for j=1:M) == w1[i])
end
@objective(m, Min, sum(T .* C))
print(m)
@time status = solve(m)
println("Objective value: ", getobjectivevalue(m))
getvalue(T)


text:
@testset "test_dot.jl" begin
    a = Dualization.CanonicalVector{Int}(1, 3)
    @test collect(a) == [1, 0, 0]
    b = Dualization.CanonicalVector{Int}(2, 3)
    @test collect(b) == [0, 1, 0]
    set = MOI.PositiveSemidefiniteConeTriangle(2)
    @test MOI.Utilities.set_dot(a, b, set) == 0
    @test MOI.Utilities.set_dot(a, a, set) == 1
    @test MOI.Utilities.set_dot(b, b, set) == 2
    a = Dualization.CanonicalVector{Int}(1, 4)
    b = Dualization.CanonicalVector{Int}(2, 4)
    c = Dualization.CanonicalVector{Int}(3, 4)
    set = MOI.RootDetConeTriangle(2)
    @test MOI.Utilities.set_dot(a, b, set) == 0
    @test MOI.Utilities.set_dot(b, c, set) == 0
    @test MOI.Utilities.set_dot(c, a, set) == 0
    @test MOI.Utilities.set_dot(a, a, set) == 1
    @test MOI.Utilities.set_dot(b, b, set) == 1
    @test MOI.Utilities.set_dot(c, c, set) == 2
    set = MOI.PositiveSemidefiniteConeSquare(2)
    @test MOI.Utilities.set_dot(a, b, set) == 0
    @test MOI.Utilities.set_dot(b, c, set) == 0
    @test MOI.Utilities.set_dot(c, a, set) == 0
    @test MOI.Utilities.set_dot(a, a, set) == 1
    @test MOI.Utilities.set_dot(b, b, set) == 1
    @test MOI.Utilities.set_dot(c, c, set) == 1
end


text:
"""
	(type) Spherical
Representation of a complex value on the Riemann sphere.
"""
struct Spherical{T<:AbstractFloat} <: Number
	lat::T
	lon::T
end

# constructors
Spherical{T}(z::Spherical{T}) where {T<:AbstractFloat} = z

latitude(z::Number) = π/2 - 2*acot(abs(z))
function Spherical{T}(z::Number) where T<:AbstractFloat
	θ,ϕ = latitude(z),angle(z)
	Spherical{T}(convert(T,θ),convert(T,ϕ))
end

# Constructors without subtype
"""
	Spherical(latitude,azimuth)
Construct a spherical representation with given `latitude` in [-π/2,π/2] and `azimuth`.
"""
function Spherical(θ::Real,ϕ::Real)
	θ,ϕ = promote(float(θ),float(ϕ))
	Spherical{typeof(θ)}(θ,ϕ)
end
"""
Spherical(z)
Construct a spherical representation of the value `z`.
"""
Spherical(z::Number) = Spherical(latitude(z),angle(z))

# one and zero
one(::Type{Spherical{T}}) where T<:Real = Spherical{T}(zero(T),zero(T))
one(::Type{Spherical}) = one(Spherical{Float64})
zero(::Type{Spherical{T}}) where T<:Real = Spherical{T}(T(-π/2),zero(T))
zero(::Type{Spherical}) = zero(Spherical{Float64})

# conversion into standard complex
function Complex(z::Spherical{S}) where S<:Real
	if iszero(z)
		zero(Complex{S})
	else
		cot(π/4-z.lat/2) * exp(complex(zero(z.lon),z.lon))
	end
end

"""
	S2coord(u::Spherical)
Convert the spherical value to a 3-vector of coordinates on the unit sphere.
"""
S2coord(u::Spherical) = [cos(u.lat)*[cos(u.lon),sin(u.lon)];sin(u.lat)]

# basic arithmetic
function +(u::Spherical,v::Spherical)
	if isinf(u)
		isinf(v) ? NaN : u
	elseif isinf(v)
		v
	else
		Spherical(Complex(u)+Complex(v))  # faster way?
	end
end
-(u::Spherical) = Spherical(u.lat,cleanangle(u.lon+π))
-(u::Spherical,v::Spherical) = u + (-v)
*(u::Spherical,v::Spherical) = Spherical(Polar(u)*Polar(v))   # faster way?
inv(u::Spherical) = Spherical(-u.lat,cleanangle(-u.lon))
/(u::Spherical,v::Spherical) = u*inv(v)

# common complex overloads
angle(u::Spherical) = cleanangle(u.lon)
function abs(z::Spherical{T}) where T
	if iszero(z)
		zero(T)
	elseif isinf(z)
		T(Inf)
	else
		cot(π/4-z.lat/2)
	end
end
abs2(u::Spherical) = abs(u)^2
real(u::Spherical) = abs(u)*cos(u.lon)
imag(u::Spherical) = abs(u)*sin(u.lon)
conj(u::Spherical) = Spherical(u.lat,-u.lon)
sign(u::Spherical) = Spherical(zero(u.lat),u.lon)

# numerical comparisons
iszero(u::Spherical) = u.lat == convert(typeof(u.lat),-π/2)
isinf(u::Spherical) = u.lat == convert(typeof(u.lat),π/2)
isfinite(u::Spherical) = ~isinf(u)
isapprox(u::Spherical,v::Spherical;args...) = S2coord(u) ≈ S2coord(v)

# pretty output
show(io::IO,z::Spherical) = print(io,"(latitude = $(z.lat/pi)⋅π, angle = $(z.lon/pi)⋅π)")
show(io::IO,::MIME"text/plain",z::Spherical) = print(io,"Complex Spherical: ",z)


text:
# function warp(img::AstroImageMat, args...; kwargs...)
#     out = warp(parentt(img), args...; kwargs...)
#     return copyheaders(img, out)
# end

# Instead of using a datatype like N0f32 to interpret integers as fixed point values in [0,1],
# we use a mappedarray to map the native data range (regardless of type) to [0,1]
ImageBase.normedview(img::AstroImageMat{<:FixedPoint}) = img
function ImageBase.normedview(img::AstroImageMat{T}) where T
    imgmin, imgmax = extrema(skipmissingnan(img))
    Δ = abs(imgmax - imgmin)
    # Do not introduce NaNs if limits are identical
    if Δ == 0
        Δ = oneunit(imgmin)
    end
    normeddata = mappedarray(
        pix -> (pix - imgmin)/Δ,
        pix_norm -> convert(T, pix_norm*Δ + imgmin),
        img
    )
    return shareheader(img, normeddata)
end

"""
    clampednormedview(arr, (min, max))

Given an AbstractArray and limits `min,max` return a view of the array
where data between [min, max] are scaled to [0, 1] and datat outside that
range are clamped to [0, 1].

See also: normedview
"""
function clampednormedview(img::AbstractArray{T}, lims) where T
    imgmin, imgmax = lims
    Δ = imgmax - imgmin
    # Do not introduce NaNs if colorlimits are identical
    if Δ == 0
        Δ = oneunit(imgmin)
    end
    normeddata = mappedarray(
        pix -> clamp((pix - imgmin)/Δ, zero(pix), oneunit(pix)),
        img
    )
    return maybe_shareheader(img, normeddata)
end

# Restrict downsizes images by roughly a factor of two.
# We want to keep the wrapper but downsize the underlying array
# TODO: correct dimensions after restrict.
ImageBase.restrict(img::AstroImage, ::Tuple{}) = img
function ImageBase.restrict(img::AstroImage, region::Dims)
    restricted = restrict(parent(img), region)
    steps = cld.(size(img), size(restricted))
    newdims = Tuple(d[begin:s:end] for (d,s) in zip(dims(img),steps))
    return rebuild(img, restricted, newdims)
end

ImageBase.pixelspacing(img::AstroImage) = step.(dims(img))

# ImageContrastAdjustment
# function ImageContrastAdjustment.adjust_histogram(::Type{T},
#     img::AstroImageMat,
#     f::Images.ImageContrastAdjustment.AbstractHistogramAdjustmentAlgorithm,
#     args...; kwargs...) where T
#     out = similar(img, axes(img))
#     adjust_histogram!(out, img, f, args...; kwargs...)
#     return out
# end

# """
#     ImageMetadata.parent(img::AstroImage)

# Returns the underlying wrapped array of `img`.
# """
# ImageMetadata.parent(img::AstroImage) = getfield(img, :data)


text:
using Test
using PropCheck
using PropCheck: getSubtypes, numTests

"""
Tests the fallback shrinking for numbers to shrink towards zero.
"""
function numsToZero(T)
    # are there other properties of the shrinkers themselves we can test?
    check(igen(T)) do x
        shrunks = shrink(x)
        if x == zero(T)
            isempty(shrunks)
        else
            !isempty(shrunks) && all( y -> zero(T) <= abs(y) < abs(x), shrunks)
        end
    end
end

"""
Tests that when a given predicate holds for the parent, it also holds for its subtrees (or at least the first 100).
"""
function predicateHoldsForSubtrees(p, T)
    g = filter(p, igen(T))
    toCheck = filter!(!isnothing, [iterate(g) for _ in 1:numTests[]])
    iszero(length(toCheck)) && @warn "no values to check for $T"
    all(toCheck) do x
        all(p ∘ root, Iterators.take(subtrees(first(x)), 100))
    end
end

function throwingProperty(x)
    if x < 5
        return true
    else
        throw(ArgumentError("x not smaller than 5"))
    end
end

const numTypes = union(getSubtypes(Integer), getSubtypes(AbstractFloat), (Float64, Float32, Float16))

@testset "All Tests" begin
    @testset "numsToZero" begin
        @testset "$T" for T in numTypes
            @test numsToZero(T) broken=(T == BigInt || T <: AbstractFloat)
        end
    end
    @testset "filter predicates hold for shrunk values" begin
        for p in (iseven, isodd)
            @testset "$p($T)" for T in (UInt8, UInt16, UInt32, UInt64)
                @test predicateHoldsForSubtrees(p, T) broken=(T == BigInt || T <: AbstractFloat)
            end
        end
    end
    @testset "random vectors are sorted" begin
        @test check(issorted, PropCheck.vector(igen(20), igen(UInt8))) == [0x1, 0x0]
    end
    @testset "all even numbers are less than 5" begin
        @test check(<(5), filter(iseven, igen(UInt8))) == 0x6
    end
    @testset "there are only even numbers" begin
        @test check(iseven, igen(UInt8)) == 0x1
    end
    @testset "throwing properties still shrink" begin
        @test check(throwingProperty, igen(UInt8)) == (0x05, ArgumentError("x not smaller than 5"))
    end
end


text:
mutable struct PlaneWaveMW{T,P}
  direction::P
  polarisation::P
  wavenumber::T
  amplitude::T
end

function PlaneWaveMW(d,p,k,a = 1)
    T = promote_type(eltype(d), eltype(p), typeof(k), typeof(a))
    P = similar_type(typeof(d), T)
    PlaneWaveMW{T,P}(d,p,k,a)
end

"""
    planewavemw3d(;direction, polarization, wavenumber[, amplitude=1])

Create a plane wave solution to Maxwell's equations.
"""
planewavemw3d(;
    direction    = error("missing arguement `direction`"),
    polarization = error("missing arguement `polarization`"),
    wavenumber   = error("missing arguement `wavenumber`"),
    amplitude    = 1,
    ) = PlaneWaveMW(direction, polarization, wavenumber, amplitude)

function (e::PlaneWaveMW)(x)
  k = e.wavenumber
  d = e.direction
  u = e.polarisation
  a = e.amplitude
  a * exp(-im * k * dot(d, x)) * u
end

function curl(field::PlaneWaveMW)
  k = field.wavenumber
  d = field.direction
  u = field.polarisation
  a = field.amplitude
  v = d × u
  b = -a * im * k
  PlaneWaveMW(d, v, k, b)
end

*(a::Number, e::PlaneWaveMW) = PlaneWaveMW(e.direction, e.polarisation, e.wavenumber, a*e.amplitude)

mutable struct CrossTraceMW{F} <: Functional
  field::F
end

mutable struct TangTraceMW{F} <: Functional
  field::F
end

cross(::NormalVector, p::Function) = CrossTraceMW(p)
cross(::NormalVector, p::PlaneWaveMW) = CrossTraceMW(p)
cross(t::CrossTraceMW, ::NormalVector) = TangTraceMW(t.field)

function (ϕ::CrossTraceMW)(p)
  F = ϕ.field
  x = cartesian(p)
  n = normal(p)
  return n × F(x)
end

function (ϕ::TangTraceMW)(p)
  F = ϕ.field
  x = cartesian(p)
  n = normal(p)
  return (n × F(x)) × n
end

integrand(::TangTraceMW, gx, ϕx) = gx[1] ⋅ ϕx
integrand(::CrossTraceMW, test_vals, field_val) = test_vals[1] ⋅ field_val

struct NDotTrace{F} <: Functional
  field::F
end

(ϕ::NDotTrace)(p) = dot(normal(p), ϕ.field(cartesian(p)))
integrand(::NDotTrace, g, ϕ) = dot(g.value, ϕ)
LinearAlgebra.dot(::NormalVector, f) = NDotTrace(f)


text:
#======================================#
#####  ┏━╸┏━┓┏━┓╺┳┓   ┏━╸┏━┓┏┳┓┏━╸ #####
#####  ┃  ┣━┫┣┳┛ ┃┃   ┃╺┓┣━┫┃┃┃┣╸  #####
#####  ┗━╸╹ ╹╹┗╸╺┻┛   ┗━┛╹ ╹╹ ╹┗━╸ #####
#======================================#

include("./card_gen.jl")
include("./utils.jl")
include("./computer.jl")
include("./moves.jl")
include("./user.jl")

#======================================#
#        Other File's Functions        #
#======================================#

# Cards tup - stores all card stats
# INFO ABOUT Cards

# FUNCTION IN card_gen
# GenCard(amntCards ::Int64) - Generate card stats of a given number of cards

# FUNCTIONS IN user
# print_starting_cards() - print user cards
# pick_card() - allows user to pick their card (returns userCardChoice (int showing the number of card user chose))

# FUNCTIONS IN computer
# comp_choosing_card() - comp_thinking() + computer picks the last stats of the Cards tup (returns compCardStats (a tup containing info about the card choosen by the computer))

# FUNCTIONS IN utils
# thinking(message ::String, end_message ::String) - show the little animation of thinking
# determine_first_player() - picks and shows the first player, also returns the first player as an int

# FUNCTIONS IN moves
# dodge(currentDodgePercentage<>) - no idea tbh
# defend()
# attack(minDmg, maxDmg, critChance, critDmgIncrease)

#======================================#
#              Functions               #
#======================================#

# play_again() - play again
# play_moves() - play moves

function playAgain()
  # play again
  print("Do you want to play again (y/N): ")
  playAgain = readline()

  if playAgain in ["Y", "y"]
    # Clear all the previous card stats
    pop_many!(x, n) = [pop!(x) for _ in 1:n]
    pop_many!(Cards, length(Cards))
    return true
  else
    println("Goodbye :)")
    return false
  end
end

#======================================#
#              Game Loop               #
#======================================#

userChances = Dict("userChanceToDodge"=>nothing, "userDefendAmount"=>nothing)
compChances = Dict("compChanceToDodge"=>nothing, "compDefendAmount"=>nothing)

try
  playing = true
  Cards = []
  while playing
    @label start_of_game

    # generate and display cards to pick from
    @async Cards = GenCard(4, Cards)
    print("testing:::: $Cards")
    @sync @async print_starting_cards(Cards)

    # get card stats for comp and user
    userCardStats = Cards[pick_card()]
    compCardStats = comp_choosing_card()

    # get the first player
    first_player = determine_first_player()

    display_info(userCardStats, compCardStats, userChances, compChances)

    # Run at end of the game
    if playAgain()
        @goto start_of_game
    else
      return playing = false
    end
  end
catch e
  println("Error: $e")
end

max_stars_repo_name: NarutoXY/austroporean-cards

max_issues_repo_name: NarutoXY/austroporean-cards

max_forks_repo_name: NarutoXY/austroporean-cards


text:
mutable struct Edges
    em::Vector{Vector{Float64}}
    Edges() = new(Vector{Vector{Float64}}())
end

size(e::Edges) = size(e.em, 1)
getindex(e::Edges, i::Union{Int,UnitRange}) = e.em[i]

function addpoint!(this::Edges, data::Vector{Float64})
    while size(data, 1) < 4
        push!(data, 1.0)
    end
    push!(this.em, data[1:4])
end

function transform!(this::Edges, transmatrix::Matrix{Float64})
    size(this) < 1 && return
    tmpmat = transmatrix * hcat(this.em...)
    this.em = [tmpmat[:, i] for i in 1:size(tmpmat, 2)]
    this
end

*(x::Edges, y::Matrix{Float64}) = transform!(x, y)

function addedge!(this::Edges, p1::Vector{Float64}, p2::Vector{Float64})
    addpoint!(this, p1)
    addpoint!(this, p2)
end

function addpolygon!(this::Edges, p1::Vector{Float64}, p2::Vector{Float64}, p3::Vector{Float64})
    addpoint!(this, p1)
    addpoint!(this, p2)
    addpoint!(this, p3)
end

function addbox!(this::Edges, tl::Vector{Float64}, width::Float64, height::Float64,
                 depth::Float64)
    br = tl + [width, -height, -depth] # topleft, bottomright
    # front, back right, left, top, bottom
    addpolygon!(this, [tl[1], tl[2], tl[3]], [br[1], br[2], tl[3]], [br[1], tl[2], tl[3]])
    addpolygon!(this, [tl[1], tl[2], tl[3]], [tl[1], br[2], tl[3]], [br[1], br[2], tl[3]])

    addpolygon!(this, [br[1], tl[2], br[3]], [tl[1], br[2], br[3]], [tl[1], tl[2], br[3]])
    addpolygon!(this, [br[1], tl[2], br[3]], [br[1], br[2], br[3]], [tl[1], br[2], br[3]])

    addpolygon!(this, [br[1], tl[2], tl[3]], [br[1], br[2], br[3]], [br[1], tl[2], br[3]])
    addpolygon!(this, [br[1], tl[2], tl[3]], [br[1], br[2], tl[3]], [br[1], br[2], br[3]])

    addpolygon!(this, [tl[1], tl[2], br[3]], [tl[1], br[2], tl[3]], [tl[1], tl[2], tl[3]])
    addpolygon!(this, [tl[1], tl[2], br[3]], [tl[1], br[2], br[3]], [tl[1], br[2], tl[3]])

    addpolygon!(this, [tl[1], tl[2], br[3]], [br[1], tl[2], tl[3]], [br[1], tl[2], br[3]])
    addpolygon!(this, [tl[1], tl[2], br[3]], [tl[1], tl[2], tl[3]], [br[1], tl[2], tl[3]])

    addpolygon!(this, [tl[1], br[2], tl[3]], [br[1], br[2], br[3]], [br[1], br[2], tl[3]])
    addpolygon!(this, [tl[1], br[2], tl[3]], [tl[1], br[2], br[3]], [br[1], br[2], br[3]])
end

function addcircle!(this::Edges, center::Vector{Float64}, radius::Float64, steps::Int)
    segstart = center + [radius, 0, 0]
    for t in [i/steps for i in 1:steps]
        segend = center + radius * [cos(2pi * t), sin(2pi * t), 0]
        addedge!(this, segstart, segend)
        segstart = segend
    end
end

function addsphere!(this::Edges, center::Vector{Float64}, radius::Float64, steps::Int)
    points = Vector{Vector{Float64}}()

    for rott in [i/steps for i in 0:steps], cirt in [i/steps for i in 0:steps]
        push!(points, center + radius * [cos(pi * cirt),
                                         sin(pi * cirt) * cos(2pi * rott),
                                         sin(pi * cirt) * sin(2pi * rott)])
    end

    ns = steps + 1
    for latt in 0:steps - 1, long in 0:steps - 1
        i = latt * ns + long
        p = [i, i + 1, (i + 1 + ns) % length(points), (i + ns) % length(points)] + 1
        addpolygon!(this, points[p[1]], points[p[2]], points[p[3]])
        addpolygon!(this, points[p[1]], points[p[3]], points[p[4]])
    end
end

function addtorus!(this::Edges, center::Vector{Float64}, irad::Float64, orad::Float64,
                   steps::Int)
    # irad: radius of a "slice"         orad: radius of the entire torus
    points = Vector{Vector{Float64}}()

    for rott in [i/steps for i in 0:steps], cirt in [i/steps for i in 0:steps]
        push!(points, center + [cos(2pi * rott) * (irad * cos(2pi * cirt) + orad),
                                irad * sin(2pi * cirt),
                                -sin(2pi * rott) * (irad * cos(2pi * cirt) + orad)])
    end

    ns = steps + 1
    for latt in 0:steps - 1, long in 0:steps - 1
        i = latt * ns + long
        p = [i, i + 1, (i + 1 + ns) % ns^2, (i + ns) % ns^2] + 1

        addpolygon!(this, points[p[1]], points[p[4]], points[p[3]])
        addpolygon!(this, points[p[1]], points[p[3]], points[p[2]])
    end
end

function addcurve!(this::Edges, p1::Vector{Float64}, p2::Vector{Float64},
                   p3::Vector{Float64}, p4::Vector{Float64}, steps::Int, ctype::Symbol)
    xcoefs = CURVES[ctype] * [p1[1]; p2[1]; p3[1]; p4[1]]
    ycoefs = CURVES[ctype] * [p1[2]; p2[2]; p3[2]; p4[2]]

    for t in [i/steps for i in 1:steps]
        # at^3 + bt^2 + ct + d
        p2 = [xcoefs[1] * t^3 + xcoefs[2] * t^2  + xcoefs[3] * t + xcoefs[4],
              ycoefs[1] * t^3 + ycoefs[2] * t^2  + ycoefs[3] * t + ycoefs[4],
              0]
        addedge!(this, p1, p2)
        p1 = p2
    end
end

function drawem!(this::Edges, display::IBuffer, color::Vector{Int})
    for i in 1:2:(size(this) - 1)
        drawline!(display, this[i], this[i + 1], color)
    end
end

function drawpm!(this::Edges, display::IBuffer, view::Vector{Float64},
                 cambient::Vector{Int}, lights::Vector{Dict{Symbol, Any}},
                 reflect::Dict{Symbol, Vector{Float64}})
    # preallocate arrays
    v     = normalize(view)
    n     = Vector{Float64}(3)
    l     = Vector{Float64}(3)
    pb    = Vector{Float64}(4)
    pm    = Vector{Float64}(4)
    pt    = Vector{Float64}(4)
    color = Vector{Float64}(3)
    dleft = Vector{Float64}(4)
    drig1 = Vector{Float64}(4)
    drig2 = Vector{Float64}(4)
    tmpa  = Vector{Float64}(4)

    # normalize all lights
    for lig in lights
        normalize!(lig[:location])
    end

    for i in 1:3:(size(this) - 2)
        # calculate the normal to see if we need to draw
        n[1:3] = cross(this[i + 1] - this[i], this[i + 2] - this[i])
        vecdot(n, view) > 0 || continue

        # sort by and round y values for cleaner spheres
        pb[1:4], pm[1:4], pt[1:4] = sort(this[i:i + 2], by = x -> x[2])
        pb[2] = round(Int, pb[2])
        pm[2] = round(Int, pm[2])
        pt[2] = round(Int, pt[2])

        # determine shading of a polygon
        normalize!(n)
        am, di, sp = cambient .* reflect[:ambient], 0, 0

        for lig in lights
            l[1:3] = lig[:location]
            vnl = vecdot(n, l)
            di += lig[:color] .* reflect[:diffuse] * vnl
            sp += lig[:color] .* reflect[:specular] * vecdot(v, 2n * vnl - l)^SPEC_EXP
        end
        @. color[1:3] = min(255, max(0, am) + max(0, di) + max(0 + sp))

        # fill in polygons with scanline conversion
        f(x, y) = y < 0.1 ? 0 : x / y   # avoid division by 0
        dleft[1:4] = [f(pt[1] - pb[1], pt[2] - pb[2]), 1, f(pt[3] - pb[3], pt[2] - pb[2]), 0]
        drig1[1:4] = [f(pm[1] - pb[1], pm[2] - pb[2]), 1, f(pm[3] - pb[3], pm[2] - pb[2]), 0]
        drig2[1:4] = [f(pt[1] - pm[1], pt[2] - pm[2]), 1, f(pm[3] - pb[3], pm[2] - pb[2]), 0]

        tmp = pm[2] - pb[2]
        for y in 1:tmp
            drawline!(display, pb + y * dleft, pb + y * drig1, color)
        end

        tmpa[1:4] = pb + tmp * dleft
        for y in 0:pt[2] - pm[2]
            drawline!(display, tmpa + y * dleft, pm + y * drig2, color)
        end
    end
end

max_stars_repo_name: AnAverageHuman/jlgraphics

max_issues_repo_name: AnAverageHuman/jlgraphics

max_forks_repo_name: AnAverageHuman/jlgraphics


text:
@defcomp grosseconomy begin
    regions = Index()                           #Note that a regional index is defined here

    YGROSS  = Variable(index=[time, regions])   #Gross output
    K       = Variable(index=[time, regions])   #Capital
    l       = Parameter(index=[time, regions])  #Labor
    tfp     = Parameter(index=[time, regions])  #Total factor productivity
    s       = Parameter(index=[time, regions])  #Savings rate
    depk    = Parameter(index=[regions])        #Depreciation rate on capital - Note that it only has a region index
    k0      = Parameter(index=[regions])        #Initial level of capital
    share   = Parameter()                       #Capital share

    tester  = Parameter(index=[time, regions])  # test for reshaped Dirichlet distribution

    function run_timestep(p, v, d, t)
    # Note that the regional dimension is defined in d and parameters and variables are indexed by 'r'

        # Define an equation for K
        for r in d.regions
            if is_first(t)
                v.K[t,r] = p.k0[r]
            else
                v.K[t,r] = (1 - p.depk[r])^5 * v.K[t-1,r] + v.YGROSS[t-1,r] * p.s[t-1,r] * 5
            end
        end

        # Define an equation for YGROSS
        for r in d.regions
            v.YGROSS[t,r] = p.tfp[t,r] * v.K[t,r]^p.share * p.l[t,r]^(1-p.share)
        end
    end
end


text:
"""
Common interface for all modes and key lengths
"""
function decrypt end

decrypt(ciphertext::CT{IV,CBC}, cipher::Cipher{CBC,ct,kt}; remove_pad=true) where {IV,ct,kt} = AESCBC_D(ciphertext.data, ciphertext.iv, cipher.key, cipher.cache; remove_pad=remove_pad)
decrypt(ciphertext::CT{IV,CTR}, cipher::Cipher{CTR,ct,kt}) where {IV,ct,kt} = AESCTR_D(ciphertext.data, ciphertext.iv, cipher.key, cipher.cache)
decrypt(ciphertext::CT{IV,ECB}, cipher::Cipher{ECB,ct,kt}) where {IV,ct,kt} = AESECB_D(ciphertext.data, cipher.key, cipher.cache)

decrypt!(plaintext, ciphertext::CT{IV,CBC}, cipher::Cipher{CBC,ct,kt}; remove_pad=true) where {IV,ct,kt} = AESCBC_D!(plaintext, ciphertext.data, ciphertext.iv, cipher.key, cipher.cache; remove_pad=remove_pad)
decrypt!(plaintext, ciphertext::CT{IV,CTR}, cipher::Cipher{CTR,ct,kt}) where {IV,ct,kt} = AESCTR_D!(plaintext, ciphertext.data, ciphertext.iv, cipher.key, cipher.cache)
decrypt!(plaintext, ciphertext::CT{IV,ECB}, cipher::Cipher{ECB,ct,kt}) where {IV,ct,kt} = AESECB_D!(plaintext, ciphertext.data, cipher.key, cipher.cache)


text:
struct ConstMean{T}
    c::T
end

(m::ConstMean)(x) = m.c.*ones(size(x, 2))

struct LinearMean{T}
    a::T
    b::T
end

function (m::LinearMean)(x)
    F, N = size(x, 1), size(x, 2)
    y = Array{eltype(m.a)}(undef, N)
    for i in 1:N
        tmp = zero(eltype(m.a))
        for j in 1:F
            tmp += m.a[j]*x[j,i]
        end
        y[i] = tmp
    end
    return y
end

max_stars_repo_name: HamletWantToCode/DiffGPFlux

max_issues_repo_name: HamletWantToCode/DiffGPFlux


text:
# example 2: string hashing with a seed

include("../src/Murmur3.jl")
using Murmur3

value = "Hello world"
seed = 42

hash = Murmur3.x86.hash32(value, seed)
println("x86  32-bit (seed: $seed): \"$value\" => $(hex(hash))")

hash = Murmur3.x86.hash128(value, seed)
println("x86 128-bit (seed: $seed): \"$value\" => $(hex(hash))")

hash = Murmur3.x64.hash128(value, seed)
println("x64 128-bit (seed: $seed): \"$value\" => $(hex(hash))")


text:
# this doesn't work yet, but we can load images using other tools anyway
## "Read image from file"
## function read_image(path::AbstractString, typ::Int)
##     # @assert rows*cols == length(data) "rows * cols != length of data!"
##     T = Float64
##     # typ = CCV_IO_RGB_COLOR | CCV_IO_ANY_FILE
##     typ = CCV_IO_GRAY | CCV_IO_ANY_FILE
##     matout = Ref{Ptr{CCVDenseMatrixType{T}}}()
##     path = "minions.jpg"
##     # scanline = sizeof(T)
##     # typ = ccv_type(T) | CCV_IO_ANY_FILE
##     code = ccall((:ccv_read_impl, libccv), Cint,
##                  (Ptr{Cchar}, Ptr{Ptr{CCVDenseMatrixType{T}}},
##                   Cint, Cint, Cint, Cint),
##                  path, matout, typ, 0, 0, 0)
##     if code != CCV_IO_FINAL
##         throw(ErrorException("ccv_read_impl() exited with error: $code"))
##     else
##         ptr = matout[]
##         return CCVDenseMatrix{T}(ptr)
##     end
## end

if Pkg.installed("Images") != nothing

    using Images, Colors

    Base.convert{T}(::Type{Matrix{Float64}}, img::Image{Gray{T}}) =
        convert(Array{Float64, 2}, convert(Array, img))

    Base.convert(::Type{Matrix{Float64}}, img::Image) =
        convert(Array{Float64, 2}, convert(Image{Gray}, img))

    Base.convert(::Type{Array{Float64, 3}}, img::Image) =
        convert(Array{Float64}, separate(img).data)

    function CCVDenseMatrix(img::Image)
        arr = convert(Matrix{Float64}, img)
        return CCVDenseMatrix(arr)
    end

end


text:
##############
# Issue type #
##############

type Issue <: GitLabType
    id::Nullable{Int}
    iid::Nullable{Int}
    project_id::Nullable{Int}
    title::Nullable{GitLabString}
    description::Nullable{GitLabString}
    state::Nullable{GitLabString}
    created_at::Nullable{Dates.DateTime}
    updated_at::Nullable{Dates.DateTime}
    ## labels::Nullable{Vector{Dict}}
    labels::Nullable{Vector{GitLabString}}
    milestone::Nullable{GitLabString}
    assignee::Nullable{Owner}
    author::Nullable{Owner}
    subscribed::Nullable{Bool}
    user_notes_count::Nullable{Int}

#=
    closed_by::Nullable{Owner}
    closed_at::Nullable{Dates.DateTime}
    pull_request::Nullable{PullRequest}
    url::Nullable{HttpCommon.URI}
    html_url::Nullable{HttpCommon.URI}
    labels_url::Nullable{HttpCommon.URI}
    comments_url::Nullable{HttpCommon.URI}
    events_url::Nullable{HttpCommon.URI}
    locked::Nullable{Bool}
=#
end

Issue(data::Dict) = json2gitlab(Issue, data)
Issue(id::Int) = Issue(Dict("id" => id))

namefield(issue::Issue) = issue.id

###############
# API Methods #
###############

function issue(repo::Repo, issue::Int; options...)
    result = gh_get_json("/api/v3/projects/$(get(repo.id))/issues/$(issue)"; options...)
    return Issue(result)
end

function issues(repo::Repo; options...)
    results, page_data = gh_get_paged_json("/api/v3/projects/$(get(repo.id))/issues"; options...)
    return map(Issue, results), page_data
end

function create_issue(repo::Repo; options...)
    result = gh_post_json("/api/v3/projects/$(get(repo.id))/issues"; options...)
    return Issue(result)
end

function edit_issue(repo::Repo, issue::Int; options...)
    result = gh_put_json("/api/v3/projects/$(get(repo.id))/issues/$(issue)"; options...)
    return Issue(result)
end

function delete_issue(repo::Repo, issue::Int; options...)
    result = gh_delete_json("/api/v3/projects/$(get(repo.id))/issues/$(issue)"; options...)
    return Issue(result)
end


text:
using Statistics
using Serialization
using Printf
using ArgParse

const JULIAVER = Base.julia_cmd()[1]

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--runs", "-n"
            help = "number of threads"
            arg_type = Int
            default = 10
        "--threads", "-t"
            help = "another option with an argument"
            arg_type = Int
            default = 1
        "--bench"
            help = "if specified, runs a single benchmark"
            arg_type = String
            default = "all"
    end

    return parse_args(s)
end

gctime(stat) = stat.total_time

function run_one_bench(runs, threads, file)
    value = []
    times = []
    stats = []
    for _ in 1:runs
        r = open(deserialize, `$JULIAVER --project=. --threads=$threads $file SERIALIZE`)
        push!(value, r.value)
        push!(times, r.times)
        push!(stats, r.stats)
    end
    @printf("run time: %0.0fms min, %0.0fms max %0.0fms median\n",
       minimum(times) / 1_000_000,
       maximum(times) / 1_000_000,
       median(times) / 1_000_000)
    time = map(gctime, stats)
    @printf("gc time: %0.0fms min, %0.0fms max, %0.0fms median\n",
       minimum(time) / 1_000_000,
       maximum(time) / 1_000_000,
       median(time) / 1_000_000)
end

function run_all_benches(runs, threads)
    for category in readdir()
        @show category
        cd(category)
        for file in readdir()
            endswith(file, ".jl") || continue
            @show file
            run_one_bench(runs, threads, file)
        end
    end
    cd("..")
end

function main()
    args = parse_commandline()
    JULIAVER = Base.julia_cmd()[1]
    bench, runs, threads = args["bench"], args["runs"], args["threads"]

    dir = joinpath(@__DIR__, "benches")
    cd(dir)
    if bench == "all"
        run_all_benches(runs, threads)
    else
        path, file = splitpath(bench)
        cd(path)
        run_one_bench(runs, threads, file)
    end
end

main()

max_stars_repo_name: chflood/GCBenchmarks

max_issues_repo_name: chflood/GCBenchmarks

max_forks_repo_name: chflood/GCBenchmarks


text:
module HarmonicOscillatorProblem

    using GeometricIntegrators.Equations

    export harmonic_oscillator_ode, harmonic_oscillator_iode, harmonic_oscillator_pode, harmonic_oscillator_sode,
           harmonic_oscillator_dae, harmonic_oscillator_idae, harmonic_oscillator_pdae

    const Δt = 0.1
    const nt = 10

    const k = 0.5
    const ω = √k

    t₀=0.0
    q₀=[0.5, 0.0]
    z₀=[0.5, 0.0, 0.5]

    function ϑ(q)
        p = zero(q)
        p[1] = q[2]
        p[2] = 0
        return p
    end

    function hamiltonian(t, q)
        q[2]^2 / 2 + k * q[1]^2 / 2
    end

    p₀=ϑ(q₀)

    A = sqrt(q₀[2]^2 / k + q₀[1]^2)
    ϕ = asin(q₀[1] / A)

    refq = A * sin(ω * Δt * nt + ϕ)
    refp = ω * Δt * nt * A * cos(ω * Δt * nt + ϕ)
    refx = [refq, refp]

    function oscillator_ode_v(t, x, v)
        v[1] = x[2]
        v[2] = -k*x[1]
        nothing
    end

    function harmonic_oscillator_ode(x₀=q₀)
        @assert size(x₀,1) == 2
        ODE(oscillator_ode_v, x₀; h=hamiltonian)
    end

    function oscillator_pode_v(t, q, p, v)
        v[1] = p[1]
        nothing
    end

    function oscillator_pode_f(t, q, p, f)
        f[1] = -k*q[1]
        nothing
    end

    function harmonic_oscillator_pode(q₀=[q₀[1]], p₀=[p₀[1]])
        @assert size(q₀,1) == size(p₀,1) == 1
        PODE(oscillator_pode_v, oscillator_pode_f, q₀, p₀)
    end

    function oscillator_sode_v_1(t, q, v, h)
        v[1] = q[1] + h * q[2]
        v[2] = q[2]
        nothing
    end

    function oscillator_sode_v_2(t, q, v, h)
        v[1] = q[1]
        v[2] = q[2] - h * k*q[1]
        nothing
    end

    function harmonic_oscillator_sode(q₀=q₀)
        SODE((oscillator_sode_v_1, oscillator_sode_v_2), q₀)
    end

    function oscillator_iode_ϑ(t, q, p)
        p[1] = q[2]
        p[2] = 0
        nothing
    end

    function oscillator_iode_ϑ(t, q, v, p)
        oscillator_iode_ϑ(t, q, p)
    end

    function oscillator_iode_f(t, q, v, f)
        f[1] = -k*q[1]
        f[2] = v[1] - q[2]
        nothing
    end

    function oscillator_iode_g(t, q, λ, g)
        g[1] = 0
        g[2] = λ[1]
        nothing
    end

    function oscillator_iode_v(t, q, v)
        v[1] = q[2]
        v[2] = -k*q[1]
        nothing
    end

    function oscillator_iode_v(t, q, p, v)
        oscillator_iode_v(t, q, v)
    end

    function harmonic_oscillator_iode(q₀=q₀, p₀=ϑ(q₀))
        @assert size(q₀,1) == size(p₀,1) == 2
        IODE(oscillator_iode_ϑ, oscillator_iode_f,
             oscillator_iode_g, q₀, p₀;
             h=hamiltonian, v=oscillator_iode_v)
    end

    function oscillator_dae_v(t, z, v)
        v[1] = z[2]
        v[2] = -k*z[1]
        v[3] = z[2] - k*z[1]
        nothing
    end

    function oscillator_dae_u(t, z, λ, u)
        u[1] = -λ[1]
        u[2] = -λ[1]
        u[3] = +λ[1]
    end

    function oscillator_dae_ϕ(t, z, λ, ϕ)
        ϕ[1] = z[3] - z[1] - z[2]
    end

    function harmonic_oscillator_dae(z₀=z₀, λ₀=[zero(eltype(z₀))])
        @assert size(z₀,1) == 3
        @assert size(λ₀,1) == 1
        DAE(oscillator_dae_v, oscillator_dae_u, oscillator_dae_ϕ, z₀, λ₀)
    end

    function oscillator_idae_u(t, q, p, λ, u)
        u[1] = λ[1]
        u[2] = λ[2]
        nothing
    end

    function oscillator_idae_g(t, q, p, λ, g)
        g[1] = 0
        g[2] = λ[1]
        nothing
    end

    function oscillator_idae_ϕ(t, q, p, ϕ)
        ϕ[1] = p[1] - q[2]
        ϕ[2] = p[2]
        nothing
    end

    function harmonic_oscillator_idae(q₀=q₀, p₀=ϑ(q₀), λ₀=zero(q₀))
        @assert size(q₀,1) == size(p₀,1) == size(λ₀,1) == 2
        IDAE(oscillator_iode_f, oscillator_iode_ϑ,
             oscillator_idae_u, oscillator_idae_g,
             oscillator_idae_ϕ, q₀, p₀, λ₀;
             v=oscillator_iode_v)
    end

    function oscillator_pdae_v(t, q, p, v)
        v[1] = q[2]
        v[2] = -k*q[1]
        nothing
    end

    function oscillator_pdae_f(t, q, p, f)
        f[1] = -k*q[1]
        f[2] = p[1] - q[2]
        nothing
    end

    function harmonic_oscillator_pdae(q₀=q₀, p₀=ϑ(q₀), λ₀=zero(q₀))
        @assert size(q₀,1) == size(p₀,1) == 2
        PDAE(oscillator_pdae_v, oscillator_pdae_f,
             oscillator_idae_u, oscillator_idae_g,
             oscillator_idae_ϕ, q₀, p₀, λ₀)
    end

end


text:
module ValidateEntropies
using Test
using Statistics

using Diversity
using Diversity.ShortNames
using StatsBase

@testset "renyientropy" begin
    @testset "Random renyientropy $i" for i in 1:20
        types = rand(1:(i*10))
        pop = rand(types)
        pop[pop .< median(pop)/2] .= 0.0
        pop ./= sum(pop)
        meta = Metacommunity(pop)
        g = Γ(meta)
        na = ᾱ(meta)
        ra = α(meta)
        qs = [rand(7)*10..., 0, 1, Inf]
        entropies = map(q -> renyientropy(pop, q), qs)
        @test log.(subdiv(na, qs)[!,:diversity]) ≈ entropies
        @test log.(subdiv(ra, qs)[!,:diversity]) ≈ entropies
        @test log.(subdiv(g, qs)[!,:diversity]) ≈ entropies
        @test log.(qD(pop, qs)) ≈ entropies
        @test log.(qDZ(pop, qs, UniqueTypes(types))) ≈ entropies
    end
end

end


text:
abstract Matrix

type LowRankMatrix{T<:Number} <: Matrix
    ## a rank-1 matrix type
    # A = u*v'
    # m is size of matrix
    U::Array{T, 2}
    V::Array{T, 2}
    n::Int64
    m::Int64
    k::Int64
end

function LowRankMatrix{T<:Number}(U::Array{T, 2}, V::Array{T, 2})
    @assert size(U)[2] == size(V)[1] # we're only dealing with square matrices
    return LowRankMatrix(U, V, size(U)[1], size(V)[2], size(V)[1])
end

function randSVD{T<:Number}(A::Array{T,2}, k::Int64, oversample = 5)
    # function to compute the randomized SVD of the matrix A
    # we hope to extract only the first k biggest eigenvalues
    (n,m) = size(A)
    #number of random vectors
    nRandVecs = min(k+oversample,min(n,m))
    randVecs = randn(nRandVecs, n);
    # we multiply A by the random vectors
    R = randVecs*A;
    # we compute the row space of A using a skinny svd
    Q = svd(R')[1]
    P = A*Q;
    (U, S, W) = svd(P)
    V = Q*W
    return (U[:, 1:k], S[1:k], V[:,1:k])
end

function compressLowRankMatrix{T<:Number}(A::Array{T, 2}, epsilon::Float64 )
    # function to compress into a low rank matrix
    (U, S,V) = svd(A)
    # computing the epsilon rank
    epsRank = length(find(S.>epsilon));
    return LowRankMatrix(U[:,1:epsRank], spdiagm(S[1:epsRank])*(V[:,1:epsRank])',
                         size(A)[1], size(A)[2], epsRank)
end

import Base.*
function *{T<:Number}(M::LowRankMatrix{T}, v::Array{T,1})
    # Matrix vector multiplication overload
    @assert M.m == length(v) # make sure vector is right size for matrix
    c = M.V*v
    return M.U*c
end
function *{T<:Number}(M::LowRankMatrix{T}, N::Array{T,2})
    # Matrix-Matrix product overload
    @assert M.m == size(N)[1] # make sure vector is right size for matrix
    C = M.V*N
    return M.U*C
end

;


text:
_diffusion_jacobian(model::AbstractSDE{D,M}, t0, x0) where {D,M} = ForwardDiff.jacobian(x->diffusion(model, t0, x), x0)
_diffusion_jacobian(model::AbstractSDE{1,M}, t0, x0) where {M} = ForwardDiff.derivative(x->diffusion(model, t0, x), x0)

function step(model::AbstractSDE{D,1}, scheme::Milstein, t0, x0, Δw) where D
  μ = drift(model, t0, x0)
  σ = diffusion(model, t0, x0)
  ∂σ = _diffusion_jacobian(model, t0, x0)
  x = x0 + μ * scheme.Δt + σ * Δw + ∂σ * σ * (Δw.^2 - scheme.Δt) / 2
  x
end


text:
using MClassification, Test, DataFrames

@testset "Predict Array{Number, 2}" begin
    initial_samples = 4
    X = [5.2603 0.40807;
    0.51512 -0.7065;
    4.721 4.5322;
    0.36285 1.0328;
    0.89978 1.6674;
    2.7659 0.74216;
    1.4156 2.2328;
    1.0293 0.12319;
    3.3528 3.261]
    y = MClassification.categorical([2, 1, 2, 1, 1, 1, 1, 1, 2])
    model = MClassification.MClassifier(r_limit = 0.1)

    fitresult, cache, report = MClassification.fit(model, 0, X[1:initial_samples, :], y[1:initial_samples])

    y_predicted = MClassification.update_predict(model, fitresult, X[initial_samples+1:end, :])

    for i in 1:length(y)-initial_samples
        @test y_predicted[i] == y[initial_samples+i]
    end

    @test initial_samples == length(fitresult)
end

@testset "Prediction - DataFrame" begin
    initial_samples = 4
    X = [5.2603 0.40807;
    0.51512 -0.7065;
    4.721 4.5322;
    0.36285 1.0328;
    0.89978 1.6674;
    2.7659 0.74216;
    1.4156 2.2328;
    1.0293 0.12319;
    3.3528 3.261] |> DataFrame

    y = MClassification.categorical([2, 1, 2, 1, 1, 1, 1, 1, 2])
    model = MClassification.MClassifier()

    fitresult, cache, report = MClassification.fit(model, 0, X[1:initial_samples, :], y[1:initial_samples])

    y_predicted = MClassification.update_predict(model, fitresult, X[initial_samples+1:end, :])

    for i in 1:length(y)-initial_samples
        @test y_predicted[i] == y[initial_samples+i]
    end

    @test initial_samples == length(fitresult)
end

@testset "Append MC" begin
    samples = [[3, 3, 3], [2, 2, 2]]
    mc = MClassification.MCluster.MC([1, 1, 1], 1)
    centroid = [1, 1, 1]
    ls = [1, 1, 1]
    ss = [1, 1, 1]
    n = 1
    @test mc.centroid == centroid
    @test mc.ls == ls
    @test mc.ss == ss
    @test mc.n == n

    for sample in samples
        r_predicted = MClassification.MCluster.predict_r(mc, sample)

        MClassification.MCluster.append!(mc, sample)
        n += 1
        ls = ls .+ sample
        ss = ss + sample.*sample
        centroid = ls ./ n

        @test r_predicted == abs(sum((ss / n) - (ls / n).^ 2)) ^ (1/2)
        @test mc.centroid == centroid
        @test mc.ls == ls
        @test mc.ss == ss
        @test mc.n == n
    end
end

@testset "Evaluate - MLJ" begin
    # Evaluating model using evaluation of MLJ Library

    X        = [5.2603 0.40807; 0.51512 -0.7065; 4.721 4.5322;
                0.36285 1.0328; 0.89978 1.6674; 2.7659 0.74216;
                1.4156 2.2328; 1.0293 0.12319; 3.3528 3.261]
    y        = MClassification.categorical([2, 1, 2, 1, 1, 1, 1, 1, 2])

    train    = 1:4
    test     = 5:MClassification.nrows(X)

    model    = MClassification.MClassifier(r_limit=0.1)
    model    = MClassification.machine(model, X, y)

    @test MClassification.evaluate!(model, resampling=[(train, test)], measure=MClassification.accuracy)[:measurement][1] == 1.

end


text:
using MixedModelsExtras
using Test

@testset "ICC" begin
    include("icc.jl")
end


text:
struct Prior <: AbstractPrior
    G::DiscreteDistribution
    xi::ContinuousDistribution
    beta::UnivariateDistribution
    beta_mv::MultivariateDistribution

    mu::UnivariateDistribution
    mu_mv::MultivariateDistribution

    S::UnivariateDistribution
    Sigma::ContinuousDistribution
    sigma::UnivariateDistribution

    bl_shape::UnivariateDistribution
    bl_scale::UnivariateDistribution
    delta::UnivariateDistribution
    delta_mv::MultivariateDistribution
end

struct PriorSpline <: AbstractPrior
    G::DiscreteDistribution
    xi::ContinuousDistribution
    beta::UnivariateDistribution
    beta_mv::MultivariateDistribution

    mu::UnivariateDistribution
    mu_mv::MultivariateDistribution

    S::UnivariateDistribution
    Sigma::ContinuousDistribution
    sigma::UnivariateDistribution

    w::UnivariateDistribution
    delta::UnivariateDistribution
    delta_mv::MultivariateDistribution
end

Prior(data::Data) = Prior(size(data.Xp,2), size(data.Z,3), size(data.Xl,3), size(data.Xe,2))

Prior(p=4, q=3, l=9, r=3) = Prior(Poisson(3), ## G
                                  ## xi (center pi_g at 1/G and away from 0 and 1 (Garrett/Zeger [2000], Elliot [2005], Neelono [2011])
                                  Normal(0, 3/2),
                                  Normal(0, sqrt(10)), ## beta_uni
                                  MvNormal(diagm(0 => 10*ones(Float64, l))), ## beta_mv
                                  Normal(0, sqrt(10)), ## mu
                                  MvNormal(diagm(0 => 10*ones(Float64, q))), ## mu_mv
                                  Cauchy(0, 2.5), ## S
                                  InverseWishart(3, diagm(0 => ones(Float64, q))), ## C
                                  Cauchy(0, 2.5), ## sigma
                                  Gamma(1,1), ##Uniform(0.01, 100), ## bl_shape
                                  Gamma(1,1), ##Uniform(0.01, 100), ## bl_scale
                                  Normal(0, sqrt(10)), ## delta
                                  MvNormal(diagm(0 => 10*ones(Float64, r)))) ## delta_mv

PriorSpline(data::Data) = Prior(size(data.Xp,2), size(data.Z,3), size(data.Xl,3), size(data.Xe,2), size(data.haz_basis,2))

PriorSpline(p=4, q=3, l=9, r=3, k=9) = PriorSpline(Poisson(3), ## G
                                                   Normal(0, 3/2),
                                                   Normal(0, sqrt(10)), ## beta_uni
                                                   MvNormal(diagm(0 => 10*ones(Float64, l))), ## beta_mv
                                                   Normal(0, sqrt(10)), ## mu
                                                   MvNormal(diagm(0 => 10*ones(Float64, q))), ## mu_mv
                                                   Cauchy(0, 2.5), ## S
                                                   InverseWishart(3, diagm(0 => ones(Float64, q))), ## C
                                                   Cauchy(0, 2.5), ## sigma
                                                   Normal(0, sqrt(50)), ## w
                                                   Normal(0, sqrt(10)), ## delta
                                                   MvNormal(diagm(0 => 10*ones(Float64, r)))) ## delta_mv


text:
include("config.jl")
Pkg.add(PackageSpec(name="CUTEst",rev="main")) # will be removed once the new CUTEst version is released

@everywhere using CUTEst

if SOLVER == "master" || SOLVER == "current"
    @everywhere begin
        using MadNLP, MadNLPHSL
        solver = nlp -> madnlp(nlp,linear_solver=MadNLPMa57,max_wall_time=900., print_level=PRINT_LEVEL)
        function get_status(code::MadNLP.Status)
            if code == MadNLP.SOLVE_SUCCEEDED
                return 1
            elseif code == MadNLP.SOLVED_TO_ACCEPTABLE_LEVEL
                return 2
            else
                return 3
            end
        end
    end
elseif SOLVER == "ipopt"
    @everywhere begin
        solver = nlp -> ipopt(nlp,linear_solver="ma57",max_cpu_time=900., print_level=PRINT_LEVEL)
        using NLPModelsIpopt
        function get_status(code::Symbol)
            if code == :first_order
                return 1
            elseif code == :acceptable
                return 2
            else
                return 3
            end
        end
    end
elseif SOLVER == "knitro"
    # TODO
else
    error("Proper SOLVER should be given")
end

@everywhere function decodemodel(name)
    println("Decoding $name")
    finalize(CUTEstModel(name))
end

@everywhere function evalmodel(name,solver;gcoff=false)
    println("Solving $name")
    nlp = CUTEstModel(name; decode=false)
    try
        gcoff && GC.enable(false);
        mem = @allocated begin
            t = @elapsed begin
                retval = solver(nlp)
            end
        end
        gcoff && GC.enable(true);
        finalize(nlp)
        return (status=get_status(retval.status),time=t,mem=mem,iter=retval.iter)
    catch e
        finalize(nlp)
        throw(e)
    end
end

function benchmark(solver,probs;warm_up_probs = [], decode = false)
    println("Warming up (forcing JIT compile)")
    decode && broadcast(decodemodel,warm_up_probs)
    r = [remotecall.(prob->evalmodel(prob,solver;gcoff=GCOFF),i,warm_up_probs) for i in procs() if i!= 1]
    fetch.(r)

    println("Decoding problems")
    decode && broadcast(decodemodel,probs)

    println("Solving problems")
    retvals = pmap(prob->evalmodel(prob,solver;gcoff=GCOFF),probs)
    status = [retval.status for retval in retvals]
    time   = [retval.time for retval in retvals]
    mem    = [retval.mem for retval in retvals]
    iter   = [retval.iter for retval in retvals]
    status,time,mem,iter
end

exclude = [
    # MadNLP running into error
    # Ipopt running into error
    "EG3", # lfact blows up
    # Problems that are hopelessly large
    "TAX213322","TAXR213322","TAX53322","TAXR53322",
    "YATP1LS","YATP2LS","YATP1CLS","YATP2CLS",
    "CYCLOOCT","CYCLOOCF",
    "LIPPERT1",
    "GAUSSELM",
    "BA-L52LS","BA-L73LS","BA-L21LS"
]

if QUICK
    probs = readdlm("cutest-quick-names.csv")[:]
else
    probs = CUTEst.select()
end

filter!(e->!(e in exclude),probs)

status,time,mem,iter = benchmark(solver,probs;warm_up_probs = ["EIGMINA"], decode = DECODE)

f = !QUICK ? "" : "-quick"
writedlm("name-cutest$f.csv",probs,',')
writedlm("status-cutest-$(SOLVER)$f.csv",status,',')
writedlm("time-cutest-$(SOLVER)$f.csv",time,',')
writedlm("mem-cutest-$(SOLVER)$f.csv",mem,',')
writedlm("iter-cutest-$(SOLVER)$f.csv",iter,',')


text:
module EulerGravity
  export EulerGravityLaw, γ, grav

  import ..Atum
  using ..Atum: avg, logavg, roe_avg
  using StaticArrays
  using LinearAlgebra: I

  struct EulerGravityLaw{γ, grav, FT, D, S} <: Atum.AbstractBalanceLaw{FT, D, S}
    function EulerGravityLaw{FT, D}(; γ = 7 // 5, grav = 981 // 100) where {FT, D}
      S = 2 + D
      new{FT(γ), FT(grav), FT, D, S}()
    end
  end

  γ(::EulerGravityLaw{_γ}) where {_γ} = _γ
  grav(::EulerGravityLaw{_γ, _grav}) where {_γ, _grav} = _grav

  function geopotential(law, x⃗)
    z = last(x⃗)
    grav(law) * z
  end

  function varsindices(law::EulerGravityLaw)
    S = Atum.numberofstates(law)
    ix_ρ = 1
    ix_ρu⃗ = StaticArrays.SUnitRange(2, S - 1)
    ix_ρe = S
    return ix_ρ, ix_ρu⃗, ix_ρe
  end

  function unpackstate(law::EulerGravityLaw, q)
    ix_ρ, ix_ρu⃗, ix_ρe = varsindices(law)
    @inbounds q[ix_ρ], q[ix_ρu⃗], q[ix_ρe]
  end

  function pressure(law::EulerGravityLaw, ρ, ρu⃗, ρe, Φ)
    (γ(law) - 1) * (ρe - ρu⃗' * ρu⃗ / 2ρ - ρ * Φ)
  end
  function energy(law::EulerGravityLaw, ρ, ρu⃗, p, Φ)
    p / (γ(law) - 1) + ρu⃗' * ρu⃗ / 2ρ + ρ * Φ
  end
  function soundspeed(law::EulerGravityLaw, ρ, p)
    sqrt(γ(law) * p / ρ)
  end
  function soundspeed(law::EulerGravityLaw, ρ, ρu⃗, ρe, Φ)
    soundspeed(law, ρ, pressure(law, ρ, ρu⃗, ρe, Φ))
  end

  function Atum.flux(law::EulerGravityLaw, q, x⃗)
    ρ, ρu⃗, ρe = unpackstate(law, q)

    Φ = geopotential(law, x⃗)

    u⃗ = ρu⃗ / ρ
    p = pressure(law, ρ, ρu⃗, ρe, Φ)

    fρ = ρu⃗
    fρu⃗ = ρu⃗ * u⃗' + p * I
    fρe = u⃗ * (ρe + p)

    hcat(fρ, fρu⃗, fρe)
  end

  function Atum.nonconservative_term!(law::EulerGravityLaw, dq, q, x⃗)
    ix_ρ, ix_ρu⃗, _ = varsindices(law)
    @inbounds dq[ix_ρu⃗[end]] -= q[ix_ρ] * grav(law)
  end

  function Atum.wavespeed(law::EulerGravityLaw, n⃗, q, x⃗)
    ρ, ρu⃗, ρe = unpackstate(law, q)

    Φ = geopotential(law, x⃗)

    u⃗ = ρu⃗ / ρ
    abs(n⃗' * u⃗) + soundspeed(law, ρ, ρu⃗, ρe, Φ)
  end

  function Atum.surfaceflux(::Atum.RoeFlux, law::EulerGravityLaw, n⃗, x⃗, q⁻, q⁺)
    Φ = geopotential(law, x⃗)

    f⁻ = Atum.flux(law, q⁻, x⃗)
    f⁺ = Atum.flux(law, q⁺, x⃗)

    ρ⁻, ρu⃗⁻, ρe⁻ = unpackstate(law, q⁻)
    u⃗⁻ = ρu⃗⁻ / ρ⁻
    e⁻ = ρe⁻ / ρ⁻
    p⁻ = pressure(law, ρ⁻, ρu⃗⁻, ρe⁻, Φ)
    h⁻ = e⁻ + p⁻ / ρ⁻
    c⁻ = soundspeed(law, ρ⁻, p⁻)

    ρ⁺, ρu⃗⁺, ρe⁺ = unpackstate(law, q⁺)
    u⃗⁺ = ρu⃗⁺ / ρ⁺
    e⁺ = ρe⁺ / ρ⁺
    p⁺ = pressure(law, ρ⁺, ρu⃗⁺, ρe⁺, Φ)
    h⁺ = e⁺ + p⁺ / ρ⁺
    c⁺ = soundspeed(law, ρ⁺, p⁺)

    ρ = sqrt(ρ⁻ * ρ⁺)
    u⃗ = roe_avg(ρ⁻, ρ⁺, u⃗⁻, u⃗⁺)
    h = roe_avg(ρ⁻, ρ⁺, h⁻, h⁺)
    c = roe_avg(ρ⁻, ρ⁺, c⁻, c⁺)

    uₙ = u⃗' * n⃗

    Δρ = ρ⁺ - ρ⁻
    Δp = p⁺ - p⁻
    Δu⃗ = u⃗⁺ - u⃗⁻
    Δuₙ = Δu⃗' * n⃗

    c⁻² = 1 / c^2
    w1 = abs(uₙ - c) * (Δp - ρ * c * Δuₙ) * c⁻² / 2
    w2 = abs(uₙ + c) * (Δp + ρ * c * Δuₙ) * c⁻² / 2
    w3 = abs(uₙ) * (Δρ - Δp * c⁻²)
    w4 = abs(uₙ) * ρ

    fp_ρ = (w1 + w2 + w3) / 2
    fp_ρu⃗ = (w1 * (u⃗ - c * n⃗) +
             w2 * (u⃗ + c * n⃗) +
             w3 * u⃗ +
             w4 * (Δu⃗ - Δuₙ * n⃗)) / 2
    fp_ρe = (w1 * (h - c * uₙ) +
             w2 * (h + c * uₙ) +
             w3 * (u⃗' * u⃗ / 2 + Φ) +
             w4 * (u⃗' * Δu⃗ - uₙ * Δuₙ)) / 2

    (f⁻ + f⁺)' * n⃗ / 2 - SVector(fp_ρ, fp_ρu⃗..., fp_ρe)
  end

  function Atum.twopointflux(::Atum.EntropyConservativeFlux,
                             law::EulerGravityLaw,
                             q₁, x⃗₁, q₂, x⃗₂)
      FT = eltype(law)
      ρ₁, ρu⃗₁, ρe₁ = unpackstate(law, q₁)
      ρ₂, ρu⃗₂, ρe₂ = unpackstate(law, q₂)

      Φ₁ = geopotential(law, x⃗₁)
      u⃗₁ = ρu⃗₁ / ρ₁
      p₁ = pressure(law, ρ₁, ρu⃗₁, ρe₁, Φ₁)
      b₁ = ρ₁ / 2p₁

      Φ₂ = geopotential(law, x⃗₂)
      u⃗₂ = ρu⃗₂ / ρ₂
      p₂ = pressure(law, ρ₂, ρu⃗₂, ρe₂, Φ₂)
      b₂ = ρ₂ / 2p₂

      ρ_avg = avg(ρ₁, ρ₂)
      u⃗_avg = avg(u⃗₁, u⃗₂)
      b_avg = avg(b₁, b₂)
      Φ_avg = avg(Φ₁, Φ₂)

      u²_avg = avg(u⃗₁' * u⃗₁, u⃗₂' * u⃗₂)
      ρ_log = logavg(ρ₁, ρ₂)
      b_log = logavg(b₁, b₂)

      fρ = u⃗_avg * ρ_log
      fρu⃗ = u⃗_avg * fρ' + ρ_avg / 2b_avg * I
      fρe = (1 / (2 * (γ(law) - 1) * b_log) - u²_avg / 2 + Φ_avg) * fρ + fρu⃗ * u⃗_avg

      # fluctuation
      α = b_avg * ρ_log / 2b₁
      fρu⃗ -= α * (Φ₁ - Φ₂) * I

      hcat(fρ, fρu⃗, fρe)
  end
end