text:
"""
    stack(df::AbstractDataFrame[, measure_vars[, id_vars] ];
          variable_name=:variable, value_name=:value,
          view::Bool=false, variable_eltype::Type=String)

Stack a data frame `df`, i.e. convert it from wide to long format.

Return the long-format `DataFrame` with: columns for each of the `id_vars`,
column `value_name` (`:value` by default)
holding the values of the stacked columns (`measure_vars`), and
column `variable_name` (`:variable` by default) a vector holding
the name of the corresponding `measure_vars` variable.

If `view=true` then return a stacked view of a data frame (long format).
The result is a view because the columns are special `AbstractVectors`
that return views into the original data frame.

# Arguments
- `df` : the AbstractDataFrame to be stacked
- `measure_vars` : the columns to be stacked (the measurement variables),
  as a column selector ($COLUMNINDEX_STR; $MULTICOLUMNINDEX_STR).
  If neither `measure_vars` or `id_vars` are given, `measure_vars`
  defaults to all floating point columns.
- `id_vars` : the identifier columns that are repeated during stacking,
  as a column selector ($COLUMNINDEX_STR; $MULTICOLUMNINDEX_STR).
  Defaults to all variables that are not `measure_vars`
- `variable_name` : the name (`Symbol` or string) of the new stacked column that
  shall hold the names of each of `measure_vars`
- `value_name` : the name (`Symbol` or string) of the new stacked column containing
  the values from each of `measure_vars`
- `view` : whether the stacked data frame should be a view rather than contain
  freshly allocated vectors.
- `variable_eltype` : determines the element type of column `variable_name`.
  By default a `PooledArray{String}` is created.
  If `variable_eltype=Symbol` a `PooledVector{Symbol}` is created,
  and if `variable_eltype=CategoricalValue{String}`
  a `CategoricalArray{String}` is produced (call `using CategoricalArrays` first if needed)
  Passing any other type `T` will produce a `PooledVector{T}` column
  as long as it supports conversion from `String`.
  When `view=true`, a `RepeatedVector{T}` is produced.

# Examples
```jldoctest
julia> df = DataFrame(a=repeat(1:3, inner=2),
                      b=repeat(1:2, inner=3),
                      c=repeat(1:1, inner=6),
                      d=repeat(1:6, inner=1),
                      e=string.('a':'f'))
6×5 DataFrame
 Row │ a      b      c      d      e
     │ Int64  Int64  Int64  Int64  String
─────┼────────────────────────────────────
   1 │     1      1      1      1  a
   2 │     1      1      1      2  b
   3 │     2      1      1      3  c
   4 │     2      2      1      4  d
   5 │     3      2      1      5  e
   6 │     3      2      1      6  f

julia> stack(df, [:c, :d])
12×5 DataFrame
 Row │ a      b      e       variable  value
     │ Int64  Int64  String  String    Int64
─────┼───────────────────────────────────────
   1 │     1      1  a       c             1
   2 │     1      1  b       c             1
   3 │     2      1  c       c             1
   4 │     2      2  d       c             1
   5 │     3      2  e       c             1
   6 │     3      2  f       c             1
   7 │     1      1  a       d             1
   8 │     1      1  b       d             2
   9 │     2      1  c       d             3
  10 │     2      2  d       d             4
  11 │     3      2  e       d             5
  12 │     3      2  f       d             6

julia> stack(df, [:c, :d], [:a])
12×3 DataFrame
 Row │ a      variable  value
     │ Int64  String    Int64
─────┼────────────────────────
   1 │     1  c             1
   2 │     1  c             1
   3 │     2  c             1
   4 │     2  c             1
   5 │     3  c             1
   6 │     3  c             1
   7 │     1  d             1
   8 │     1  d             2
   9 │     2  d             3
  10 │     2  d             4
  11 │     3  d             5
  12 │     3  d             6

julia> stack(df, Not([:a, :b, :e]))
12×5 DataFrame
 Row │ a      b      e       variable  value
     │ Int64  Int64  String  String    Int64
─────┼───────────────────────────────────────
   1 │     1      1  a       c             1
   2 │     1      1  b       c             1
   3 │     2      1  c       c             1
   4 │     2      2  d       c             1
   5 │     3      2  e       c             1
   6 │     3      2  f       c             1
   7 │     1      1  a       d             1
   8 │     1      1  b       d             2
   9 │     2      1  c       d             3
  10 │     2      2  d       d             4
  11 │     3      2  e       d             5
  12 │     3      2  f       d             6

julia> stack(df, Not([:a, :b, :e]), variable_name=:somemeasure)
12×5 DataFrame
 Row │ a      b      e       somemeasure  value
     │ Int64  Int64  String  String       Int64
─────┼──────────────────────────────────────────
   1 │     1      1  a       c                1
   2 │     1      1  b       c                1
   3 │     2      1  c       c                1
   4 │     2      2  d       c                1
   5 │     3      2  e       c                1
   6 │     3      2  f       c                1
   7 │     1      1  a       d                1
   8 │     1      1  b       d                2
   9 │     2      1  c       d                3
  10 │     2      2  d       d                4
  11 │     3      2  e       d                5
  12 │     3      2  f       d                6
```
"""
function stack(df::AbstractDataFrame,
               measure_vars = findall(col -> eltype(col) <: Union{AbstractFloat, Missing},
                                      eachcol(df)),
               id_vars = Not(measure_vars);
               variable_name::SymbolOrString=:variable,
               value_name::SymbolOrString=:value, view::Bool=false,
               variable_eltype::Type=String)
    variable_name_s = Symbol(variable_name)
    value_name_s = Symbol(value_name)
    # getindex from index returns either Int or AbstractVector{Int}
    mv_tmp = index(df)[measure_vars]
    ints_measure_vars = mv_tmp isa Int ? [mv_tmp] : mv_tmp
    idv_tmp = index(df)[id_vars]
    ints_id_vars = idv_tmp isa Int ? [idv_tmp] : idv_tmp
    if view
        return _stackview(df, ints_measure_vars, ints_id_vars,
                          variable_name=variable_name_s,
                          value_name=value_name_s,
                          variable_eltype=variable_eltype)
    end
    N = length(ints_measure_vars)
    cnames = _names(df)[ints_id_vars]
    push!(cnames, variable_name_s)
    push!(cnames, value_name_s)
    if variable_eltype === Symbol
        catnms = PooledArray(_names(df)[ints_measure_vars])
    elseif variable_eltype === String
        catnms = PooledArray(names(df, ints_measure_vars))
    else
        # this covers CategoricalArray{String} in particular
        # (note that copyto! inserts levels in their order of appearance)
        nms = names(df, ints_measure_vars)
        simnms = similar(nms, variable_eltype)
        catnms = simnms isa Vector ? PooledArray(catnms) : simnms
        copyto!(catnms, nms)
    end
    return DataFrame(AbstractVector[[repeat(df[!, c], outer=N) for c in ints_id_vars]..., # id_var columns
                                    repeat(catnms, inner=nrow(df)),                       # variable
                                    vcat([df[!, c] for c in ints_measure_vars]...)],      # value
                     cnames, copycols=false)
end

function _stackview(df::AbstractDataFrame, measure_vars::AbstractVector{Int},
                    id_vars::AbstractVector{Int}; variable_name::Symbol,
                    value_name::Symbol, variable_eltype::Type)
    N = length(measure_vars)
    cnames = _names(df)[id_vars]
    push!(cnames, variable_name)
    push!(cnames, value_name)
    if variable_eltype === Symbol
        catnms = _names(df)[measure_vars]
    elseif variable_eltype === String
        catnms = names(df, measure_vars)
    else
        # this covers CategoricalArray{String} in particular,
        # as copyto! inserts levels in their order of appearance
        nms = names(df, measure_vars)
        catnms = copyto!(similar(nms, variable_eltype), nms)
    end
    return DataFrame(AbstractVector[[RepeatedVector(df[!, c], 1, N) for c in id_vars]..., # id_var columns
                                    RepeatedVector(catnms, nrow(df), 1),                  # variable
                                    StackedVector(Any[df[!, c] for c in measure_vars])],  # value
                     cnames, copycols=false)
end

"""
    unstack(df::AbstractDataFrame, rowkeys, colkey, value;
            renamecols::Function=identity, allowmissing::Bool=false,
            allowduplicates::Bool=false, valuestransform=nothing,
            fill=missing)
    unstack(df::AbstractDataFrame, colkey, value;
            renamecols::Function=identity, allowmissing::Bool=false,
            allowduplicates::Bool=false, valuestransform=nothing,
            fill=missing)
    unstack(df::AbstractDataFrame;
            renamecols::Function=identity, allowmissing::Bool=false,
            allowduplicates::Bool=false, valuestransform=nothing,
            fill=missing)

Unstack data frame `df`, i.e. convert it from long to wide format.

Row and column keys will be ordered in the order of their first appearance.

# Positional arguments
- `df` : the AbstractDataFrame to be unstacked
- `rowkeys` : the columns with a unique key for each row, if not given, find a
  key by grouping on anything not a `colkey` or `value`. Can be any column
  selector ($COLUMNINDEX_STR; $MULTICOLUMNINDEX_STR). If `rowkeys` contains no
  columns all rows are assumed to have the same key.
- `colkey` : the column ($COLUMNINDEX_STR) holding the column names in wide
  format, defaults to `:variable`
- `values` : the column storing values ($COLUMNINDEX_STR), defaults to `:value`

# Keyword arguments

- `renamecols`: a function called on each unique value in `colkey`; it must
  return the name of the column to be created (typically as a string or a
  `Symbol`). Duplicates in resulting names when converted to `Symbol` are not
  allowed. By default no transformation is performed.
- `allowmissing`: if `false` (the default) then an error will be thrown if
  `colkey` contains `missing` values; if `true` then a column referring to
  `missing` value will be created.
- `allowduplicates`: if `false` (the default) then an error an error will be
  thrown if combination of `rowkeys` and `colkey` contains duplicate entries; if
  `true` then the last encountered `value` will be retained;
  this keyword argument is ignored if `valuestransform` keyword argument is passed.
- `valuestransform`: if passed then `allowduplicates` is ignored and instead
   the passed function will be called on a vector view containing all elements
   for each combination of `rowkeys` and `colkey` present in the data.
- `fill`: missing row/column combinations are filled with this value. The
  default is `missing`. If the `value` column is a `CategoricalVector` and
  `fill` is not `missing` then in order to keep unstacked value columns also
  `CategoricalVector` the `fill` must be passed as `CategoricalValue`

# Examples

```jldoctest
julia> wide = DataFrame(id=1:6,
                        a=repeat(1:3, inner=2),
                        b=repeat(1.0:2.0, inner=3),
                        c=repeat(1.0:1.0, inner=6),
                        d=repeat(1.0:3.0, inner=2))
6×5 DataFrame
 Row │ id     a      b        c        d
     │ Int64  Int64  Float64  Float64  Float64
─────┼─────────────────────────────────────────
   1 │     1      1      1.0      1.0      1.0
   2 │     2      1      1.0      1.0      1.0
   3 │     3      2      1.0      1.0      2.0
   4 │     4      2      2.0      1.0      2.0
   5 │     5      3      2.0      1.0      3.0
   6 │     6      3      2.0      1.0      3.0

julia> long = stack(wide)
18×4 DataFrame
 Row │ id     a      variable  value
     │ Int64  Int64  String    Float64
─────┼─────────────────────────────────
   1 │     1      1  b             1.0
   2 │     2      1  b             1.0
   3 │     3      2  b             1.0
   4 │     4      2  b             2.0
   5 │     5      3  b             2.0
   6 │     6      3  b             2.0
   7 │     1      1  c             1.0
   8 │     2      1  c             1.0
  ⋮  │   ⋮      ⋮       ⋮         ⋮
  12 │     6      3  c             1.0
  13 │     1      1  d             1.0
  14 │     2      1  d             1.0
  15 │     3      2  d             2.0
  16 │     4      2  d             2.0
  17 │     5      3  d             3.0
  18 │     6      3  d             3.0
                         3 rows omitted

julia> unstack(long)
6×5 DataFrame
 Row │ id     a      b         c         d
     │ Int64  Int64  Float64?  Float64?  Float64?
─────┼────────────────────────────────────────────
   1 │     1      1       1.0       1.0       1.0
   2 │     2      1       1.0       1.0       1.0
   3 │     3      2       1.0       1.0       2.0
   4 │     4      2       2.0       1.0       2.0
   5 │     5      3       2.0       1.0       3.0
   6 │     6      3       2.0       1.0       3.0

julia> unstack(long, :variable, :value)
6×5 DataFrame
 Row │ id     a      b         c         d
     │ Int64  Int64  Float64?  Float64?  Float64?
─────┼────────────────────────────────────────────
   1 │     1      1       1.0       1.0       1.0
   2 │     2      1       1.0       1.0       1.0
   3 │     3      2       1.0       1.0       2.0
   4 │     4      2       2.0       1.0       2.0
   5 │     5      3       2.0       1.0       3.0
   6 │     6      3       2.0       1.0       3.0

julia> unstack(long, :id, :variable, :value)
6×4 DataFrame
 Row │ id     b         c         d
     │ Int64  Float64?  Float64?  Float64?
─────┼─────────────────────────────────────
   1 │     1       1.0       1.0       1.0
   2 │     2       1.0       1.0       1.0
   3 │     3       1.0       1.0       2.0
   4 │     4       2.0       1.0       2.0
   5 │     5       2.0       1.0       3.0
   6 │     6       2.0       1.0       3.0

julia> unstack(long, [:id, :a], :variable, :value)
6×5 DataFrame
 Row │ id     a      b         c         d
     │ Int64  Int64  Float64?  Float64?  Float64?
─────┼────────────────────────────────────────────
   1 │     1      1       1.0       1.0       1.0
   2 │     2      1       1.0       1.0       1.0
   3 │     3      2       1.0       1.0       2.0
   4 │     4      2       2.0       1.0       2.0
   5 │     5      3       2.0       1.0       3.0
   6 │     6      3       2.0       1.0       3.0

julia> unstack(long, :id, :variable, :value, renamecols=x->Symbol(:_, x))
6×4 DataFrame
 Row │ id     _b        _c        _d
     │ Int64  Float64?  Float64?  Float64?
─────┼─────────────────────────────────────
   1 │     1       1.0       1.0       1.0
   2 │     2       1.0       1.0       1.0
   3 │     3       1.0       1.0       2.0
   4 │     4       2.0       1.0       2.0
   5 │     5       2.0       1.0       3.0
   6 │     6       2.0       1.0       3.0
```
Note that there are some differences between the widened results above.

```jldoctest
julia> df = DataFrame(id=["1", "1", "2"],
                      variable=["Var1", "Var2", "Var1"],
                      value=[1, 2, 3])
3×3 DataFrame
 Row │ id      variable  value
     │ String  String    Int64
─────┼─────────────────────────
   1 │ 1       Var1          1
   2 │ 1       Var2          2
   3 │ 2       Var1          3

julia> unstack(df, :variable, :value, fill=0)
2×3 DataFrame
 Row │ id      Var1   Var2
     │ String  Int64  Int64
─────┼──────────────────────
   1 │ 1           1      2
   2 │ 2           3      0

julia> df = DataFrame(cols=["a", "a", "b"], values=[1, 2, 4])
3×2 DataFrame
 Row │ cols    values
     │ String  Int64
─────┼────────────────
   1 │ a            1
   2 │ a            2
   3 │ b            4

julia> unstack(df, :cols, :values, valuestransform=copy)
1×2 DataFrame
 Row │ a        b
     │ Array…?  Array…?
─────┼──────────────────
   1 │ [1, 2]   [4]

julia> unstack(df, :cols, :values, valuestransform=sum)
1×2 DataFrame
 Row │ a       b
     │ Int64?  Int64?
─────┼────────────────
   1 │      3       4
```
"""
function unstack(df::AbstractDataFrame, rowkeys, colkey::ColumnIndex,
                 values::ColumnIndex; renamecols::Function=identity,
                 allowmissing::Bool=false, allowduplicates::Bool=false,
                 valuestransform=nothing, fill=missing)
    if !isnothing(valuestransform)
        groupcols = vcat(index(df)[rowkeys], index(df)[colkey])
        @assert groupcols isa AbstractVector{Int}
        gdf = groupby(df, groupcols)
        if check_aggregate(valuestransform, df[!, values]) isa AbstractAggregate
            # if valuestransform function is AbstractAggregate
            # then we are sure it will return a scalar number so we can
            # leave it as is and be sure we use fast path in combine
            agg_fun = valuestransform
        else
            # in general valuestransform function could return e.g. a vector,
            # which would get expanded to multiple rows so we protect it with
            # Ref that will get unwrapped by combine
            agg_fun = Ref∘valuestransform
        end
        df_op = combine(gdf, values => agg_fun, renamecols=false)
        group_rows = find_group_row(gdf)
        if !issorted(group_rows)
            df_op = df_op[sortperm(group_rows), :]
        end
        # set allowduplicates to true as we should not have any duplicates now
        # and allowduplicates=true is a bit faster
        allowduplicates = true
    else
        df_op = df
    end
    # use df_op below to make sure it is type stable
    rowkey_ints = vcat(index(df_op)[rowkeys])
    @assert rowkey_ints isa AbstractVector{Int}
    g_rowkey = groupby(df_op, rowkey_ints)
    g_colkey = groupby(df_op, colkey)
    valuecol = df_op[!, values]
    return _unstack(df_op, rowkey_ints, index(df_op)[colkey], g_colkey,
                    valuecol, g_rowkey, renamecols, allowmissing, allowduplicates, fill)
end

function unstack(df::AbstractDataFrame, colkey::ColumnIndex, values::ColumnIndex;
                 renamecols::Function=identity,
                 allowmissing::Bool=false, allowduplicates::Bool=false,
                 valuestransform=nothing, fill=missing)
    colkey_int = index(df)[colkey]
    value_int = index(df)[values]
    return unstack(df, Not(colkey_int, value_int), colkey_int, value_int,
            renamecols=renamecols, allowmissing=allowmissing,
            allowduplicates=allowduplicates, valuestransform=valuestransform, fill=fill)
end

unstack(df::AbstractDataFrame; renamecols::Function=identity,
        allowmissing::Bool=false, allowduplicates::Bool=false,
        valuestransform=nothing, fill=missing) =
    unstack(df, :variable, :value, renamecols=renamecols, allowmissing=allowmissing,
            allowduplicates=allowduplicates, valuestransform=valuestransform, fill=fill)

# we take into account the fact that idx, starts and ends are computed lazily
# so we rather directly reference the gdf.groups
# this function is tailor made for unstack so it does assume that no groups were
# dropped (i.e. gdf.groups does not contain 0 entries)
function find_group_row(gdf::GroupedDataFrame)
    rows = zeros(Int, length(gdf))
    isempty(rows) && return rows

    filled = 0
    i = 1
    groups = gdf.groups
    while filled < length(gdf)
        group = groups[i]
        if rows[group] == 0
            rows[group] = i
            filled += 1
        end
        i += 1
    end
    return rows # return row index of first occurrence of each group in gdf.groups
end

function _unstack(df::AbstractDataFrame, rowkeys::AbstractVector{Int},
                  colkey::Int, g_colkey::GroupedDataFrame,
                  valuecol::AbstractVector, g_rowkey::GroupedDataFrame,
                  renamecols::Function, allowmissing::Bool,
                  allowduplicates::Bool, fill)
    rowref = g_rowkey.groups
    row_group_row_idxs = find_group_row(g_rowkey)
    Nrow = length(g_rowkey)

    @assert groupcols(g_colkey) == _names(df)[colkey:colkey]
    colref = g_colkey.groups
    Ncol = length(g_colkey)
    col_group_row_idxs = find_group_row(g_colkey)
    colref_map = df[col_group_row_idxs, colkey]
    if any(ismissing, colref_map) && !allowmissing
        throw(ArgumentError("Missing value in variable :$(_names(df)[colkey]). " *
                            "Pass `allowmissing=true` to skip missings."))
    end
    @assert length(rowref) == length(colref) == length(valuecol)

    unstacked_val = [fill!(similar(valuecol,
                           promote_type(eltype(valuecol), typeof(fill)),
                           Nrow),
                     fill) for _ in 1:Ncol]

    # use a separate path for allowduplicates to reduce memory use and increase speed
    if allowduplicates
        for (k, (row_id, col_id, val)) in enumerate(zip(rowref, colref, valuecol))
            unstacked_val[col_id][row_id] = val
        end
    else
        mask_filled = falses(Nrow, Ncol)
        for (k, (row_id, col_id, val)) in enumerate(zip(rowref, colref, valuecol))
            if mask_filled[row_id, col_id]
                bad_key = tuple((df[k, s] for s in rowkeys)...)
                bad_var = colref_map[col_id]
                throw(ArgumentError("Duplicate entries in unstack at row $k for key "*
                                    "$bad_key and variable $bad_var. " *
                                    "Pass allowduplicates=true to allow them."))
            end
            unstacked_val[col_id][row_id] = val
            mask_filled[row_id, col_id] = true
        end
    end

    # note that Symbol(renamecols(x)) must produce unique column names
    # and names between df1 and df2 must be unique
    df1 = df[row_group_row_idxs, g_rowkey.cols]
    df2 = DataFrame(unstacked_val, Symbol[Symbol(renamecols(x)) for x in colref_map],
                    copycols=false)

    @assert length(col_group_row_idxs) == ncol(df2)
    # avoid reordering when col_group_row_idxs was already ordered
    if !issorted(col_group_row_idxs)
        df2 = df2[!, sortperm(col_group_row_idxs)]
    end

    res_df = hcat(df1, df2, copycols=false)

    @assert length(row_group_row_idxs) == nrow(res_df)
    # avoid reordering when row_group_row_idxs was already ordered
    if !issorted(row_group_row_idxs)
        res_df = res_df[sortperm(row_group_row_idxs), :]
    end

    return res_df
end

"""
    StackedVector <: AbstractVector

An `AbstractVector` that is a linear, concatenated view into
another set of AbstractVectors

NOTE: Not exported.

# Constructor
```julia
StackedVector(d::AbstractVector)
```

# Arguments
- `d...` : one or more AbstractVectors

# Examples
```julia
StackedVector(Any[[1, 2], [9, 10], [11, 12]])  # [1, 2, 9, 10, 11, 12]
```
"""
struct StackedVector{T} <: AbstractVector{T}
    components::Vector{Any}
end

StackedVector(d::AbstractVector) =
    StackedVector{promote_type(map(eltype, d)...)}(d)

function Base.getindex(v::StackedVector{T}, i::Int)::T where T
    lengths = [length(x)::Int for x in v.components]
    cumlengths = [0; cumsum(lengths)]
    j = searchsortedlast(cumlengths .+ 1, i)
    if j > length(cumlengths)
        error("indexing bounds error")
    end
    k = i - cumlengths[j]
    if k < 1 || k > length(v.components[j])
        error("indexing bounds error")
    end
    return v.components[j][k]
end

Base.IndexStyle(::Type{StackedVector}) = Base.IndexLinear()
Base.size(v::StackedVector) = (length(v),)
Base.length(v::StackedVector) = sum(map(length, v.components))
Base.eltype(v::Type{StackedVector{T}}) where {T} = T
Base.similar(v::StackedVector, T::Type, dims::Union{Integer, AbstractUnitRange}...) =
    similar(v.components[1], T, dims...)

"""
    RepeatedVector{T} <: AbstractVector{T}

An AbstractVector that is a view into another AbstractVector with
repeated elements

NOTE: Not exported.

# Constructor
```julia
RepeatedVector(parent::AbstractVector, inner::Int, outer::Int)
```

# Arguments
- `parent` : the AbstractVector that's repeated
- `inner` : the numer of times each element is repeated
- `outer` : the numer of times the whole vector is repeated after
  expanded by `inner`

`inner` and `outer` have the same meaning as similarly named arguments
to `repeat`.

# Examples
```julia
RepeatedVector([1, 2], 3, 1)   # [1, 1, 1, 2, 2, 2]
RepeatedVector([1, 2], 1, 3)   # [1, 2, 1, 2, 1, 2]
RepeatedVector([1, 2], 2, 2)   # [1, 2, 1, 2, 1, 2, 1, 2]
```
"""
struct RepeatedVector{T} <: AbstractVector{T}
    parent::AbstractVector{T}
    inner::Int
    outer::Int
end

Base.parent(v::RepeatedVector) = v.parent

function Base.getindex(v::RepeatedVector, i::Int)
    N = length(parent(v))
    idx = Base.fld1(mod1(i, v.inner*N), v.inner)
    parent(v)[idx]
end

Base.IndexStyle(::Type{<:RepeatedVector}) = Base.IndexLinear()
Base.size(v::RepeatedVector) = (length(v),)
Base.length(v::RepeatedVector) = v.inner * v.outer * length(parent(v))
Base.eltype(v::Type{RepeatedVector{T}}) where {T} = T
Base.reverse(v::RepeatedVector) = RepeatedVector(reverse(parent(v)), v.inner, v.outer)
Base.similar(v::RepeatedVector, T::Type, dims::Dims) = similar(parent(v), T, dims)
Base.unique(v::RepeatedVector) = unique(parent(v))

Base.transpose(::AbstractDataFrame, args...; kwargs...) =
    throw(ArgumentError("`transpose` not defined for `AbstractDataFrame`s. Try `permutedims` instead"))

"""
    permutedims(df::AbstractDataFrame, src_namescol::Union{Int, Symbol, AbstractString},
                [dest_namescol::Union{Symbol, AbstractString}];
                makeunique::Bool=false, strict::Bool=true)

Turn `df` on its side such that rows become columns
and values in the column indexed by `src_namescol` become the names of new columns.
In the resulting `DataFrame`, column names of `df` will become the first column
with name specified by `dest_namescol`.

# Arguments
- `df` : the `AbstractDataFrame`
- `src_namescol` : the column that will become the new header.
- `dest_namescol` : the name of the first column in the returned `DataFrame`.
  Defaults to the same name as `src_namescol`.
- `makeunique` : if `false` (the default), an error will be raised
  if duplicate names are found; if `true`, duplicate names will be suffixed
  with `_i` (`i` starting at 1 for the first duplicate).
- `strict` : if `true` (the default), an error will be raised if the values
  contained in the `src_namescol` are not all `Symbol` or all `AbstractString`,
  or can all be converted to `String` using `convert`. If `false`
  then any values are accepted and the will be changed to strings using
  the `string` function.

Note: The element types of columns in resulting `DataFrame`
(other than the first column, which always has element type `String`)
will depend on the element types of _all_ input columns
based on the result of `promote_type`.
That is, if the source data frame contains `Int` and `Float64` columns,
resulting columns will have element type `Float64`. If the source has
`Int` and `String` columns, resulting columns will have element type `Any`.

# Examples

```jldoctest
julia> df1 = DataFrame(a=["x", "y"], b=[1.0, 2.0], c=[3, 4], d=[true, false])
2×4 DataFrame
 Row │ a       b        c      d
     │ String  Float64  Int64  Bool
─────┼───────────────────────────────
   1 │ x           1.0      3   true
   2 │ y           2.0      4  false

julia> permutedims(df1, 1) # note the column types
3×3 DataFrame
 Row │ a       x        y
     │ String  Float64  Float64
─────┼──────────────────────────
   1 │ b           1.0      2.0
   2 │ c           3.0      4.0
   3 │ d           1.0      0.0

julia> df2 = DataFrame(a=["x", "y"], b=[1, "two"], c=[3, 4], d=[true, false])
2×4 DataFrame
 Row │ a       b    c      d
     │ String  Any  Int64  Bool
─────┼───────────────────────────
   1 │ x       1        3   true
   2 │ y       two      4  false

julia> permutedims(df2, 1, "different_name")
3×3 DataFrame
 Row │ different_name  x     y
     │ String          Any   Any
─────┼─────────────────────────────
   1 │ b               1     two
   2 │ c               3     4
   3 │ d               true  false
```
"""
function Base.permutedims(df::AbstractDataFrame, src_namescol::ColumnIndex,
                          dest_namescol::Union{Symbol, AbstractString};
                          makeunique::Bool=false, strict::Bool=true)

    if src_namescol isa Integer
        1 <= src_namescol <= ncol(df) || throw(BoundsError(index(df), src_namescol))
    end
    src_col_names = df[!, src_namescol]
    local new_col_names
    if eltype(src_col_names) <: SymbolOrString
        new_col_names = src_col_names
    elseif all(x -> x isa Symbol, src_col_names)
        new_col_names = collect(Symbol, src_col_names)
    elseif !strict
        new_col_names = string.(src_col_names)
    else
        try
            new_col_names = collect(String, src_col_names)
        catch e
            if e isa MethodError && e.f === convert
                throw(ArgumentError("all elements of src_namescol must support " *
                                    "conversion to String"))
            else
                rethrow(e)
            end
        end
    end

    df_notsrc = df[!, Not(src_namescol)]
    df_permuted = DataFrame(dest_namescol => names(df_notsrc))

    if ncol(df_notsrc) == 0
        df_tmp = DataFrame(AbstractVector[[] for _ in 1:nrow(df)], new_col_names,
                           makeunique=makeunique, copycols=false)
    else
        m = permutedims(Matrix(df_notsrc))
        df_tmp = rename!(DataFrame(Tables.table(m)), new_col_names, makeunique=makeunique)
    end
    return hcat!(df_permuted, df_tmp, makeunique=makeunique, copycols=false)
end

function Base.permutedims(df::AbstractDataFrame, src_namescol::ColumnIndex;
                          makeunique::Bool=false, strict::Bool=true)
    if src_namescol isa Integer
        1 <= src_namescol <= ncol(df) || throw(BoundsError(index(df), src_namescol))
        dest_namescol = _names(df)[src_namescol]
    else
        dest_namescol = src_namescol
    end
    return permutedims(df, src_namescol, dest_namescol;
                       makeunique=makeunique, strict=strict)
end


text:
using Printf
using LinearAlgebra

include("../common/PWGrid_v01.jl")
include("../common/gen_lattice.jl")
include("../common/print_matrix.jl")

function write_XSF( filnam, LL_in, atpos; molecule=false )
    #
    f = open(filnam, "w")
    Natoms = size(atpos)[2]
    #
    if molecule
        @printf(f, "MOLECULE\n")
    else
        @printf(f, "CRYSTAL\n")
    end
    LL = LL_in'  # Xcrysden uses lattice vectors by rows
    @printf(f, "PRIMVEC\n")
    @printf(f, "%18.10f %18.10f %18.10f\n", LL[1,1], LL[1,2], LL[1,3])
    @printf(f, "%18.10f %18.10f %18.10f\n", LL[2,1], LL[2,2], LL[2,3])
    @printf(f, "%18.10f %18.10f %18.10f\n", LL[3,1], LL[3,2], LL[3,3])
    @printf(f, "PRIMCOORD\n")
    @printf(f, "%8d %8d\n", Natoms, 1)
    for ia = 1:Natoms
        @printf(f, "X %18.10f %18.10f %18.10f\n", atpos[1,ia], atpos[2,ia], atpos[3,ia])
    end
    close(f)
end

function test_main_hexagonal()

    #LL = 16.0*diagm([1.0, 1.0, 1.0])
    #LL = gen_lattice_fcc(16.0)
    #LL = gen_lattice_bcc(16.0)

    Ns = [10, 10, 20]
    LL = gen_lattice_hexagonal(10.0, coa=2.0)
    pw = PWGrid( Ns, LL )
    atpos = pw.r

    println("Lattice vectors (by column)")
    print_matrix(pw.LatVecs)
    println("Reciprocal lattice vectors (by column)")
    print_matrix(pw.RecVecs)

    write_XSF("R_grid_hexagonal.xsf", LL, atpos)

    Rec = pw.RecVecs*Ns[1]/2.0
    atpos = pw.G
    write_XSF("G_grid_hexagonal.xsf", LL, atpos, molecule=true)

    for ii = 1:3
        @printf("LatVecLen %d %18.10f\n", ii, norm(pw.LatVecs[:,ii]))
    end
    @printf("Ratio coa: %18.10f\n", norm(pw.LatVecs[:,1])/norm(pw.LatVecs[:,3]))

    @printf("\n")
    for ii = 1:3
        @printf("RecVecLen %d %18.10f\n", ii, norm(pw.RecVecs[ii,:]))
    end
    @printf("Ratio coa: %18.10f\n", norm(pw.RecVecs[:,3])/norm(pw.RecVecs[:,3]))

end

test_main_hexagonal()

ext: jl


text:
# generated
# length(dict) = 829
const repl_emoji_symbols_below_170 = Dict{String, String}(
    "\\:+1:" => "👍",
    "\\:-1:" => "👎",
    "\\:100:" => "💯",
    "\\:1234:" => "🔢",
    "\\:8ball:" => "🎱",
    "\\:a:" => "🅰",
    "\\:ab:" => "🆎",
    "\\:abc:" => "🔤",
    "\\:abcd:" => "🔡",
    "\\:accept:" => "🉑",
    "\\:aerial_tramway:" => "🚡",
    "\\:airplane:" => "✈",
    "\\:alarm_clock:" => "⏰",
    "\\:alien:" => "👽",
    "\\:ambulance:" => "🚑",
    "\\:anchor:" => "⚓",
    "\\:angel:" => "👼",
    "\\:anger:" => "💢",
    "\\:angry:" => "😠",
    "\\:anguished:" => "😧",
    "\\:ant:" => "🐜",
    "\\:apple:" => "🍎",
    "\\:aquarius:" => "♒",
    "\\:aries:" => "♈",
    "\\:arrow_backward:" => "◀",
    "\\:arrow_double_down:" => "⏬",
    "\\:arrow_double_up:" => "⏫",
    "\\:arrow_down:" => "⬇",
    "\\:arrow_down_small:" => "🔽",
    "\\:arrow_forward:" => "▶",
    "\\:arrow_heading_down:" => "⤵",
    "\\:arrow_heading_up:" => "⤴",
    "\\:arrow_left:" => "⬅",
    "\\:arrow_lower_left:" => "↙",
    "\\:arrow_lower_right:" => "↘",
    "\\:arrow_right:" => "➡",
    "\\:arrow_right_hook:" => "↪",
    "\\:arrow_up:" => "⬆",
    "\\:arrow_up_down:" => "↕",
    "\\:arrow_up_small:" => "🔼",
    "\\:arrow_upper_left:" => "↖",
    "\\:arrow_upper_right:" => "↗",
    "\\:arrows_clockwise:" => "🔃",
    "\\:arrows_counterclockwise:" => "🔄",
    "\\:art:" => "🎨",
    "\\:articulated_lorry:" => "🚛",
    "\\:astonished:" => "😲",
    "\\:athletic_shoe:" => "👟",
    "\\:atm:" => "🏧",
    "\\:b:" => "🅱",
    "\\:baby:" => "👶",
    "\\:baby_bottle:" => "🍼",
    "\\:baby_chick:" => "🐤",
    "\\:baby_symbol:" => "🚼",
    "\\:back:" => "🔙",
    "\\:baggage_claim:" => "🛄",
    "\\:balloon:" => "🎈",
    "\\:ballot_box_with_check:" => "☑",
    "\\:bamboo:" => "🎍",
    "\\:banana:" => "🍌",
    "\\:bangbang:" => "‼",
    "\\:bank:" => "🏦",
    "\\:bar_chart:" => "📊",
    "\\:barber:" => "💈",
    "\\:baseball:" => "⚾",
    "\\:basketball:" => "🏀",
    "\\:bath:" => "🛀",
    "\\:bathtub:" => "🛁",
    "\\:battery:" => "🔋",
    "\\:bear:" => "🐻",
    "\\:bee:" => "🐝",
    "\\:beer:" => "🍺",
    "\\:beers:" => "🍻",
    "\\:beetle:" => "🐞",
    "\\:beginner:" => "🔰",
    "\\:bell:" => "🔔",
    "\\:bento:" => "🍱",
    "\\:bicyclist:" => "🚴",
    "\\:bike:" => "🚲",
    "\\:bikini:" => "👙",
    "\\:bird:" => "🐦",
    "\\:birthday:" => "🎂",
    "\\:black_circle:" => "⚫",
    "\\:black_joker:" => "🃏",
    "\\:black_large_square:" => "⬛",
    "\\:black_medium_small_square:" => "◾",
    "\\:black_medium_square:" => "◼",
    "\\:black_nib:" => "✒",
    "\\:black_small_square:" => "▪",
    "\\:black_square_button:" => "🔲",
    "\\:blossom:" => "🌼",
    "\\:blowfish:" => "🐡",
    "\\:blue_book:" => "📘",
    "\\:blue_car:" => "🚙",
    "\\:blue_heart:" => "💙",
    "\\:blush:" => "😊",
    "\\:boar:" => "🐗",
    "\\:boat:" => "⛵",
    "\\:bomb:" => "💣",
    "\\:book:" => "📖",
    "\\:bookmark:" => "🔖",
    "\\:bookmark_tabs:" => "📑",
    "\\:books:" => "📚",
    "\\:boom:" => "💥",
    "\\:boot:" => "👢",
    "\\:bouquet:" => "💐",
    "\\:bow:" => "🙇",
    "\\:bowling:" => "🎳",
    "\\:boy:" => "👦",
    "\\:bread:" => "🍞",
    "\\:bride_with_veil:" => "👰",
    "\\:bridge_at_night:" => "🌉",
    "\\:briefcase:" => "💼",
    "\\:broken_heart:" => "💔",
    "\\:bug:" => "🐛",
    "\\:bulb:" => "💡",
    "\\:bullettrain_front:" => "🚅",
    "\\:bullettrain_side:" => "🚄",
    "\\:bus:" => "🚌",
    "\\:busstop:" => "🚏",
    "\\:bust_in_silhouette:" => "👤",
    "\\:busts_in_silhouette:" => "👥",
    "\\:cactus:" => "🌵",
    "\\:cake:" => "🍰",
    "\\:calendar:" => "📆",
    "\\:calling:" => "📲",
    "\\:camel:" => "🐫",
    "\\:camera:" => "📷",
    "\\:cancer:" => "♋",
    "\\:candy:" => "🍬",
    "\\:capital_abcd:" => "🔠",
    "\\:capricorn:" => "♑",
    "\\:car:" => "🚗",
    "\\:card_index:" => "📇",
    "\\:carousel_horse:" => "🎠",
    "\\:cat2:" => "🐈",
    "\\:cat:" => "🐱",
    "\\:cd:" => "💿",
    "\\:chart:" => "💹",
    "\\:chart_with_downwards_trend:" => "📉",
    "\\:chart_with_upwards_trend:" => "📈",
    "\\:checkered_flag:" => "🏁",
    "\\:cherries:" => "🍒",
    "\\:cherry_blossom:" => "🌸",
    "\\:chestnut:" => "🌰",
    "\\:chicken:" => "🐔",
    "\\:children_crossing:" => "🚸",
    "\\:chocolate_bar:" => "🍫",
    "\\:christmas_tree:" => "🎄",
    "\\:church:" => "⛪",
    "\\:cinema:" => "🎦",
    "\\:circus_tent:" => "🎪",
    "\\:city_sunrise:" => "🌇",
    "\\:city_sunset:" => "🌆",
    "\\:cl:" => "🆑",
    "\\:clap:" => "👏",
    "\\:clapper:" => "🎬",
    "\\:clipboard:" => "📋",
    "\\:clock1030:" => "🕥",
    "\\:clock10:" => "🕙",
    "\\:clock1130:" => "🕦",
    "\\:clock11:" => "🕚",
    "\\:clock1230:" => "🕧",
    "\\:clock12:" => "🕛",
    "\\:clock130:" => "🕜",
    "\\:clock1:" => "🕐",
    "\\:clock230:" => "🕝",
    "\\:clock2:" => "🕑",
    "\\:clock330:" => "🕞",
    "\\:clock3:" => "🕒",
    "\\:clock430:" => "🕟",
    "\\:clock4:" => "🕓",
    "\\:clock530:" => "🕠",
    "\\:clock5:" => "🕔",
    "\\:clock630:" => "🕡",
    "\\:clock6:" => "🕕",
    "\\:clock730:" => "🕢",
    "\\:clock7:" => "🕖",
    "\\:clock830:" => "🕣",
    "\\:clock8:" => "🕗",
    "\\:clock930:" => "🕤",
    "\\:clock9:" => "🕘",
    "\\:closed_book:" => "📕",
    "\\:closed_lock_with_key:" => "🔐",
    "\\:closed_umbrella:" => "🌂",
    "\\:cloud:" => "☁",
    "\\:clubs:" => "♣",
    "\\:cocktail:" => "🍸",
    "\\:coffee:" => "☕",
    "\\:cold_sweat:" => "😰",
    "\\:computer:" => "💻",
    "\\:confetti_ball:" => "🎊",
    "\\:confounded:" => "😖",
    "\\:confused:" => "😕",
    "\\:congratulations:" => "㊗",
    "\\:construction:" => "🚧",
    "\\:construction_worker:" => "👷",
    "\\:convenience_store:" => "🏪",
    "\\:cookie:" => "🍪",
    "\\:cool:" => "🆒",
    "\\:cop:" => "👮",
    "\\:copyright:" => "©",
    "\\:corn:" => "🌽",
    "\\:couple:" => "👫",
    "\\:couple_with_heart:" => "💑",
    "\\:couplekiss:" => "💏",
    "\\:cow2:" => "🐄",
    "\\:cow:" => "🐮",
    "\\:credit_card:" => "💳",
    "\\:crescent_moon:" => "🌙",
    "\\:crocodile:" => "🐊",
    "\\:crossed_flags:" => "🎌",
    "\\:crown:" => "👑",
    "\\:cry:" => "😢",
    "\\:crying_cat_face:" => "😿",
    "\\:crystal_ball:" => "🔮",
    "\\:cupid:" => "💘",
    "\\:curly_loop:" => "➰",
    "\\:currency_exchange:" => "💱",
    "\\:curry:" => "🍛",
    "\\:custard:" => "🍮",
    "\\:customs:" => "🛃",
    "\\:cyclone:" => "🌀",
    "\\:dancer:" => "💃",
    "\\:dancers:" => "👯",
    "\\:dango:" => "🍡",
    "\\:dart:" => "🎯",
    "\\:dash:" => "💨",
    "\\:date:" => "📅",
    "\\:deciduous_tree:" => "🌳",
    "\\:department_store:" => "🏬",
    "\\:diamond_shape_with_a_dot_inside:" => "💠",
    "\\:diamonds:" => "♦",
    "\\:disappointed:" => "😞",
    "\\:disappointed_relieved:" => "😥",
    "\\:dizzy:" => "💫",
    "\\:dizzy_face:" => "😵",
    "\\:do_not_litter:" => "🚯",
    "\\:dog2:" => "🐕",
    "\\:dog:" => "🐶",
    "\\:dollar:" => "💵",
    "\\:dolls:" => "🎎",
    "\\:dolphin:" => "🐬",
    "\\:door:" => "🚪",
    "\\:doughnut:" => "🍩",
    "\\:dragon:" => "🐉",
    "\\:dragon_face:" => "🐲",
    "\\:dress:" => "👗",
    "\\:dromedary_camel:" => "🐪",
    "\\:droplet:" => "💧",
    "\\:dvd:" => "📀",
    "\\:e-mail:" => "📧",
    "\\:ear:" => "👂",
    "\\:ear_of_rice:" => "🌾",
    "\\:earth_africa:" => "🌍",
    "\\:earth_americas:" => "🌎",
    "\\:earth_asia:" => "🌏",
    "\\:egg:" => "🍳",
    "\\:eggplant:" => "🍆",
    "\\:eight_pointed_black_star:" => "✴",
    "\\:eight_spoked_asterisk:" => "✳",
    "\\:electric_plug:" => "🔌",
    "\\:elephant:" => "🐘",
    "\\:email:" => "✉",
    "\\:end:" => "🔚",
    "\\:envelope_with_arrow:" => "📩",
    "\\:euro:" => "💶",
    "\\:european_castle:" => "🏰",
    "\\:european_post_office:" => "🏤",
    "\\:evergreen_tree:" => "🌲",
    "\\:exclamation:" => "❗",
    "\\:expressionless:" => "😑",
    "\\:eyeglasses:" => "👓",
    "\\:eyes:" => "👀",
    "\\:facepunch:" => "👊",
    "\\:factory:" => "🏭",
    "\\:fallen_leaf:" => "🍂",
    "\\:family:" => "👪",
    "\\:fast_forward:" => "⏩",
    "\\:fax:" => "📠",
    "\\:fearful:" => "😨",
    "\\:feet:" => "🐾",
    "\\:ferris_wheel:" => "🎡",
    "\\:file_folder:" => "📁",
    "\\:fire:" => "🔥",
    "\\:fire_engine:" => "🚒",
    "\\:fireworks:" => "🎆",
    "\\:first_quarter_moon:" => "🌓",
    "\\:first_quarter_moon_with_face:" => "🌛",
    "\\:fish:" => "🐟",
    "\\:fish_cake:" => "🍥",
    "\\:fishing_pole_and_fish:" => "🎣",
    "\\:fist:" => "✊",
    "\\:flags:" => "🎏",
    "\\:flashlight:" => "🔦",
    "\\:floppy_disk:" => "💾",
    "\\:flower_playing_cards:" => "🎴",
    "\\:flushed:" => "😳",
    "\\:foggy:" => "🌁",
    "\\:football:" => "🏈",
    "\\:footprints:" => "👣",
    "\\:fork_and_knife:" => "🍴",
    "\\:fountain:" => "⛲",
    "\\:four_leaf_clover:" => "🍀",
    "\\:free:" => "🆓",
    "\\:fried_shrimp:" => "🍤",
    "\\:fries:" => "🍟",
    "\\:frog:" => "🐸",
    "\\:frowning:" => "😦",
    "\\:fuelpump:" => "⛽",
    "\\:full_moon:" => "🌕",
    "\\:full_moon_with_face:" => "🌝",
    "\\:game_die:" => "🎲",
    "\\:gem:" => "💎",
    "\\:gemini:" => "♊",
    "\\:ghost:" => "👻",
    "\\:gift:" => "🎁",
    "\\:gift_heart:" => "💝",
    "\\:girl:" => "👧",
    "\\:globe_with_meridians:" => "🌐",
    "\\:goat:" => "🐐",
    "\\:golf:" => "⛳",
    "\\:grapes:" => "🍇",
    "\\:green_apple:" => "🍏",
    "\\:green_book:" => "📗",
    "\\:green_heart:" => "💚",
    "\\:grey_exclamation:" => "❕",
    "\\:grey_question:" => "❔",
    "\\:grimacing:" => "😬",
    "\\:grin:" => "😁",
    "\\:grinning:" => "😀",
    "\\:guardsman:" => "💂",
    "\\:guitar:" => "🎸",
    "\\:gun:" => "🔫",
    "\\:haircut:" => "💇",
    "\\:hamburger:" => "🍔",
    "\\:hammer:" => "🔨",
    "\\:hamster:" => "🐹",
    "\\:hand:" => "✋",
    "\\:handbag:" => "👜",
    "\\:hankey:" => "💩",
    "\\:hatched_chick:" => "🐥",
    "\\:hatching_chick:" => "🐣",
    "\\:headphones:" => "🎧",
    "\\:hear_no_evil:" => "🙉",
    "\\:heart:" => "❤",
    "\\:heart_decoration:" => "💟",
    "\\:heart_eyes:" => "😍",
    "\\:heart_eyes_cat:" => "😻",
    "\\:heartbeat:" => "💓",
    "\\:heartpulse:" => "💗",
    "\\:hearts:" => "♥",
    "\\:heavy_check_mark:" => "✔",
    "\\:heavy_division_sign:" => "➗",
    "\\:heavy_dollar_sign:" => "💲",
    "\\:heavy_minus_sign:" => "➖",
    "\\:heavy_multiplication_x:" => "✖",
    "\\:heavy_plus_sign:" => "➕",
    "\\:helicopter:" => "🚁",
    "\\:herb:" => "🌿",
    "\\:hibiscus:" => "🌺",
    "\\:high_brightness:" => "🔆",
    "\\:high_heel:" => "👠",
    "\\:hocho:" => "🔪",
    "\\:honey_pot:" => "🍯",
    "\\:horse:" => "🐴",
    "\\:horse_racing:" => "🏇",
    "\\:hospital:" => "🏥",
    "\\:hotel:" => "🏨",
    "\\:hotsprings:" => "♨",
    "\\:hourglass:" => "⌛",
    "\\:hourglass_flowing_sand:" => "⏳",
    "\\:house:" => "🏠",
    "\\:house_with_garden:" => "🏡",
    "\\:hushed:" => "😯",
    "\\:ice_cream:" => "🍨",
    "\\:icecream:" => "🍦",
    "\\:id:" => "🆔",
    "\\:ideograph_advantage:" => "🉐",
    "\\:imp:" => "👿",
    "\\:inbox_tray:" => "📥",
    "\\:incoming_envelope:" => "📨",
    "\\:information_desk_person:" => "💁",
    "\\:information_source:" => "ℹ",
    "\\:innocent:" => "😇",
    "\\:interrobang:" => "⁉",
    "\\:iphone:" => "📱",
    "\\:izakaya_lantern:" => "🏮",
    "\\:jack_o_lantern:" => "🎃",
    "\\:japan:" => "🗾",
    "\\:japanese_castle:" => "🏯",
    "\\:japanese_goblin:" => "👺",
    "\\:japanese_ogre:" => "👹",
    "\\:jeans:" => "👖",
    "\\:joy:" => "😂",
    "\\:joy_cat:" => "😹",
    "\\:key:" => "🔑",
    "\\:keycap_ten:" => "🔟",
    "\\:kimono:" => "👘",
    "\\:kiss:" => "💋",
    "\\:kissing:" => "😗",
    "\\:kissing_cat:" => "😽",
    "\\:kissing_closed_eyes:" => "😚",
    "\\:kissing_heart:" => "😘",
    "\\:kissing_smiling_eyes:" => "😙",
    "\\:koala:" => "🐨",
    "\\:koko:" => "🈁",
    "\\:large_blue_circle:" => "🔵",
    "\\:large_blue_diamond:" => "🔷",
    "\\:large_orange_diamond:" => "🔶",
    "\\:last_quarter_moon:" => "🌗",
    "\\:last_quarter_moon_with_face:" => "🌜",
    "\\:laughing:" => "😆",
    "\\:leaves:" => "🍃",
    "\\:ledger:" => "📒",
    "\\:left_luggage:" => "🛅",
    "\\:left_right_arrow:" => "↔",
    "\\:leftwards_arrow_with_hook:" => "↩",
    "\\:lemon:" => "🍋",
    "\\:leo:" => "♌",
    "\\:leopard:" => "🐆",
    "\\:libra:" => "♎",
    "\\:light_rail:" => "🚈",
    "\\:link:" => "🔗",
    "\\:lips:" => "👄",
    "\\:lipstick:" => "💄",
    "\\:lock:" => "🔒",
    "\\:lock_with_ink_pen:" => "🔏",
    "\\:lollipop:" => "🍭",
    "\\:loop:" => "➿",
    "\\:loud_sound:" => "🔊",
    "\\:loudspeaker:" => "📢",
    "\\:love_hotel:" => "🏩",
    "\\:love_letter:" => "💌",
    "\\:low_brightness:" => "🔅",
    "\\:m:" => "Ⓜ",
    "\\:mag:" => "🔍",
    "\\:mag_right:" => "🔎",
    "\\:mahjong:" => "🀄",
    "\\:mailbox:" => "📫",
    "\\:mailbox_closed:" => "📪",
    "\\:mailbox_with_mail:" => "📬",
    "\\:mailbox_with_no_mail:" => "📭",
    "\\:man:" => "👨",
    "\\:man_with_gua_pi_mao:" => "👲",
    "\\:man_with_turban:" => "👳",
    "\\:mans_shoe:" => "👞",
    "\\:maple_leaf:" => "🍁",
    "\\:mask:" => "😷",
    "\\:massage:" => "💆",
    "\\:meat_on_bone:" => "🍖",
    "\\:mega:" => "📣",
    "\\:melon:" => "🍈",
    "\\:memo:" => "📝",
    "\\:mens:" => "🚹",
    "\\:metro:" => "🚇",
    "\\:microphone:" => "🎤",
    "\\:microscope:" => "🔬",
    "\\:milky_way:" => "🌌",
    "\\:minibus:" => "🚐",
    "\\:minidisc:" => "💽",
    "\\:mobile_phone_off:" => "📴",
    "\\:money_with_wings:" => "💸",
    "\\:moneybag:" => "💰",
    "\\:monkey:" => "🐒",
    "\\:monkey_face:" => "🐵",
    "\\:monorail:" => "🚝",
    "\\:moon:" => "🌔",
    "\\:mortar_board:" => "🎓",
    "\\:mount_fuji:" => "🗻",
    "\\:mountain_bicyclist:" => "🚵",
    "\\:mountain_cableway:" => "🚠",
    "\\:mountain_railway:" => "🚞",
    "\\:mouse2:" => "🐁",
    "\\:mouse:" => "🐭",
    "\\:movie_camera:" => "🎥",
    "\\:moyai:" => "🗿",
    "\\:muscle:" => "💪",
    "\\:mushroom:" => "🍄",
    "\\:musical_keyboard:" => "🎹",
    "\\:musical_note:" => "🎵",
    "\\:musical_score:" => "🎼",
    "\\:mute:" => "🔇",
    "\\:nail_care:" => "💅",
    "\\:name_badge:" => "📛",
    "\\:necktie:" => "👔",
    "\\:negative_squared_cross_mark:" => "❎",
    "\\:neutral_face:" => "😐",
    "\\:new:" => "🆕",
    "\\:new_moon:" => "🌑",
    "\\:new_moon_with_face:" => "🌚",
    "\\:newspaper:" => "📰",
    "\\:ng:" => "🆖",
    "\\:night_with_stars:" => "🌃",
    "\\:no_bell:" => "🔕",
    "\\:no_bicycles:" => "🚳",
    "\\:no_entry:" => "⛔",
    "\\:no_entry_sign:" => "🚫",
    "\\:no_good:" => "🙅",
    "\\:no_mobile_phones:" => "📵",
    "\\:no_mouth:" => "😶",
    "\\:no_pedestrians:" => "🚷",
    "\\:no_smoking:" => "🚭",
    "\\:non-potable_water:" => "🚱",
    "\\:nose:" => "👃",
    "\\:notebook:" => "📓",
    "\\:notebook_with_decorative_cover:" => "📔",
    "\\:notes:" => "🎶",
    "\\:nut_and_bolt:" => "🔩",
    "\\:o2:" => "🅾",
    "\\:o:" => "⭕",
    "\\:ocean:" => "🌊",
    "\\:octopus:" => "🐙",
    "\\:oden:" => "🍢",
    "\\:office:" => "🏢",
    "\\:ok:" => "🆗",
    "\\:ok_hand:" => "👌",
    "\\:ok_woman:" => "🙆",
    "\\:older_man:" => "👴",
    "\\:older_woman:" => "👵",
    "\\:on:" => "🔛",
    "\\:oncoming_automobile:" => "🚘",
    "\\:oncoming_bus:" => "🚍",
    "\\:oncoming_police_car:" => "🚔",
    "\\:oncoming_taxi:" => "🚖",
    "\\:open_file_folder:" => "📂",
    "\\:open_hands:" => "👐",
    "\\:open_mouth:" => "😮",
    "\\:ophiuchus:" => "⛎",
    "\\:orange_book:" => "📙",
    "\\:outbox_tray:" => "📤",
    "\\:ox:" => "🐂",
    "\\:package:" => "📦",
    "\\:page_facing_up:" => "📄",
    "\\:page_with_curl:" => "📃",
    "\\:pager:" => "📟",
    "\\:palm_tree:" => "🌴",
    "\\:panda_face:" => "🐼",
    "\\:paperclip:" => "📎",
    "\\:parking:" => "🅿",
    "\\:part_alternation_mark:" => "〽",
    "\\:partly_sunny:" => "⛅",
    "\\:passport_control:" => "🛂",
    "\\:peach:" => "🍑",
    "\\:pear:" => "🍐",
    "\\:pencil2:" => "✏",
    "\\:penguin:" => "🐧",
    "\\:pensive:" => "😔",
    "\\:performing_arts:" => "🎭",
    "\\:persevere:" => "😣",
    "\\:person_frowning:" => "🙍",
    "\\:person_with_blond_hair:" => "👱",
    "\\:person_with_pouting_face:" => "🙎",
    "\\:phone:" => "☎",
    "\\:pig2:" => "🐖",
    "\\:pig:" => "🐷",
    "\\:pig_nose:" => "🐽",
    "\\:pill:" => "💊",
    "\\:pineapple:" => "🍍",
    "\\:pisces:" => "♓",
    "\\:pizza:" => "🍕",
    "\\:point_down:" => "👇",
    "\\:point_left:" => "👈",
    "\\:point_right:" => "👉",
    "\\:point_up:" => "☝",
    "\\:point_up_2:" => "👆",
    "\\:police_car:" => "🚓",
    "\\:poodle:" => "🐩",
    "\\:post_office:" => "🏣",
    "\\:postal_horn:" => "📯",
    "\\:postbox:" => "📮",
    "\\:potable_water:" => "🚰",
    "\\:pouch:" => "👝",
    "\\:poultry_leg:" => "🍗",
    "\\:pound:" => "💷",
    "\\:pouting_cat:" => "😾",
    "\\:pray:" => "🙏",
    "\\:princess:" => "👸",
    "\\:purple_heart:" => "💜",
    "\\:purse:" => "👛",
    "\\:pushpin:" => "📌",
    "\\:put_litter_in_its_place:" => "🚮",
    "\\:question:" => "❓",
    "\\:rabbit2:" => "🐇",
    "\\:rabbit:" => "🐰",
    "\\:racehorse:" => "🐎",
    "\\:radio:" => "📻",
    "\\:radio_button:" => "🔘",
    "\\:rage:" => "😡",
    "\\:railway_car:" => "🚃",
    "\\:rainbow:" => "🌈",
    "\\:raised_hands:" => "🙌",
    "\\:raising_hand:" => "🙋",
    "\\:ram:" => "🐏",
    "\\:ramen:" => "🍜",
    "\\:rat:" => "🐀",
    "\\:recycle:" => "♻",
    "\\:red_circle:" => "🔴",
    "\\:registered:" => "®",
    "\\:relaxed:" => "☺",
    "\\:relieved:" => "😌",
    "\\:repeat:" => "🔁",
    "\\:repeat_one:" => "🔂",
    "\\:restroom:" => "🚻",
    "\\:revolving_hearts:" => "💞",
    "\\:rewind:" => "⏪",
    "\\:ribbon:" => "🎀",
    "\\:rice:" => "🍚",
    "\\:rice_ball:" => "🍙",
    "\\:rice_cracker:" => "🍘",
    "\\:rice_scene:" => "🎑",
    "\\:ring:" => "💍",
    "\\:rocket:" => "🚀",
    "\\:roller_coaster:" => "🎢",
    "\\:rooster:" => "🐓",
    "\\:rose:" => "🌹",
    "\\:rotating_light:" => "🚨",
    "\\:round_pushpin:" => "📍",
    "\\:rowboat:" => "🚣",
    "\\:rugby_football:" => "🏉",
    "\\:runner:" => "🏃",
    "\\:running_shirt_with_sash:" => "🎽",
    "\\:sa:" => "🈂",
    "\\:sagittarius:" => "♐",
    "\\:sake:" => "🍶",
    "\\:sandal:" => "👡",
    "\\:santa:" => "🎅",
    "\\:satellite:" => "📡",
    "\\:saxophone:" => "🎷",
    "\\:school:" => "🏫",
    "\\:school_satchel:" => "🎒",
    "\\:scissors:" => "✂",
    "\\:scorpius:" => "♏",
    "\\:scream:" => "😱",
    "\\:scream_cat:" => "🙀",
    "\\:scroll:" => "📜",
    "\\:seat:" => "💺",
    "\\:secret:" => "㊙",
    "\\:see_no_evil:" => "🙈",
    "\\:seedling:" => "🌱",
    "\\:shaved_ice:" => "🍧",
    "\\:sheep:" => "🐑",
    "\\:shell:" => "🐚",
    "\\:ship:" => "🚢",
    "\\:shirt:" => "👕",
    "\\:shower:" => "🚿",
    "\\:signal_strength:" => "📶",
    "\\:six_pointed_star:" => "🔯",
    "\\:ski:" => "🎿",
    "\\:skin-tone-2:" => "🏻",
    "\\:skin-tone-3:" => "🏼",
    "\\:skin-tone-4:" => "🏽",
    "\\:skin-tone-5:" => "🏾",
    "\\:skin-tone-6:" => "🏿",
    "\\:skull:" => "💀",
    "\\:sleeping:" => "😴",
    "\\:sleepy:" => "😪",
    "\\:slot_machine:" => "🎰",
    "\\:small_blue_diamond:" => "🔹",
    "\\:small_orange_diamond:" => "🔸",
    "\\:small_red_triangle:" => "🔺",
    "\\:small_red_triangle_down:" => "🔻",
    "\\:smile:" => "😄",
    "\\:smile_cat:" => "😸",
    "\\:smiley:" => "😃",
    "\\:smiley_cat:" => "😺",
    "\\:smiling_imp:" => "😈",
    "\\:smirk:" => "😏",
    "\\:smirk_cat:" => "😼",
    "\\:smoking:" => "🚬",
    "\\:snail:" => "🐌",
    "\\:snake:" => "🐍",
    "\\:snowboarder:" => "🏂",
    "\\:snowflake:" => "❄",
    "\\:snowman:" => "⛄",
    "\\:sob:" => "😭",
    "\\:soccer:" => "⚽",
    "\\:soon:" => "🔜",
    "\\:sos:" => "🆘",
    "\\:sound:" => "🔉",
    "\\:space_invader:" => "👾",
    "\\:spades:" => "♠",
    "\\:spaghetti:" => "🍝",
    "\\:sparkle:" => "❇",
    "\\:sparkler:" => "🎇",
    "\\:sparkles:" => "✨",
    "\\:sparkling_heart:" => "💖",
    "\\:speak_no_evil:" => "🙊",
    "\\:speaker:" => "🔈",
    "\\:speech_balloon:" => "💬",
    "\\:speedboat:" => "🚤",
    "\\:star2:" => "🌟",
    "\\:star:" => "⭐",
    "\\:stars:" => "🌠",
    "\\:station:" => "🚉",
    "\\:statue_of_liberty:" => "🗽",
    "\\:steam_locomotive:" => "🚂",
    "\\:stew:" => "🍲",
    "\\:straight_ruler:" => "📏",
    "\\:strawberry:" => "🍓",
    "\\:stuck_out_tongue:" => "😛",
    "\\:stuck_out_tongue_closed_eyes:" => "😝",
    "\\:stuck_out_tongue_winking_eye:" => "😜",
    "\\:sun_with_face:" => "🌞",
    "\\:sunflower:" => "🌻",
    "\\:sunglasses:" => "😎",
    "\\:sunny:" => "☀",
    "\\:sunrise:" => "🌅",
    "\\:sunrise_over_mountains:" => "🌄",
    "\\:surfer:" => "🏄",
    "\\:sushi:" => "🍣",
    "\\:suspension_railway:" => "🚟",
    "\\:sweat:" => "😓",
    "\\:sweat_drops:" => "💦",
    "\\:sweat_smile:" => "😅",
    "\\:sweet_potato:" => "🍠",
    "\\:swimmer:" => "🏊",
    "\\:symbols:" => "🔣",
    "\\:syringe:" => "💉",
    "\\:tada:" => "🎉",
    "\\:tanabata_tree:" => "🎋",
    "\\:tangerine:" => "🍊",
    "\\:taurus:" => "♉",
    "\\:taxi:" => "🚕",
    "\\:tea:" => "🍵",
    "\\:telephone_receiver:" => "📞",
    "\\:telescope:" => "🔭",
    "\\:tennis:" => "🎾",
    "\\:tent:" => "⛺",
    "\\:thought_balloon:" => "💭",
    "\\:ticket:" => "🎫",
    "\\:tiger2:" => "🐅",
    "\\:tiger:" => "🐯",
    "\\:tired_face:" => "😫",
    "\\:tm:" => "™",
    "\\:toilet:" => "🚽",
    "\\:tokyo_tower:" => "🗼",
    "\\:tomato:" => "🍅",
    "\\:tongue:" => "👅",
    "\\:top:" => "🔝",
    "\\:tophat:" => "🎩",
    "\\:tractor:" => "🚜",
    "\\:traffic_light:" => "🚥",
    "\\:train2:" => "🚆",
    "\\:train:" => "🚋",
    "\\:tram:" => "🚊",
    "\\:triangular_flag_on_post:" => "🚩",
    "\\:triangular_ruler:" => "📐",
    "\\:trident:" => "🔱",
    "\\:triumph:" => "😤",
    "\\:trolleybus:" => "🚎",
    "\\:trophy:" => "🏆",
    "\\:tropical_drink:" => "🍹",
    "\\:tropical_fish:" => "🐠",
    "\\:truck:" => "🚚",
    "\\:trumpet:" => "🎺",
    "\\:tulip:" => "🌷",
    "\\:turtle:" => "🐢",
    "\\:tv:" => "📺",
    "\\:twisted_rightwards_arrows:" => "🔀",
    "\\:two_hearts:" => "💕",
    "\\:two_men_holding_hands:" => "👬",
    "\\:two_women_holding_hands:" => "👭",
    "\\:u5272:" => "🈹",
    "\\:u5408:" => "🈴",
    "\\:u55b6:" => "🈺",
    "\\:u6307:" => "🈯",
    "\\:u6708:" => "🈷",
    "\\:u6709:" => "🈶",
    "\\:u6e80:" => "🈵",
    "\\:u7121:" => "🈚",
    "\\:u7533:" => "🈸",
    "\\:u7981:" => "🈲",
    "\\:u7a7a:" => "🈳",
    "\\:umbrella:" => "☔",
    "\\:unamused:" => "😒",
    "\\:underage:" => "🔞",
    "\\:unlock:" => "🔓",
    "\\:up:" => "🆙",
    "\\:v:" => "✌",
    "\\:vertical_traffic_light:" => "🚦",
    "\\:vhs:" => "📼",
    "\\:vibration_mode:" => "📳",
    "\\:video_camera:" => "📹",
    "\\:video_game:" => "🎮",
    "\\:violin:" => "🎻",
    "\\:virgo:" => "♍",
    "\\:volcano:" => "🌋",
    "\\:vs:" => "🆚",
    "\\:walking:" => "🚶",
    "\\:waning_crescent_moon:" => "🌘",
    "\\:waning_gibbous_moon:" => "🌖",
    "\\:warning:" => "⚠",
    "\\:watch:" => "⌚",
    "\\:water_buffalo:" => "🐃",
    "\\:watermelon:" => "🍉",
    "\\:wave:" => "👋",
    "\\:wavy_dash:" => "〰",
    "\\:waxing_crescent_moon:" => "🌒",
    "\\:wc:" => "🚾",
    "\\:weary:" => "😩",
    "\\:wedding:" => "💒",
    "\\:whale2:" => "🐋",
    "\\:whale:" => "🐳",
    "\\:wheelchair:" => "♿",
    "\\:white_check_mark:" => "✅",
    "\\:white_circle:" => "⚪",
    "\\:white_flower:" => "💮",
    "\\:white_large_square:" => "⬜",
    "\\:white_medium_small_square:" => "◽",
    "\\:white_medium_square:" => "◻",
    "\\:white_small_square:" => "▫",
    "\\:white_square_button:" => "🔳",
    "\\:wind_chime:" => "🎐",
    "\\:wine_glass:" => "🍷",
    "\\:wink:" => "😉",
    "\\:wolf:" => "🐺",
    "\\:woman:" => "👩",
    "\\:womans_clothes:" => "👚",
    "\\:womans_hat:" => "👒",
    "\\:womens:" => "🚺",
    "\\:worried:" => "😟",
    "\\:wrench:" => "🔧",
    "\\:x:" => "❌",
    "\\:yellow_heart:" => "💛",
    "\\:yen:" => "💴",
    "\\:yum:" => "😋",
    "\\:zap:" => "⚡",
    "\\:zzz:" => "💤",
)


text:
# Because `deepcopy(mod::Module)` throws an error, we need a safe approach.
# Strategy: wrap the IdDict so that our methods get called rather than Base's.
# It's not guaranteed to work for user types that specialize `deepcopy_internal`,
# but hopefully that's rare.

struct WrappedIdDict
    dict::IdDict
end
Base.getindex(w::WrappedIdDict, key) = w.dict[key]
Base.setindex!(w::WrappedIdDict, val, key) = w.dict[key] = val
Base.haskey(w::WrappedIdDict, k) = haskey(w.dict, k)

function safe_deepcopy(a, args...)
    stackdict = WrappedIdDict(IdDict())
    _safe_deepcopy(stackdict, a, args...)
end
safe_deepcopy() = ()
_safe_deepcopy(stackdict, a, args...) =
    (Base.deepcopy_internal(a, stackdict), _safe_deepcopy(stackdict, args...)...)
_safe_deepcopy(stackdict) = ()

# This is the one method we want to override
Base.deepcopy_internal(x::Module, stackdict::WrappedIdDict) = x

# But the rest are necessary to make it work. This is just a direct copy from Base.
Base.deepcopy_internal(x::Union{Symbol,Core.MethodInstance,Method,GlobalRef,DataType,Union,Task},
                  stackdict::WrappedIdDict) = x
Base.deepcopy_internal(x::Tuple, stackdict::WrappedIdDict) =
    ntuple(i->Base.deepcopy_internal(x[i], stackdict), length(x))

function Base.deepcopy_internal(x::Base.SimpleVector, stackdict::WrappedIdDict)
    if haskey(stackdict, x)
        return stackdict[x]
    end
    y = Core.svec(Any[Base.deepcopy_internal(x[i], stackdict) for i = 1:length(x)]...)
    stackdict[x] = y
    return y
end

function Base.deepcopy_internal(x::String, stackdict::WrappedIdDict)
    if haskey(stackdict, x)
        return stackdict[x]
    end
    y = GC.@preserve x unsafe_string(pointer(x), sizeof(x))
    stackdict[x] = y
    return y
end

function Base.deepcopy_internal(@nospecialize(x), stackdict::WrappedIdDict)
    T = typeof(x)::DataType
    nf = nfields(x)
    (isbitstype(T) || nf == 0) && return x
    if haskey(stackdict, x)
        return stackdict[x]
    end
    if T.mutable
        y = ccall(:jl_new_struct_uninit, Any, (Any,), T)
        stackdict[x] = y
        for i in 1:nf
            if isdefined(x,i)
                ccall(:jl_set_nth_field, Cvoid, (Any, Csize_t, Any), y, i-1,
                      Base.deepcopy_internal(getfield(x,i), stackdict))
            end
        end
    else
        flds = Vector{Any}(undef, nf)
        for i in 1:nf
            if isdefined(x, i)
                xi = getfield(x, i)
                xi = Base.deepcopy_internal(xi, stackdict)::typeof(xi)
                flds[i] = xi
            else
                nf = i - 1 # rest of tail must be undefined values
                break
            end
        end
        y = ccall(:jl_new_structv, Any, (Any, Ptr{Any}, UInt32), T, flds, nf)
    end
    return y::T
end

function Base.deepcopy_internal(x::Array, stackdict::WrappedIdDict)
    if haskey(stackdict, x)
        return stackdict[x]
    end
    _deepcopy_array_t(x, eltype(x), stackdict)
end

function _deepcopy_array_t(@nospecialize(x), T, stackdict::WrappedIdDict)
    if isbitstype(T)
        return (stackdict[x]=copy(x))
    end
    dest = similar(x)
    stackdict[x] = dest
    for i = 1:(length(x)::Int)
        if ccall(:jl_array_isassigned, Cint, (Any, Csize_t), x, i-1) != 0
            xi = ccall(:jl_arrayref, Any, (Any, Csize_t), x, i-1)
            if !isbits(xi)
                xi = Base.deepcopy_internal(xi, stackdict)
            end
            ccall(:jl_arrayset, Cvoid, (Any, Any, Csize_t), dest, xi, i-1)
        end
    end
    return dest
end

function Base.deepcopy_internal(x::Union{Dict,IdDict}, stackdict::WrappedIdDict)
    if haskey(stackdict, x)
        return stackdict[x]::typeof(x)
    end

    if isbitstype(eltype(x))
        return (stackdict[x] = copy(x))
    end

    dest = empty(x)
    stackdict[x] = dest
    for (k, v) in x
        dest[Base.deepcopy_internal(k, stackdict)] = Base.deepcopy_internal(v, stackdict)
    end
    dest
end


text:
cdir = pwd(); ii = findlast("MDP", cdir); MDPpath = cdir[1:ii[end]] * "/";
include(MDPpath * "Installation/setpath.jl");

# Read xyz file to obtain configurations
filename = "/Users/ngoccuongnguyen/Dropbox (MIT)/MDP/data/LJCluster/curated_lj_cluster.xyz";
species = ["Ar"];
config = Preprocessing.readEXTXYZ(filename, species);

# shift atom positions to put them in the box [0,10]x[0,10]x[0,10]
config.x = config.x .+ 5.0

# cumalative sum of numbers of atoms
natom = [0; cumsum(config.natom[:])];

ci = 1; # configuration index in the data
x = config.x[:,(natom[ci]+1):natom[ci+1]] # atom positions of configuration ci
t = config.t[:,(natom[ci]+1):natom[ci+1]]                        # atom types of configuration ci
q = []

# 3 vectors define the simulation box
a = [10.0; 0.0; 0.0]
b = [0.0; 10.0; 0.0]
c = [0.0; 0.0; 10.0]

# periodic boundary conditions
pbc = [0; 0; 0]

# LJ parameters
epsilon = 1.0;
sigma = 1.0;
A = 2*epsilon*sigma^12;
B = 2*epsilon*sigma^6;
rcut = 4.0;

# define a list of potentials
potential = Array{Any}(undef, 1)
potential[1] = Potential.addpotential("pair", "nonbonded", [0], "ljpotential", [A; B], rcut);
for i = 1:length(potential)
    include(potential[i].potentialfunction * ".jl");  # include the potential file
end

# compute total energy, peratom energies, and forces
rcutmax = rcut; eta = [0]; kappa = [0];
etot, e, f = Potential.empiricalpotential(config.x[:,1:13], q, t, a, b, c, pbc, rcutmax, eta, kappa, potential);

# snap parameters
twojmax = Int32(6)
ntypes = Int32(1)
ncoeffall = Int32(31)
rcutfac = 4.0;
rfac0 = 1.0;
rmin0 = 0;
bzeroflag = 0;
quadraticflag = 0;
switchflag = 1;
chemflag = 0;
bnormflag = 0;
wselfallflag = 0;

snaparam = [ntypes, twojmax, rcutfac, rfac0, rmin0, bzeroflag, switchflag, quadraticflag, chemflag, bnormflag, wselfallflag]
elemradius = 0.5*ones(Float64, ntypes)
elemweight = ones(Float64, ntypes)

# intialize sna struct to store snap parameters
sna = Potential.initsna(snaparam, elemradius, elemweight)

# bispectrum, bispectrum derivatives, bispectrum virial
bi, bd, bv = Potential.snapdescriptors(x, t, a, b, c, pbc, sna)

# using LAMMPS.jl to compute energy and forces
unitstyle = "lj"
pair_style = "lj/cut $rcut"
pair_coeff = ["* * $epsilon $sigma"]
lbi, lbd, lbv, re, rf, rv = Potential.lammpssnapdescriptors(x, t, a, b, c, pbc,
            unitstyle, pair_style, pair_coeff, snaparam, elemradius, elemweight)

# compare with LAMMPS.jl
print("Maximum error in bispectrum components: "); display(maximum(abs.(bi[:] .- lbi[:])));
print("Maximum error in bispectrum derivatives: "); display(maximum(abs.(bd[:] .- lbd[:])));
print("Maximum error in bispectrum virials: "); display(maximum(abs.(bv[:] .- lbv[:])));
print("Maximum error in reference energy: "); display(abs(etot - re));
print("Maximum error in reference forces: "); display(maximum(abs.(f[:] .- rf[:])));

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
mutable struct Client
    id::Symbol
    stream::Any
    stream_accesstoken::Channel{Nothing}
    connected_notebook::Union{Notebook,Nothing}
    pendingupdates::Channel
end

Client(id::Symbol, stream) = let
    at = Channel{Nothing}(1)
    put!(at, nothing)
    Client(id, stream, at, nothing, Channel(1024))
end

struct Initiator
    clientID::Symbol
    requestID::Symbol
end

struct UpdateMessage
    type::Symbol
    message::Any
    notebook::Union{Notebook,Nothing}
    cell::Union{Cell,Nothing}
    initiator::Union{Initiator,Missing}
end

UpdateMessage(type::Symbol, message::Any) = UpdateMessage(type, message, nothing, nothing, missing)
UpdateMessage(type::Symbol, message::Any, notebook::Notebook) = UpdateMessage(type, message, notebook, nothing, missing)

function clientupdate_cell_output(notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)
    payload, mime = cell.output_repr, cell.repr_mime

    return UpdateMessage(:cell_output,
            Dict(:mime => mime,
             :output => payload,
             :errormessage => cell.error_repr,
             :runtime => cell.runtime,
            ),
            notebook, cell, initiator)
end

function clientupdate_cell_input(notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)
    return UpdateMessage(:cell_input,
        Dict(:code => cell.code), notebook, cell, initiator)
end

function clientupdate_cell_added(notebook::Notebook, cell::Cell, new_index::Integer; initiator::Union{Initiator,Missing}=missing)
    return UpdateMessage(:cell_added,
        Dict(:index => new_index - 1, # 1-based index (julia) to 0-based index (js)
            ), notebook, cell, initiator)
end

function clientupdate_cell_deleted(notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)
    return UpdateMessage(:cell_deleted,
        Dict(), notebook, cell, initiator)
end

function clientupdate_cell_moved(notebook::Notebook, cell::Cell, new_index::Integer; initiator::Union{Initiator,Missing}=missing)
    return UpdateMessage(:cell_moved,
        Dict(:index => new_index - 1, # 1-based index (julia) to 0-based index (js)
            ), notebook, cell, initiator)
end

function clientupdate_cell_running(notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)
    return UpdateMessage(:cell_running,
        Dict(), notebook, cell, initiator)
end

function clientupdate_notebook_list(notebooks; initiator::Union{Initiator,Missing}=missing)

    short_paths = Dict()

    notebookpaths = map(values(notebooks)) do notebook
        pathsep = Sys.iswindows() ? '\\' : '/'
        path_split = split(notebook.path, pathsep)
        if path_split[1] == ""
            path_split = path_split[2:end]
        end
        NotebookPath(notebook.uuid, path_split, "", -1)
    end

    make_paths_distinct!(Set(notebookpaths))

    short_paths = Dict(map(notebookpaths) do np
        np.uuid => np.current_path
    end...)

    update = UpdateMessage(:notebook_list,
        Dict(:notebooks => [Dict(
                :uuid => string(notebook.uuid),
                :path => notebook.path,
                :shortpath => short_paths[notebook.uuid]
                ) for notebook in values(notebooks)]), nothing, nothing, initiator)
end


text:
using UnitCellModelling
using Test

@testset "UnitCellModelling.jl" begin
    # Write your tests here.
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    @assert precompile(Tuple{typeof(maybe_evaluate_builtin), Frame, Expr, Bool})
    @assert precompile(Tuple{typeof(getargs), Vector{Any}, Frame})
    @assert precompile(Tuple{typeof(get_call_framecode), Vector{Any}, FrameCode, Int})
    @assert precompile(evaluate_call_recurse!, (Function, Frame, Expr))
    @assert precompile(evaluate_call_compiled!, (Compiled, Frame, Expr))
    for f in (evaluate_structtype,
              evaluate_abstracttype,
              evaluate_primitivetype)
        @assert precompile(Tuple{typeof(f), Any, Frame, Expr})
    end
    @assert precompile(Tuple{typeof(evaluate_foreigncall), Frame, Expr})
    @assert precompile(Tuple{typeof(evaluate_methoddef), Frame, Expr})
    @assert precompile(Tuple{typeof(lookup_global_refs!), Expr})
    @assert precompile(Tuple{typeof(lookup_or_eval), Any, Frame, Any})
    @assert precompile(Tuple{typeof(eval_rhs), Any, Frame, Expr})
    @assert precompile(Tuple{typeof(step_expr!), Any, Frame, Any, Bool})
    for f in (finish!, finish_and_return!, finish_stack!, next_call!, maybe_next_call!, next_line!)
        @assert precompile(Tuple{typeof(f), Any, Frame, Bool})
    end
    @assert precompile(Tuple{typeof(through_methoddef_or_done!), Any, Frame})
    @assert precompile(Tuple{typeof(split_expressions), Module, Expr})
    @assert precompile(Tuple{typeof(split_expressions!), Vector{Tuple{Module,Expr}}, Dict{Module,Vector{Expr}}, Expr, Module, Expr})
    @assert precompile(Tuple{typeof(prepare_thunk), Module, Expr})
    @assert precompile(Tuple{typeof(prepare_thunk), Module, Expr, Bool})
    @assert precompile(Tuple{typeof(prepare_framedata), FrameCode, Vector{Any}, SimpleVector, Bool})
    @assert precompile(Tuple{typeof(prepare_args), Any, Vector{Any}, Vector{Any}})
    @assert precompile(Tuple{typeof(prepare_call), Any, Vector{Any}})
    @assert precompile(Tuple{typeof(Core.kwfunc(prepare_call)), NamedTuple{(:enter_generated,),Tuple{Bool}}, typeof(prepare_call), Function, Vector{Any}})
    @assert precompile(Tuple{typeof(Core.kwfunc(prepare_framecode)), NamedTuple{(:enter_generated,),Tuple{Bool}}, typeof(prepare_framecode), Method, Tuple{Int}})
    @assert precompile(Tuple{typeof(prepare_frame), FrameCode, Vector{Any}, Core.SimpleVector})
    @assert precompile(Tuple{typeof(extract_args), Module, Expr})
    @assert precompile(Tuple{typeof(enter_call), Int, Int})
    @assert precompile(Tuple{typeof(enter_call_expr), Expr})
    @assert precompile(Tuple{typeof(copy_codeinfo), Core.CodeInfo})
    @assert precompile(Tuple{typeof(optimize!), Core.CodeInfo, Module})
    @assert precompile(Tuple{typeof(optimize!), Core.CodeInfo, Method})
    @assert precompile(Tuple{typeof(set_structtype_const), Module, Symbol})
    @assert precompile(Tuple{typeof(namedtuple), Vector{Any}})
    @assert precompile(Tuple{typeof(resolvefc), Frame, Any})
    @assert precompile(Tuple{typeof(check_isdefined), Frame, Any})
    @assert precompile(Tuple{typeof(find_used), Core.CodeInfo})
    @assert precompile(Tuple{typeof(do_assignment!), Frame, Any, Any})
    @assert precompile(Tuple{typeof(pc_expr), Frame})
end


text:
"""
`Quarter(v) -> Quarter`

Construct a `Quarter` object with the given `v` value. Input must be losslessly convertible to an Int64.
"""
struct Quarter <: DatePeriod
	value::Int64
    Quarter(v::Number) = new(v)
end

#arithmetic
Base.:+(x::Date, y::Quarter) = x + Month(y)
Base.:-(x::Date, y::Quarter) = x - Month(y)
Base.:+(x::DateTime, y::Quarter) = x + Month(y)
Base.:-(x::DateTime, y::Quarter) = x - Month(y)
#periods
let vmax = typemax(Int64) ÷ 4, vmin = typemin(Int64) ÷ 4
    @eval function Base.convert(::Type{Quarter}, x::Year)
        $vmin ≤ value(x) ≤ $vmax || throw(InexactError(:convert, Quarter, x))
        Quarter(value(x) * 4)
    end
end
Base.convert(::Type{Year}, x::Quarter) = Year(divexact(value(x), 4))
Base.promote_rule(::Type{Year}, ::Type{Quarter}) = Quarter

let vmax = typemax(Int64) ÷ 3, vmin = typemin(Int64) ÷ 3
    @eval function Base.convert(::Type{Month}, x::Quarter)
        $vmin ≤ value(x) ≤ $vmax || throw(InexactError(:convert, Month, x))
        Month(value(x) * 3)
    end
end
Base.convert(::Type{Quarter}, x::Month) = Quarter(Dates.divexact(value(x), 3))
Base.promote_rule(::Type{Quarter}, ::Type{Month}) = Month
Base.hash(x::Quarter, h::UInt) = hash(Month(x), h)
Dates.toms(x::Quarter) = Dates.toms(Month(x))
Dates.days(x::Quarter) = Dates.days(Month(x))
Dates.periodisless(::Year, ::Quarter) = false
Dates.periodisless(::Period, ::Quarter) = true
Dates.periodisless(::Quarter, ::Month) = false
# rounding
Base.floor(dt::Date, p::Quarter) = floor(dt, Month(p))

function Base.print(io::IO, q::Quarter)
    if value(q) == 1
        print(io, "1 quarter")
    else
        print(io, "$(value(q)) quarters")
    end
end
Base.show(io::IO, q::Quarter) = print(io, q)

"""
`Quarter(dt::TimeType) -> Quarter`

The quarter part of a `TimeType` as a `Quarter`
"""
Quarter(dt::Date) = Quarter(quarterofyear(dt))
Quarter(dt::DateTime) = Quarter(quarterofyear(dt))

export Quarter


text:
#----------------------------------------
# @__{slide()}__
# @__{theme("simple")}__
# @__{header("""
Demo Slides for jl2nb
""")}__
# @__{footer("""
Footer
""")}__
#
#
# Title
# =====
#----------------------------------------
# @__{slide()}__
#
# ## This is a slide with  some subslides
#
# It contains an executable cell with Julia code
println("Hello World");
#----
# @__{subslide()}__
# ### This is a sub-slide, again with executable code
println("Hello World from the subslide");
#----------------------------------------
# @__{slide()}__
#=

## Another slide with fragments

Here, we use block comments ( `#=` and `=#` must be placed in the first column),
and we can nicely use LaTeX formulas.

\begin{pmatrix}
  \int_a^b sin(x^2) dx & X\cap\emptyset=\emptyset & X\cup\emptyset=X\\
  \Rightarrow& x^{y^z}& \sin(2\pi)
\end{pmatrix}

=#
#-
# @__{fragment()}__
# - First item
#-
# @__{fragment()}__
# - Second item
#-
# @__{fragment()}__
# - Third item


text:
#
# Boundaries and Boundary Conditions
#

const coordinates = (:x, :y, :z)
const dims = length(coordinates)
const solution_fields = (:u, :v, :w, :T, :S)
const nsolution = length(solution_fields)

abstract type BCType end
struct Default <: BCType end

abstract type NonDefault <: BCType end
struct Flux <: NonDefault end
struct Gradient <: NonDefault end
struct Value <: NonDefault end

"""
    BoundaryCondition(BCType, condition)

Construct a boundary condition of `BCType` with `condition`,
where `BCType` is `Flux` or `Gradient`. `condition` may be a
number, array, or a function with signature:

    condition(t, Δx, Δy, Δz, Nx, Ny, Nz, u, v, w, T, S, iteration, i, j) = # function definition

`i` and `j` are indices along the boundary.
"""
struct BoundaryCondition{C<:BCType, T}
    condition::T
end

# Constructors
BoundaryCondition(Tbc, c) = BoundaryCondition{Tbc, typeof(c)}(c)

Adapt.adapt_structure(to, b::BoundaryCondition{C, A}) where {C<:BCType, A<:AbstractArray} = BoundaryCondition(C, Adapt.adapt(to, parent(b.condition)))

"""
    TimeVaryingBoundaryCondition(T, f)

Construct a boundary condition f(t) as a function of time only.
"""
function TimeVaryingBoundaryCondition(T::BCType, f)
  @inline condition(t, args...) = f(t)
  BoundaryCondition{T, Function}(condition)
end

DefaultBC() = BoundaryCondition{Default, Nothing}(nothing)

"""
    CoordinateBoundaryConditions(c)

Construct `CoordinateBoundaryCondition` to be applied along coordinate `c`, where
`c` is `:x`, `:y`, or `:z`. A CoordinateBoundaryCondition has two fields
`left` and `right` that store boundary conditions on the 'left' (negative side)
and 'right' (positive side) of a given coordinate.
"""
mutable struct CoordinateBoundaryConditions
  left::BoundaryCondition
  right::BoundaryCondition
end

CoordinateBoundaryConditions() = CoordinateBoundaryConditions(DefaultBC(), DefaultBC())

"""
    FieldBoundaryConditions()

Construct `FieldBoundaryConditions` for a field.
A FieldBoundaryCondition has `CoordinateBoundaryConditions` in
`x`, `y`, and `z`.
"""
struct FieldBoundaryConditions <: FieldVector{dims, CoordinateBoundaryConditions}
  x::CoordinateBoundaryConditions
  y::CoordinateBoundaryConditions
  z::CoordinateBoundaryConditions
end

"""
    BoundaryConditions()

Construct a boundary condition type full of default
`FieldBoundaryConditions` for u, v, w, T, S.
"""
struct ModelBoundaryConditions <: FieldVector{nsolution, FieldBoundaryConditions}
  u::FieldBoundaryConditions
  v::FieldBoundaryConditions
  w::FieldBoundaryConditions
  T::FieldBoundaryConditions
  S::FieldBoundaryConditions
end

FieldBoundaryConditions() = FieldBoundaryConditions(CoordinateBoundaryConditions(),
                                                    CoordinateBoundaryConditions(),
                                                    CoordinateBoundaryConditions())

function ModelBoundaryConditions()
  bcs = (FieldBoundaryConditions() for i = 1:length(solution_fields))
  return ModelBoundaryConditions(bcs...)
end

#
# User API
#
# Note:
#
# The syntax model.boundary_conditions.u.x.left = bc works, out of the box.
# How can we make it easier for users to set boundary conditions?

const BC = BoundaryCondition
const FBCs = FieldBoundaryConditions

#
# Physics goes here.
#

#=
Currently we support flux and gradient boundary conditions
at the top and bottom of the domain.

Notes:

- The boundary condition on a z-boundary is a callable object with arguments

      (t, Δx, Δy, Δz, Nx, Ny, Nz, u, v, w, T, S, iteration, i, j),

  where i and j are the x and y indices, respectively. No other function signature will work.
  We do not have abstractions that generalize to non-uniform grids.

- We assume that the boundary tendency has been previously calculated for
  a 'no-flux' boundary condition.

  This means that boudnary conditions take the form of
  an addition/subtraction to the tendency associated with a flux at point (A, A, I) below the bottom cell.
  This paradigm holds as long as consider boundary conditions on (A, A, C) variables only, where A is
  "any" of C or I.

 - We use the physics-based convention that

        flux = -κ * gradient,

    and that

        tendency = ∂ϕ/∂t = Gϕ = - ∇ ⋅ flux

=#

# Do nothing in default case. These functions are called in cases where one of the
# z-boundaries is set, but not the other.
@inline apply_z_top_bc!(args...) = nothing
@inline apply_z_bottom_bc!(args...) = nothing

# These functions compute vertical fluxes for (A, A, C) quantities. They are not currently used.
@inline ∇κ∇ϕ_t(κ, ϕt, ϕt₋₁, flux, ΔzC, ΔzF) = (      -flux        - κ*(ϕt - ϕt₋₁)/ΔzC ) / ΔzF
@inline ∇κ∇ϕ_b(κ, ϕb, ϕb₊₁, flux, ΔzC, ΔzF) = ( κ*(ϕb₊₁ - ϕb)/ΔzC +       flux        ) / ΔzF

"Add flux divergence to ∂ϕ/∂t associated with a top boundary condition on ϕ."
@inline apply_z_top_bc!(top_flux::BC{<:Flux, <:Function}, ϕ, Gϕ, κ, t, grid, u, v, w, T, S, iteration, i, j) =
    Gϕ[i, j, 1] -= top_flux.condition(t, grid, u, v, w, T, S, iteration, i, j) / grid.Δz

@inline apply_z_top_bc!(top_flux::BC{<:Flux, <:Number}, ϕ, Gϕ, κ, t, grid, u, v, w, T, S, iteration, i, j) =
    Gϕ[i, j, 1] -= top_flux.condition / grid.Δz

@inline apply_z_top_bc!(top_flux::BC{<:Flux, <:AbstractArray}, ϕ, Gϕ, κ, t, grid, u, v, w, T, S, iteration, i, j) =
    Gϕ[i, j, 1] -= top_flux.condition[i, j] / grid.Δz

@inline apply_z_top_bc!(top_gradient::BC{<:Gradient, <:Function}, ϕ, Gϕ, κ, t, grid, u, v, w, T, S, iteration, i, j) =
    Gϕ[i, j, 1] += κ*top_gradient.condition(t, grid, u, v, w, T, S, iteration, i, j) / grid.Δz

@inline apply_z_top_bc!(top_gradient::BC{<:Gradient, <:Number}, ϕ, Gϕ, κ, t, grid, u, v, w, T, S, iteration, i, j) =
    Gϕ[i, j, 1] += κ*top_gradient.condition / grid.Δz

@inline apply_z_top_bc!(top_gradient::BC{<:Gradient, <:AbstractArray}, ϕ, Gϕ, κ, t, grid, u, v, w, T, S, iteration, i, j) =
    Gϕ[i, j, 1] += κ*top_gradient.condition[i, j] / grid.Δz

"Add flux divergence to ∂ϕ/∂t associated with a bottom boundary condition on ϕ."
@inline apply_z_bottom_bc!(bottom_flux::BC{<:Flux, <:Function}, ϕ, Gϕ, κ, t, grid, u, v, w, T, S, iteration, i, j) =
    Gϕ[i, j, grid.Nz] += bottom_flux.condition(t, grid, u, v, w, T, S, iteration, i, j) / grid.Δz

@inline apply_z_bottom_bc!(bottom_flux::BC{<:Flux, <:Number}, ϕ, Gϕ, κ, t, grid, u, v, w, T, S, iteration, i, j) =
    Gϕ[i, j, grid.Nz] += bottom_flux.condition / grid.Δz

@inline apply_z_bottom_bc!(bottom_flux::BC{<:Flux, <:AbstractArray}, ϕ, Gϕ, κ, t, grid, u, v, w, T, S, iteration, i, j) =
    Gϕ[i, j, grid.Nz] += bottom_flux.condition[i, j] / grid.Δz

@inline apply_z_bottom_bc!(bottom_gradient::BC{<:Gradient, <:Function}, ϕ, Gϕ, κ, t, grid, u, v, w, T, S, iteration, i, j) =
    Gϕ[i, j, grid.Nz] -= κ*bottom_gradient.condition(t, grid, u, v, w, T, S, iteration, i, j) / grid.Δz

@inline apply_z_bottom_bc!(bottom_gradient::BC{<:Gradient, <:Number}, ϕ, Gϕ, κ, t, grid, u, v, w, T, S, iteration, i, j) =
    Gϕ[i, j, grid.Nz] -= κ*bottom_gradient.condition / grid.Δz

@inline apply_z_bottom_bc!(bottom_gradient::BC{<:Gradient, <:AbstractArray}, ϕ, Gϕ, κ, t, grid, u, v, w, T, S, iteration, i, j) =
    Gϕ[i, j, grid.Nz] -= κ*bottom_gradient.condition[i, j] / grid.Δz


text:
function Base.show(io::IO, m::CPLEX.Model)

    @printf "%0s \n" "CPLEX Model"
    @printf "   %-15s : %10s \n" "sense " CPLEX.get_sense(m) == :Max ? "maximize" : "minimize"
    @printf "   %-15s = %10s \n" "variables" CPLEX.num_var(m)
    @printf "   %-15s = %10s \n" "constraints" CPLEX.num_constr(m)
    @printf "   %-15s = %10s \n" "coefficients" CPLEX.get_nnz(m)
end

# To change settings for CPLEX, find the parameter you wish
# to change, uncomment its line, and change the X to the desired
# value.
# For quick (necessarily) reference for parameters, see
# http://www-eio.upc.edu/lceio/manuals/cplex-11/pdf/refparameterscplex.pdf
# for detailed specification of parameters
#
# Note: You can also, instead of numbers, use fx
#
# CPLEX.set_param!(env, "CPX_PARAM_SCRIND", CPX_ON)
#
# CPLEX.set_param!(env, "CPX_PARAM_SCRIND", 1)
#
function settings_cplex()
    env = CPLEX.Env()

    CPLEX.set_param!(env, "CPX_PARAM_LPMETHOD", 1) # 0=automatic
    CPLEX.set_param!(env, "CPX_PARAM_SCRIND", 1) # 1=display messages
    CPLEX.set_param!(env, "CPX_PARAM_SIMDISPLAY", 0) # 0=simplex msg off

    #
        #CPLEX.set_param!(env, CPX_PARAM_ADVIND, X)
        #CPLEX.set_param!(env, CPX_PARAM_AGGCUTLIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_AGGFILL, X)
        #CPLEX.set_param!(env, CPX_PARAM_AGGIND, X)
        #CPLEX.set_param!(env, CPX_PARAM_ALL_MAX, X)
        #CPLEX.set_param!(env, CPX_PARAM_ALL_MIN, X)
        #CPLEX.set_param!(env, CPX_PARAM_APIENCODING, X)
        #CPLEX.set_param!(env, CPX_PARAM_BARALG, X)
        #CPLEX.set_param!(env, CPX_PARAM_BARCOLNZ, X)
        #CPLEX.set_param!(env, CPX_PARAM_BARCROSSALG, X)
        #CPLEX.set_param!(env, CPX_PARAM_BARDISPLAY, X)
        #CPLEX.set_param!(env, CPX_PARAM_BAREPCOMP, X)
        #CPLEX.set_param!(env, CPX_PARAM_BARGROWTH, X)
        #CPLEX.set_param!(env, CPX_PARAM_BARITLIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_BARMAXCOR, X)
        #CPLEX.set_param!(env, CPX_PARAM_BAROBJRNG, X)
        #CPLEX.set_param!(env, CPX_PARAM_BARORDER, X)
        #CPLEX.set_param!(env, CPX_PARAM_BARQCPEPCOMP, X)
        #CPLEX.set_param!(env, CPX_PARAM_BARSTARTALG, X)
        #CPLEX.set_param!(env, CPX_PARAM_BBINTERVAL, X)
        #CPLEX.set_param!(env, CPX_PARAM_BNDSTRENIND, X)
        #CPLEX.set_param!(env, CPX_PARAM_BQPCUTS, X)
        #CPLEX.set_param!(env, CPX_PARAM_BRDIR, X)
        #CPLEX.set_param!(env, CPX_PARAM_BTTOL, X)
        #CPLEX.set_param!(env, CPX_PARAM_CALCQCPDUALS, X)
        #CPLEX.set_param!(env, CPX_PARAM_CLIQUES, X)
        #CPLEX.set_param!(env, CPX_PARAM_CLOCKTYPE, X)
        #CPLEX.set_param!(env, CPX_PARAM_CLONELOG, X)
        #CPLEX.set_param!(env, CPX_PARAM_COEREDIND, X)
        #CPLEX.set_param!(env, CPX_PARAM_COLREADLIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_CONFLICTALG, X)
        #CPLEX.set_param!(env, CPX_PARAM_CONFLICTDISPLAY, X)
        #CPLEX.set_param!(env, CPX_PARAM_COVERS, X)
        #CPLEX.set_param!(env, CPX_PARAM_CPUMASK, X)
        #CPLEX.set_param!(env, CPX_PARAM_CRAIND, X)
        #CPLEX.set_param!(env, CPX_PARAM_CUTLO, X)
        #CPLEX.set_param!(env, CPX_PARAM_CUTPASS, X)
        #CPLEX.set_param!(env, CPX_PARAM_CUTSFACTOR, X)
        #CPLEX.set_param!(env, CPX_PARAM_CUTUP, X)
        #CPLEX.set_param!(env, CPX_PARAM_DATACHECK, X)
        #CPLEX.set_param!(env, CPX_PARAM_DEPIND, X)
        #CPLEX.set_param!(env, CPX_PARAM_DETTILIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_DISJCUTS, X)
        #CPLEX.set_param!(env, CPX_PARAM_DIVETYPE, X)
        #CPLEX.set_param!(env, CPX_PARAM_DPRIIND, X)
        #CPLEX.set_param!(env, CPX_PARAM_EACHCUTLIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_EPAGAP, X)
        #CPLEX.set_param!(env, CPX_PARAM_EPGAP, X)
        #CPLEX.set_param!(env, CPX_PARAM_EPINT, X)
        #CPLEX.set_param!(env, CPX_PARAM_EPLIN, X)
        #CPLEX.set_param!(env, CPX_PARAM_EPMRK, X)
        #CPLEX.set_param!(env, CPX_PARAM_EPOPT, X)
        #CPLEX.set_param!(env, CPX_PARAM_EPPER, X)
        #CPLEX.set_param!(env, CPX_PARAM_EPRELAX, X)
        #CPLEX.set_param!(env, CPX_PARAM_EPRHS, X)
        #CPLEX.set_param!(env, CPX_PARAM_FEASOPTMODE, X)
        #CPLEX.set_param!(env, CPX_PARAM_FILEENCODING, X)
        #CPLEX.set_param!(env, CPX_PARAM_FLOWCOVERS, X)
        #CPLEX.set_param!(env, CPX_PARAM_FLOWPATHS, X)
        #CPLEX.set_param!(env, CPX_PARAM_FPHEUR, X)
        #CPLEX.set_param!(env, CPX_PARAM_FRACCAND, X)
        #CPLEX.set_param!(env, CPX_PARAM_FRACCUTS, X)
        #CPLEX.set_param!(env, CPX_PARAM_FRACPASS, X)
        #CPLEX.set_param!(env, CPX_PARAM_GUBCOVERS, X)
        #CPLEX.set_param!(env, CPX_PARAM_HEURFREQ, X)
        #CPLEX.set_param!(env, CPX_PARAM_IMPLBD, X)
        #CPLEX.set_param!(env, CPX_PARAM_INTSOLFILEPREFIX, X)
        #CPLEX.set_param!(env, CPX_PARAM_INTSOLLIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_ITLIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_LANDPCUTS, X)
        #CPLEX.set_param!(env, CPX_PARAM_LBHEUR, X)
        #CPLEX.set_param!(env, CPX_PARAM_LOCALIMPLBD, X)
        #CPLEX.set_param!(env, CPX_PARAM_LPMETHOD, X)
        #CPLEX.set_param!(env, CPX_PARAM_MCFCUTS, X)
        #CPLEX.set_param!(env, CPX_PARAM_MEMORYEMPHASIS, X)
        #CPLEX.set_param!(env, CPX_PARAM_MIPCBREDLP, X)
        #CPLEX.set_param!(env, CPX_PARAM_MIPDISPLAY, X)
        #CPLEX.set_param!(env, CPX_PARAM_MIPEMPHASIS, X)
        #CPLEX.set_param!(env, CPX_PARAM_MIPINTERVAL, X)
        #CPLEX.set_param!(env, CPX_PARAM_MIPKAPPASTATS, X)
        #CPLEX.set_param!(env, CPX_PARAM_MIPORDIND, X)
        #CPLEX.set_param!(env, CPX_PARAM_MIPORDTYPE, X)
        #CPLEX.set_param!(env, CPX_PARAM_MIPSEARCH, X)
        #CPLEX.set_param!(env, CPX_PARAM_MPSLONGNUM, X)
        #CPLEX.set_param!(env, CPX_PARAM_NETDISPLAY, X)
        #CPLEX.set_param!(env, CPX_PARAM_NETEPOPT, X)
        #CPLEX.set_param!(env, CPX_PARAM_NETEPRHS, X)
        #CPLEX.set_param!(env, CPX_PARAM_NETFIND, X)
        #CPLEX.set_param!(env, CPX_PARAM_NETITLIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_NETPPRIIND, X)
        #CPLEX.set_param!(env, CPX_PARAM_NODEFILEIND, X)
        #CPLEX.set_param!(env, CPX_PARAM_NODELIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_NODESEL, X)
        #CPLEX.set_param!(env, CPX_PARAM_NUMERICALEMPHASIS, X)
        #CPLEX.set_param!(env, CPX_PARAM_NZREADLIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_OBJDIF, X)
        #CPLEX.set_param!(env, CPX_PARAM_OBJLLIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_OBJULIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_OPTIMALITYTARGET, X)
        #CPLEX.set_param!(env, CPX_PARAM_PARALLELMODE, X)
        #CPLEX.set_param!(env, CPX_PARAM_PERIND, X)
        #CPLEX.set_param!(env, CPX_PARAM_PERLIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_POLISHAFTERDETTIME, X)
        #CPLEX.set_param!(env, CPX_PARAM_POLISHAFTEREPAGAP, X)
        #CPLEX.set_param!(env, CPX_PARAM_POLISHAFTEREPGAP, X)
        #CPLEX.set_param!(env, CPX_PARAM_POLISHAFTERINTSOL, X)
        #CPLEX.set_param!(env, CPX_PARAM_POLISHAFTERNODE, X)
        #CPLEX.set_param!(env, CPX_PARAM_POLISHAFTERTIME, X)
        #CPLEX.set_param!(env, CPX_PARAM_POLISHTIME, X)
        #CPLEX.set_param!(env, CPX_PARAM_POPULATELIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_PPRIIND, X)
        #CPLEX.set_param!(env, CPX_PARAM_PREDUAL, X)
        #CPLEX.set_param!(env, CPX_PARAM_PREIND, X)
        #CPLEX.set_param!(env, CPX_PARAM_PREPASS, X)
        #CPLEX.set_param!(env, CPX_PARAM_PROBETIME, X)
        #CPLEX.set_param!(env, CPX_PARAM_QPMAKEPSDIND, X)
        #CPLEX.set_param!(env, CPX_PARAM_QPMETHOD, X)
        #CPLEX.set_param!(env, CPX_PARAM_QPNZREADLIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_QTOLININD, X)
        #CPLEX.set_param!(env, CPX_PARAM_RAMPUPDETTILIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_RAMPUPDURATION, X)
        #CPLEX.set_param!(env, CPX_PARAM_RAMPUPTILIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_REINV, X)
        #CPLEX.set_param!(env, CPX_PARAM_RELAXPREIND, X)
        #CPLEX.set_param!(env, CPX_PARAM_RELOBJDIF, X)
        #CPLEX.set_param!(env, CPX_PARAM_REPAIRTRIES, X)
        #CPLEX.set_param!(env, CPX_PARAM_RINSHEUR, X)
        #CPLEX.set_param!(env, CPX_PARAM_ROWREADLIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_SCAIND, X)
        #CPLEX.set_param!(env, CPX_PARAM_SCRIND, X)
        #CPLEX.set_param!(env, CPX_PARAM_SIFTALG, X)
        #CPLEX.set_param!(env, CPX_PARAM_SIFTDISPLAY, X)
        #CPLEX.set_param!(env, CPX_PARAM_SIFTITLIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_SIMDISPLAY, X)
        #CPLEX.set_param!(env, CPX_PARAM_SINGLIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_SOLNPOOLAGAP, X)
        #CPLEX.set_param!(env, CPX_PARAM_SOLNPOOLGAP, X)
        #CPLEX.set_param!(env, CPX_PARAM_SOLUTIONTYPE, X)
        #CPLEX.set_param!(env, CPX_PARAM_STARTALG, X)
        #CPLEX.set_param!(env, CPX_PARAM_STRONGITLIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_SUBALG, X)
        #CPLEX.set_param!(env, CPX_PARAM_SYMMETRY, X)
        #CPLEX.set_param!(env, CPX_PARAM_THREADS, X)
        #CPLEX.set_param!(env, CPX_PARAM_TILIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_TRELIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_TUNINGREPEAT, X)
        #CPLEX.set_param!(env, CPX_PARAM_TUNINGTILIM, X)
        #CPLEX.set_param!(env, CPX_PARAM_VARSEL, X)
        #CPLEX.set_param!(env, CPX_PARAM_WORKDIR, X)
        #CPLEX.set_param!(env, CPX_PARAM_WORKMEM, X)
        #CPLEX.set_param!(env, CPX_PARAM_WRITELEVEL, X)
        #CPLEX.set_param!(env, CPX_PARAM_ZEROHALFCUTS, X)
    return env
end

function change_objective_coef(lp::CPLEX.Model, objective::Number, value::Number = 1.0)
	num_colums = CPLEX.num_var(lp)

	new_objective = zeros(num_colums)

	new_objective[objective] = value

	CPLEX.set_obj!(lp, new_objective)
end

function change_objective_coef{T <: Number}(lp::CPLEX.Model, objective::Array{T}, value::Array = [])
	if value == []
		value = ones(Float64, length(objective))
	end

	num_colums = CPLEX.num_var(lp)

	new_objective = zeros(num_colums)

	map((x,y) -> new_objective[x] = y, objective,value)

	CPLEX.set_obj!(lp, new_objective)
end

function get_variable_bounds(lp::CPLEX.Model)
	lb = CPLEX.get_varLB(lp)
	ub = CPLEX.get_varUB(lp)

	return lb, ub
end

function get_variable_bounds(lp::CPLEX.Model, index::Number)
	lb = CPLEX.get_varLB(lp)[index]
	ub = CPLEX.get_varUB(lp)[index]
	return lb,ub
end

function get_constraint_bounds(lp::CPLEX.Model)
	row_lb = CPLEX.get_constrLB(lp)
	row_ub = CPLEX.get_constrUB(lp)

	return row_lb, row_ub
end

function get_constraint_matrix(lp::CPLEX.Model)
	return CPLEX.get_constr_matrix(lp)
end

function get_direction(lp::CPLEX.Model)
	direction = CPLEX.get_sense(lp)

	if CPLEX.get_sense(lp) == :Max
		return "maximize"
	else
		return "minimize"
	end
end

function get_objective_value(lp::CPLEX.Model)
	return CPLEX.get_objval(lp)
end

function get_reduced_costs(lp::CPLEX.Model)
	return CPLEX.get_reduced_costs(lp)
end

function get_solution(lp::CPLEX.Model)
	return CPLEX.get_solution(lp)
end

function get_solution_status(lp::CPLEX.Model)
	code = get_solution_status_code(lp)
	if code == 1 #1 CPX_STAT_OPTIMAL
		return "Optimal"
	elseif code == 4 #4 CPX_STAT_INForUNBD
		return "No feasible solution"
	elseif code == 3 #3 CPX_STAT_INFEASIBLE
		return "Infeasible"
	elseif code == 5 #5 CPX_STAT_OPTIMAL_INFEAS
		return "Feasible"
	elseif code == 2 #2 CPX_STAT_UNBOUNDED
		return "Unbound"
	elseif code == 6 #6 CPX_STAT_NUM_BEST
		return "Undefined"
	else
		return warn(string("Unknown status code: ", code))
	end

	#7 CPX_STAT_FEASIBLE_RELAXED
	#8 CPX_STAT_OPTIMAL_RELAXED
	#10 CPX_STAT_ABORT_IT_LIM
	#11 CPX_STAT_ABORT_TIME_LIM
	#12 CPX_STAT_ABORT_OBJ_LIM
	#13 CPX_STAT_ABORT_USER
	#20 CPX_STAT_OPTIMAL_FACE_UNBOUNDED
end

function get_solution_status_code(lp::CPLEX.Model)
	return getfield(CPLEX, CPLEX.get_status(lp))
end

function get_slack(lp::CPLEX.Model)
	return CPLEX.get_constr_duals(lp)
end

function get_objective(lp::CPLEX.Model)
	return CPLEX.get_obj(lp)
end

function set_direction(lp::CPLEX.Model, direction::AbstractString)
	if lowercase(direction) == "max"
		CPLEX.set_sense!(lp, :Max)
	else
		CPLEX.set_sense!(lp, :Min)
	end
end

function set_col_bounds(lp::CPLEX.Model, index::Number, lb::Number, ub::Number)
	new_lb = CPLEX.get_varLB(lp)
	new_lb[index] = lb
	CPLEX.set_varLB!(lp, new_lb)

	new_ub = CPLEX.get_varUB(lp)
	new_ub[index] = ub
	CPLEX.set_varUB!(lp, new_ub)
end

function set_col_bounds(lp::CPLEX.Model, index::Number, fixed_value::Number)
	new_lb = CPLEX.get_varLB(lp)
	new_lb[index] = fixed_value
	CPLEX.set_varLB!(lp, new_lb)

	new_ub = CPLEX.get_varUB(lp)
	new_ub[index] = fixed_value
	CPLEX.set_varUB!(lp, new_ub)
end

function set_objective(lp::CPLEX.Model, objective::Number, value::Number = 1.0)
	new_objective = CPLEX.get_obj(lp)
	new_objective[objective] = value

	CPLEX.set_obj!(lp, new_objective)
end

function setup_cplex(model::Model, objective::AbstractString = "max")
  # Flux Balance Analysis usingh CPLEX.jl
	nrxns = length(model.c)
	nmets = length(model.b)

	# Construct the CPLEX model
	env = settings_cplex()
	lp = CPLEX.Model(env, "CPLEX")

	if objective == "max"
		set_direction(lp, "max")
	else
		set_direction(lp, "min")
	end

    CPLEX.add_vars!(lp, model.c, model.lb, model.ub)

    ####
    if length(model.b) == length(model.csense["="])
	   CPLEX.add_constrs!(lp, model.S, fill('=', nmets), model.b)
    else
        b_lower = deepcopy(model.b)
        b_upper = deepcopy(model.b)
        b_lower[get(model.csense,"<=", [])] = -Inf
        b_upper[get(model.csense,">=", [])] = Inf
        CPLEX.add_rangeconstrs!(lp, model.S, b_lower, b_upper)
    end
    ####

	return lp
end

function setup_cplex{T <: Number}(c::Array{T}, lb::Array{T}, ub::Array{T}, b::Array{T}, S::SparseMatrixCSC, objective::AbstractString = "max")
  # Flux Balance Analysis usingh CPLEX.jl

	nmets = length(b)

	# Construct the CPLEX model
	env = CPLEX.Env()
	lp = CPLEX.Model(env, "CPLEX")
	CPLEX.set_param!(env, "CPX_PARAM_LPMETHOD", CPLEX.CPX_ALG_AUTOMATIC) # should?

	if objective == "max"
		set_direction(lp, "max")
	else
		set_direction(lp, "min")
	end

	CPLEX.add_vars!(lp, c, lb, ub)

	CPLEX.add_constrs!(lp, S, fill('=', nmets), b)

	return lp
end

function setup_cplex{T <: Number}(c::Array{T}, lb::Array{T}, ub::Array{T}, b_lower::Array{T}, b_upper::Array{T}, S::SparseMatrixCSC, objective::AbstractString = "max")
    nmets = length(b_lower)

    # Construct the CPLEX model
    env = CPLEX.Env()
    lp = CPLEX.Model(env, "CPLEX")
    CPLEX.set_param!(env, "CPX_PARAM_LPMETHOD", CPLEX.CPX_ALG_AUTOMATIC) # should?

    if objective == "max"
        set_direction(lp, "max")
    else
        set_direction(lp, "min")
    end

    CPLEX.add_vars!(lp, c, lb, ub)

    CPLEX.add_rangeconstrs!(lp, S, b_lower, b_upper)

    return lp
end

function setup_cplex{T <: Number}(c::Array{T}, lb::Array{T}, ub::Array{T}, b::Array{T}, csense::Array{String}, S::SparseMatrixCSC, objective::AbstractString = "max")
    b_lower = deepcopy(b)
    b_upper = deepcopy(b)
    b_lower[find(csense .== "<=")] = -Inf
    b_upper[find(csense .== ">=")] = Inf

    setup_cplex(c, lb, ub, b_lower, b_upper, S, objective)
end

function setup_cplex{T <: Number}(c::Array{T}, lb::Array{T}, ub::Array{T}, b::Array{T}, csense::Dict, S::SparseMatrixCSC, objective::AbstractString = "max")
    b_lower = deepcopy(b)
    b_upper = deepcopy(b)
    b_lower[get(csense,"<=", [])] = -Inf
    b_upper[get(csense,">=", [])] = Inf

    setup_cplex(c, lb, ub, b_lower, b_upper, S, objective)

end

function setup_cplex{T <: Number}(c::Array{T}, lb::Array{T}, ub::Array{T}, beq::Array{T}, Seq::SparseMatrixCSC, bineq::Array{T}, Sineq::SparseMatrixCSC, objective::AbstractString = "max")
  # Flux Balance Analysis usingh CPLEX.jl

    # Construct the CPLEX model
    env = CPLEX.Env()
    lp = CPLEX.Model(env, "CPLEX")
    CPLEX.set_param!(env, "CPX_PARAM_LPMETHOD", CPLEX.CPX_ALG_AUTOMATIC) # should?

    if objective == "max"
        set_direction(lp, "max")
    else
        set_direction(lp, "min")
    end

    CPLEX.add_vars!(lp, c, lb, ub)

    CPLEX.add_rangeconstrs!(lp, Sineq, fill(-Inf, length(bineq)), bineq)
    CPLEX.add_constrs!(lp, Seq, fill('=', length(beq)), beq)

    return lp
end

function solve(lp::CPLEX.Model)
	CPLEX.optimize!(lp)
end

function solver_status(lp::CPLEX.Model)
    string(CPLEX.status_symbols[get_solution_status_code(lp)])
end

max_stars_repo_name: isebarn/CBM

max_issues_repo_name: isebarn/CBM

max_forks_repo_name: isebarn/CBM


text:
"""
    CBIOMES_combine_files()

Take all files generated by `CBIOMES_model_climatogy.jl` and combine them into one file.
"""
function CBIOMES_combine_files()

fil_out=joinpath(tempdir(),"CBIOMES_clim.nc")
list_in=["Chl","EuphoticDepth","MLD","OceanDepth","PAR",
    "Rrs412","Rrs443","Rrs490","Rrs510","Rrs555","Rrs670",
    "SSS","SST","TKE","WindSpeed"]

##

!isfile(fil_out) ? cp(joinpath(tempdir(),"SST_clim.nc"),fil_out) : nothing

for ii in list_in
    ds = Dataset(fil_out,"a")
    if !haskey(ds,ii)
        fil_in=joinpath(tempdir(),ii*"_clim.nc")
        ds_in = Dataset(fil_in,"r")
        tmp11=ds_in[ii][:]
        u=ds_in[ii].attrib["units"]
        ln=ds_in[ii].attrib["long_name"]
        close(ds_in)

        ##

        v = defVar(ds,ii,Float64,("lon","lat","t"),
        attrib = Dict("units" => u, "long_name" => ln))
        v[:] = tmp11
    end
    close(ds)
end

ds = Dataset(fil_out,"a")
ds.attrib["description,1"]="Source: Gael Forget"
ds.attrib["description,2"]="Product: CBIOMES-global climatology"
ds.attrib["description,3"]="Version: alpha"
if haskey(ds.attrib,"A")
    delete!(ds.attrib,"description")
    delete!(ds.attrib,"A")
    delete!(ds.attrib,"B")
end
close(ds)

"all set"
end


text:
include("line_scatter.jl")

function subplots1()
    p1 = linescatter1()
    p2 = linescatter2()
    p = [p1 p2]
    p
end

function subplots2()
    p1 = linescatter1()
    p2 = linescatter2()
    p = [p1, p2]
    p
end

function subplots3()
    p1 = linescatter6()
    p2 = linescatter2()
    p3 = linescatter3()
    p4 = batman()
    p = [p1 p2; p3 p4]
    p.plot.layout["showlegend"] = false
    p.plot.layout["width"] = 1000
    p.plot.layout["height"] = 600
    p
end

function subplots_withcomprehension()
    hcat([plot(scatter(x = 1:5, y = rand(5))) for i in 1:3]...)
end

function subplots_withsharedaxes()
    data =  [
    scatter(x=1:3, y=2:4),
    scatter(x=20:10:40, y=fill(5, 3), xaxis="x2", yaxis="y"),
    scatter(x=2:4, y=600:100:800, xaxis="x", yaxis="y3"),
    scatter(x=4000:1000:6000, y=7000:1000:9000, xaxis="x4", yaxis="y4")
    ]
    layout = Layout(
        xaxis_domain=[0, 0.45],
        yaxis_domain=[0, 0.45],
        xaxis4=attr(domain=[0.55, 1.0], anchor="y4"),
        xaxis2_domain=[0.55, 1],
        yaxis3_domain=[0.55, 1],
        yaxis4=attr(domain=[0.55, 1], anchor="x4")
    )
    plot(data, layout)
end