text:
"""
An AbstractDataFrame that stores a set of named columns

The columns are normally AbstractVectors stored in memory,
particularly a Vector, NullableVector, or CategoricalVector.

**Constructors**

```julia
DataFrame(columns::Vector{Any}, names::Vector{Symbol})
DataFrame(kwargs...)
DataFrame() # an empty DataFrame
DataFrame(t::Type, nrows::Integer, ncols::Integer) # an empty DataFrame of arbitrary size
DataFrame(column_eltypes::Vector, names::Vector, nrows::Integer)
DataFrame(ds::Vector{Associative})
```

**Arguments**

* `columns` : a Vector{Any} with each column as contents
* `names` : the column names
* `kwargs` : the key gives the column names, and the value is the
  column contents
* `t` : elemental type of all columns
* `nrows`, `ncols` : number of rows and columns
* `column_eltypes` : elemental type of each column
* `ds` : a vector of Associatives

Each column in `columns` should be the same length.

**Notes**

Most of the default constructors convert columns to `NullableArray`.  The
base constructor, `DataFrame(columns::Vector{Any},
names::Vector{Symbol})` does not convert to `NullableArray`.

A `DataFrame` is a lightweight object. As long as columns are not
manipulated, creation of a DataFrame from existing AbstractVectors is
inexpensive. For example, indexing on columns is inexpensive, but
indexing by rows is expensive because copies are made of each column.

Because column types can vary, a DataFrame is not type stable. For
performance-critical code, do not index into a DataFrame inside of
loops.

**Examples**

```julia
df = DataFrame()
v = ["x","y","z"][rand(1:3, 10)]
df1 = DataFrame(Any[[1:10], v, rand(10)], [:A, :B, :C])  # columns are Arrays
df2 = DataFrame(A = 1:10, B = v, C = rand(10))           # columns are NullableArrays
dump(df1)
dump(df2)
describe(df2)
DataFrames.head(df1)
df1[:A] + df2[:C]
df1[1:4, 1:2]
df1[[:A,:C]]
df1[1:2, [:A,:C]]
df1[:, [:A,:C]]
df1[:, [1,3]]
df1[1:4, :]
df1[1:4, :C]
df1[1:4, :C] = 40. * df1[1:4, :C]
[df1; df2]  # vcat
[df1  df2]  # hcat
size(df1)
```

"""
type DataFrame <: AbstractDataFrame
    columns::Vector{Any}
    colindex::Index

    function DataFrame(columns::Vector{Any}, colindex::Index)
        ncols = length(columns)
        if ncols > 1
            nrows = length(columns[1])
            equallengths = true
            for i in 2:ncols
                equallengths &= length(columns[i]) == nrows
            end
            if !equallengths
                msg = "All columns in a DataFrame must be the same length"
                throw(ArgumentError(msg))
            end
        end
        if length(colindex) != ncols
            msg = "Columns and column index must be the same length"
            throw(ArgumentError(msg))
        end
        new(columns, colindex)
    end
end

function DataFrame(; kwargs...)
    result = DataFrame(Any[], Index())
    for (k, v) in kwargs
        result[k] = v
    end
    return result
end

function DataFrame(columns::Vector{Any},
                   cnames::Vector{Symbol} = gennames(length(columns)))
    return DataFrame(columns, Index(cnames))
end

# Initialize empty DataFrame objects of arbitrary size
function DataFrame(t::Type, nrows::Integer, ncols::Integer)
    columns = Array(Any, ncols)
    for i in 1:ncols
        columns[i] = NullableArray(t, nrows)
    end
    cnames = gennames(ncols)
    return DataFrame(columns, Index(cnames))
end

# Initialize an empty DataFrame with specific eltypes and names
function DataFrame(column_eltypes::Vector, cnames::Vector, nrows::Integer)
    p = length(column_eltypes)
    columns = Array(Any, p)
    for j in 1:p
        columns[j] = NullableArray(column_eltypes[j], nrows)
    end
    return DataFrame(columns, Index(cnames))
end
# Initialize an empty DataFrame with specific eltypes and names
# and whether a nominal array should be created
function DataFrame(column_eltypes::Vector{DataType}, cnames::Vector{Symbol},
                   nominal::Vector{Bool}, nrows::Integer)
    p = length(column_eltypes)
    columns = Array(Any, p)
    for j in 1:p
      if nominal[j]
        columns[j] = NullableCategoricalArray(column_eltypes[j], nrows)
      else
        columns[j] = NullableArray(column_eltypes[j], nrows)
      end
    end
    return DataFrame(columns, Index(cnames))
end

# Initialize an empty DataFrame with specific eltypes
function DataFrame(column_eltypes::Vector, nrows::Integer)
    p = length(column_eltypes)
    columns = Array(Any, p)
    cnames = gennames(p)
    for j in 1:p
        columns[j] = NullableArray(column_eltypes[j], nrows)
    end
    return DataFrame(columns, Index(cnames))
end

# Initialize from a Vector of Associatives (aka list of dicts)
function DataFrame{D <: Associative}(ds::Vector{D})
    ks = Set()
    for d in ds
        union!(ks, keys(d))
    end
    DataFrame(ds, [ks...])
end

# Initialize from a Vector of Associatives (aka list of dicts)
function DataFrame{D <: Associative}(ds::Vector{D}, ks::Vector)
    #get column eltypes
    col_eltypes = Type[@compat(Union{}) for _ = 1:length(ks)]
    for d in ds
        for (i,k) in enumerate(ks)
            if haskey(d, k) && !_isnull(d[k])
                col_eltypes[i] = promote_type(col_eltypes[i], typeof(d[k]))
            end
        end
    end
    col_eltypes[col_eltypes .== @compat(Union{})] = Any

    # create empty DataFrame, and fill
    df = DataFrame(col_eltypes, ks, length(ds))
    for (i,d) in enumerate(ds)
        for (j,k) in enumerate(ks)
            df[i,j] = get(d, k, Nullable())
        end
    end

    df
end

##############################################################################
##
## AbstractDataFrame interface
##
##############################################################################

index(df::DataFrame) = df.colindex
columns(df::DataFrame) = df.columns

# TODO: Remove these
nrow(df::DataFrame) = ncol(df) > 0 ? length(df.columns[1])::Int : 0
ncol(df::DataFrame) = length(index(df))

##############################################################################
##
## getindex() definitions
##
##############################################################################

# Cases:
#
# df[SingleColumnIndex] => AbstractDataVector
# df[MultiColumnIndex] => (Sub)?DataFrame
# df[SingleRowIndex, SingleColumnIndex] => Scalar
# df[SingleRowIndex, MultiColumnIndex] => (Sub)?DataFrame
# df[MultiRowIndex, SingleColumnIndex] => (Sub)?AbstractDataVector
# df[MultiRowIndex, MultiColumnIndex] => (Sub)?DataFrame
#
# General Strategy:
#
# Let getindex(index(df), col_inds) from Index() handle the resolution
#  of column indices
# Let getindex(df.columns[j], row_inds) from AbstractDataVector() handle
#  the resolution of row indices

typealias ColumnIndex @compat(Union{Real, Symbol})

# df[SingleColumnIndex] => AbstractDataVector
function Base.getindex(df::DataFrame, col_ind::ColumnIndex)
    selected_column = index(df)[col_ind]
    return df.columns[selected_column]
end

# df[MultiColumnIndex] => (Sub)?DataFrame
function Base.getindex{T <: ColumnIndex}(df::DataFrame,
                                         col_inds::Union{AbstractVector{T},
                                                         AbstractVector{Nullable{T}}})
    selected_columns = index(df)[col_inds]
    new_columns = df.columns[selected_columns]
    return DataFrame(new_columns, Index(_names(df)[selected_columns]))
end

# df[:] => (Sub)?DataFrame
Base.getindex(df::DataFrame, col_inds::Colon) = copy(df)

# df[SingleRowIndex, SingleColumnIndex] => Scalar
function Base.getindex(df::DataFrame, row_ind::Real, col_ind::ColumnIndex)
    selected_column = index(df)[col_ind]
    return df.columns[selected_column][row_ind]
end

# df[SingleRowIndex, MultiColumnIndex] => (Sub)?DataFrame
function Base.getindex{T <: ColumnIndex}(df::DataFrame,
                                         row_ind::Real,
                                         col_inds::Union{AbstractVector{T},
                                                         AbstractVector{Nullable{T}}})
    selected_columns = index(df)[col_inds]
    new_columns = Any[dv[[row_ind]] for dv in df.columns[selected_columns]]
    return DataFrame(new_columns, Index(_names(df)[selected_columns]))
end

# df[MultiRowIndex, SingleColumnIndex] => (Sub)?AbstractDataVector
function Base.getindex{T <: Real}(df::DataFrame,
                                  row_inds::Union{AbstractVector{T}, AbstractVector{Nullable{T}}},
                                  col_ind::ColumnIndex)
    selected_column = index(df)[col_ind]
    return df.columns[selected_column][row_inds]
end

# df[MultiRowIndex, MultiColumnIndex] => (Sub)?DataFrame
function Base.getindex{R <: Real, T <: ColumnIndex}(df::DataFrame,
                                                    row_inds::Union{AbstractVector{R},
                                                                    AbstractVector{Nullable{R}}},
                                                    col_inds::Union{AbstractVector{T},
                                                                    AbstractVector{Nullable{T}}})
    selected_columns = index(df)[col_inds]
    new_columns = Any[dv[row_inds] for dv in df.columns[selected_columns]]
    return DataFrame(new_columns, Index(_names(df)[selected_columns]))
end

# df[:, SingleColumnIndex] => (Sub)?AbstractVector
# df[:, MultiColumnIndex] => (Sub)?DataFrame
Base.getindex{T<:ColumnIndex}(df::DataFrame,
                              row_inds::Colon,
                              col_inds::Union{T, AbstractVector{T},
                                              AbstractVector{Nullable{T}}}) =
    df[col_inds]

# df[SingleRowIndex, :] => (Sub)?DataFrame
Base.getindex(df::DataFrame, row_ind::Real, col_inds::Colon) = df[[row_ind], col_inds]

# df[MultiRowIndex, :] => (Sub)?DataFrame
function Base.getindex{R<:Real}(df::DataFrame,
                                row_inds::Union{AbstractVector{R},
                                                AbstractVector{Nullable{R}}},
                                col_inds::Colon)
    new_columns = Any[dv[row_inds] for dv in df.columns]
    return DataFrame(new_columns, copy(index(df)))
end

# df[:, :] => (Sub)?DataFrame
Base.getindex(df::DataFrame, ::Colon, ::Colon) = copy(df)

##############################################################################
##
## setindex!()
##
##############################################################################

isnextcol(df::DataFrame, col_ind::Symbol) = true
function isnextcol(df::DataFrame, col_ind::Real)
    return ncol(df) + 1 == @compat Int(col_ind)
end

function nextcolname(df::DataFrame)
    return @compat(Symbol(string("x", ncol(df) + 1)))
end

# Will automatically add a new column if needed
function insert_single_column!(df::DataFrame,
                               dv::AbstractVector,
                               col_ind::ColumnIndex)

    if ncol(df) != 0 && nrow(df) != length(dv)
        error("New columns must have the same length as old columns")
    end
    if haskey(index(df), col_ind)
        j = index(df)[col_ind]
        df.columns[j] = dv
    else
        if typeof(col_ind) <: Symbol
            push!(index(df), col_ind)
            push!(df.columns, dv)
        else
            if isnextcol(df, col_ind)
                push!(index(df), nextcolname(df))
                push!(df.columns, dv)
            else
                error("Cannot assign to non-existent column: $col_ind")
            end
        end
    end
    return dv
end

function insert_single_entry!(df::DataFrame, v::Any, row_ind::Real, col_ind::ColumnIndex)
    if haskey(index(df), col_ind)
        df.columns[index(df)[col_ind]][row_ind] = v
        return v
    else
        error("Cannot assign to non-existent column: $col_ind")
    end
end

function insert_multiple_entries!{T <: Real}(df::DataFrame,
                                            v::Any,
                                            row_inds::AbstractVector{T},
                                            col_ind::ColumnIndex)
    if haskey(index(df), col_ind)
        df.columns[index(df)[col_ind]][row_inds] = v
        return v
    else
        error("Cannot assign to non-existent column: $col_ind")
    end
end

upgrade_vector{T<:Nullable}(v::AbstractArray{T}) = v
upgrade_vector(v::CategoricalArray) = NullableCategoricalArray(v)
upgrade_vector(v::AbstractArray) = NullableArray(v)

function upgrade_scalar(df::DataFrame, v::AbstractArray)
    msg = "setindex!(::DataFrame, ...) only broadcasts scalars, not arrays"
    throw(ArgumentError(msg))
end
function upgrade_scalar(df::DataFrame, v::Any)
    n = (ncol(df) == 0) ? 1 : nrow(df)
    NullableArray(fill(v, n))
end

# df[SingleColumnIndex] = AbstractVector
function Base.setindex!(df::DataFrame,
                v::AbstractVector,
                col_ind::ColumnIndex)
    insert_single_column!(df, upgrade_vector(v), col_ind)
end

# df[SingleColumnIndex] = Single Item (EXPANDS TO NROW(DF) if NCOL(DF) > 0)
function Base.setindex!(df::DataFrame, v, col_ind::ColumnIndex)
    if haskey(index(df), col_ind)
        fill!(df[col_ind], v)
    else
        insert_single_column!(df, upgrade_scalar(df, v), col_ind)
    end
    return df
end

# df[MultiColumnIndex] = DataFrame
function Base.setindex!(df::DataFrame,
                new_df::DataFrame,
                col_inds::AbstractVector{Bool})
    setindex!(df, new_df, find(col_inds))
end
function Base.setindex!{T <: ColumnIndex}(df::DataFrame,
                                  new_df::DataFrame,
                                  col_inds::AbstractVector{T})
    for j in 1:length(col_inds)
        insert_single_column!(df, new_df[j], col_inds[j])
    end
    return df
end

# df[MultiColumnIndex] = AbstractVector (REPEATED FOR EACH COLUMN)
function Base.setindex!(df::DataFrame,
                v::AbstractVector,
                col_inds::AbstractVector{Bool})
    setindex!(df, v, find(col_inds))
end
function Base.setindex!{T <: ColumnIndex}(df::DataFrame,
                                  v::AbstractVector,
                                  col_inds::AbstractVector{T})
    dv = upgrade_vector(v)
    for col_ind in col_inds
        df[col_ind] = dv
    end
    return df
end

# df[MultiColumnIndex] = Single Item (REPEATED FOR EACH COLUMN; EXPANDS TO NROW(DF) if NCOL(DF) > 0)
function Base.setindex!(df::DataFrame,
                val::Any,
                col_inds::AbstractVector{Bool})
    setindex!(df, val, find(col_inds))
end
function Base.setindex!{T <: ColumnIndex}(df::DataFrame,
                                  val::Any,
                                  col_inds::AbstractVector{T})
    for col_ind in col_inds
        df[col_ind] = val
    end
    return df
end

# df[:] = AbstractVector or Single Item
Base.setindex!(df::DataFrame, v, ::Colon) = (df[1:size(df, 2)] = v; df)

# df[SingleRowIndex, SingleColumnIndex] = Single Item
function Base.setindex!(df::DataFrame,
                v::Any,
                row_ind::Real,
                col_ind::ColumnIndex)
    insert_single_entry!(df, v, row_ind, col_ind)
end

# df[SingleRowIndex, MultiColumnIndex] = Single Item
function Base.setindex!(df::DataFrame,
                v::Any,
                row_ind::Real,
                col_inds::AbstractVector{Bool})
    setindex!(df, v, row_ind, find(col_inds))
end
function Base.setindex!{T <: ColumnIndex}(df::DataFrame,
                                  v::Any,
                                  row_ind::Real,
                                  col_inds::AbstractVector{T})
    for col_ind in col_inds
        insert_single_entry!(df, v, row_ind, col_ind)
    end
    return df
end

# df[SingleRowIndex, MultiColumnIndex] = 1-Row DataFrame
function Base.setindex!(df::DataFrame,
                new_df::DataFrame,
                row_ind::Real,
                col_inds::AbstractVector{Bool})
    setindex!(df, new_df, row_ind, find(col_inds))
end
function Base.setindex!{T <: ColumnIndex}(df::DataFrame,
                                  new_df::DataFrame,
                                  row_ind::Real,
                                  col_inds::AbstractVector{T})
    for j in 1:length(col_inds)
        insert_single_entry!(df, new_df[j][1], row_ind, col_inds[j])
    end
    return df
end

# df[MultiRowIndex, SingleColumnIndex] = AbstractVector
function Base.setindex!(df::DataFrame,
                v::AbstractVector,
                row_inds::AbstractVector{Bool},
                col_ind::ColumnIndex)
    setindex!(df, v, find(row_inds), col_ind)
end
function Base.setindex!{T <: Real}(df::DataFrame,
                           v::AbstractVector,
                           row_inds::AbstractVector{T},
                           col_ind::ColumnIndex)
    insert_multiple_entries!(df, v, row_inds, col_ind)
    return df
end

# df[MultiRowIndex, SingleColumnIndex] = Single Item
function Base.setindex!(df::DataFrame,
                v::Any,
                row_inds::AbstractVector{Bool},
                col_ind::ColumnIndex)
    setindex!(df, v, find(row_inds), col_ind)
end
function Base.setindex!{T <: Real}(df::DataFrame,
                           v::Any,
                           row_inds::AbstractVector{T},
                           col_ind::ColumnIndex)
    insert_multiple_entries!(df, v, row_inds, col_ind)
    return df
end

# df[MultiRowIndex, MultiColumnIndex] = DataFrame
function Base.setindex!(df::DataFrame,
                new_df::DataFrame,
                row_inds::AbstractVector{Bool},
                col_inds::AbstractVector{Bool})
    setindex!(df, new_df, find(row_inds), find(col_inds))
end
function Base.setindex!{T <: ColumnIndex}(df::DataFrame,
                                  new_df::DataFrame,
                                  row_inds::AbstractVector{Bool},
                                  col_inds::AbstractVector{T})
    setindex!(df, new_df, find(row_inds), col_inds)
end
function Base.setindex!{R <: Real}(df::DataFrame,
                           new_df::DataFrame,
                           row_inds::AbstractVector{R},
                           col_inds::AbstractVector{Bool})
    setindex!(df, new_df, row_inds, find(col_inds))
end
function Base.setindex!{R <: Real, T <: ColumnIndex}(df::DataFrame,
                                             new_df::DataFrame,
                                             row_inds::AbstractVector{R},
                                             col_inds::AbstractVector{T})
    for j in 1:length(col_inds)
        insert_multiple_entries!(df, new_df[:, j], row_inds, col_inds[j])
    end
    return df
end

# df[MultiRowIndex, MultiColumnIndex] = AbstractVector
function Base.setindex!(df::DataFrame,
                v::AbstractVector,
                row_inds::AbstractVector{Bool},
                col_inds::AbstractVector{Bool})
    setindex!(df, v, find(row_inds), find(col_inds))
end
function Base.setindex!{T <: ColumnIndex}(df::DataFrame,
                                  v::AbstractVector,
                                  row_inds::AbstractVector{Bool},
                                  col_inds::AbstractVector{T})
    setindex!(df, v, find(row_inds), col_inds)
end
function Base.setindex!{R <: Real}(df::DataFrame,
                           v::AbstractVector,
                           row_inds::AbstractVector{R},
                           col_inds::AbstractVector{Bool})
    setindex!(df, v, row_inds, find(col_inds))
end
function Base.setindex!{R <: Real, T <: ColumnIndex}(df::DataFrame,
                                             v::AbstractVector,
                                             row_inds::AbstractVector{R},
                                             col_inds::AbstractVector{T})
    for col_ind in col_inds
        insert_multiple_entries!(df, v, row_inds, col_ind)
    end
    return df
end

# df[MultiRowIndex, MultiColumnIndex] = Single Item
function Base.setindex!(df::DataFrame,
                v::Any,
                row_inds::AbstractVector{Bool},
                col_inds::AbstractVector{Bool})
    setindex!(df, v, find(row_inds), find(col_inds))
end
function Base.setindex!{T <: ColumnIndex}(df::DataFrame,
                                  v::Any,
                                  row_inds::AbstractVector{Bool},
                                  col_inds::AbstractVector{T})
    setindex!(df, v, find(row_inds), col_inds)
end
function Base.setindex!{R <: Real}(df::DataFrame,
                           v::Any,
                           row_inds::AbstractVector{R},
                           col_inds::AbstractVector{Bool})
    setindex!(df, v, row_inds, find(col_inds))
end
function Base.setindex!{R <: Real, T <: ColumnIndex}(df::DataFrame,
                                             v::Any,
                                             row_inds::AbstractVector{R},
                                             col_inds::AbstractVector{T})
    for col_ind in col_inds
        insert_multiple_entries!(df, v, row_inds, col_ind)
    end
    return df
end

# df[:] = DataFrame, df[:, :] = DataFrame
function Base.setindex!(df::DataFrame,
                                  new_df::DataFrame,
                                  row_inds::Colon,
                                  col_inds::Colon=Colon())
    df.columns = copy(new_df.columns)
    df.colindex = copy(new_df.colindex)
    df
end

# df[:, :] = ...
Base.setindex!(df::DataFrame, v, ::Colon, ::Colon) =
    (df[1:size(df, 1), 1:size(df, 2)] = v; df)

# df[Any, :] = ...
Base.setindex!(df::DataFrame, v, row_inds, ::Colon) =
    (df[row_inds, 1:size(df, 2)] = v; df)

# df[:, Any] = ...
Base.setindex!(df::DataFrame, v, ::Colon, col_inds) =
    (df[col_inds] = v; df)

# Special deletion assignment
Base.setindex!(df::DataFrame, x::Void, col_ind::Int) = delete!(df, col_ind)

##############################################################################
##
## Mutating Associative methods
##
##############################################################################

Base.empty!(df::DataFrame) = (empty!(df.columns); empty!(index(df)); df)

function Base.insert!(df::DataFrame, col_ind::Int, item::AbstractVector, name::Symbol)
    0 < col_ind <= ncol(df) + 1 || throw(BoundsError())
    size(df, 1) == length(item) || size(df, 1) == 0 || error("number of rows does not match")

    insert!(index(df), col_ind, name)
    insert!(df.columns, col_ind, item)
    df
end

# FIXME: Needed to work around a crash: JuliaLang/julia#18299
function Base.insert!(df::DataFrame, col_ind::Int, item::NullableArray, name::Symbol)
    0 < col_ind <= ncol(df) + 1 || throw(BoundsError())
    size(df, 1) == length(item) || size(df, 1) == 0 || error("number of rows does not match")

    insert!(index(df), col_ind, name)
    insert!(df.columns, col_ind, item)
    df
end

function Base.insert!(df::DataFrame, col_ind::Int, item, name::Symbol)
    insert!(df, col_ind, upgrade_scalar(df, item), name)
end

function Base.merge!(df::DataFrame, others::AbstractDataFrame...)
    for other in others
        for n in _names(other)
            df[n] = other[n]
        end
    end
    return df
end

##############################################################################
##
## Copying
##
##############################################################################

# A copy of a DataFrame points to the original column vectors but
#   gets its own Index.
Base.copy(df::DataFrame) = DataFrame(copy(columns(df)), copy(index(df)))

# Deepcopy is recursive -- if a column is a vector of DataFrames, each of
#   those DataFrames is deepcopied.
function Base.deepcopy(df::DataFrame)
    DataFrame(deepcopy(columns(df)), deepcopy(index(df)))
end

##############################################################################
##
## Deletion / Subsetting
##
##############################################################################

# delete!() deletes columns; deleterows!() deletes rows
# delete!(df, 1)
# delete!(df, :Old)
function Base.delete!(df::DataFrame, inds::Vector{Int})
    for ind in sort(inds, rev = true)
        if 1 <= ind <= ncol(df)
            splice!(df.columns, ind)
            delete!(index(df), ind)
        else
            throw(ArgumentError("Can't delete a non-existent DataFrame column"))
        end
    end
    return df
end
Base.delete!(df::DataFrame, c::Int) = delete!(df, [c])
Base.delete!(df::DataFrame, c::Any) = delete!(df, index(df)[c])

# deleterows!()
function deleterows!(df::DataFrame, ind::@compat(Union{Integer, UnitRange{Int}}))
    for i in 1:ncol(df)
        df.columns[i] = deleteat!(df.columns[i], ind)
    end
    df
end

function deleterows!(df::DataFrame, ind::AbstractVector{Int})
    ind2 = sort(ind)
    n = size(df, 1)

    idf = 1
    iind = 1
    ikeep = 1
    keep = Array(Int, n-length(ind2))
    while idf <= n && iind <= length(ind2)
        1 <= ind2[iind] <= n || error(BoundsError())
        if idf == ind2[iind]
            iind += 1
        else
            keep[ikeep] = idf
            ikeep += 1
        end
        idf += 1
    end
    keep[ikeep:end] = idf:n

    for i in 1:ncol(df)
        df.columns[i] = df.columns[i][keep]
    end
    df
end

##############################################################################
##
## Hcat specialization
##
##############################################################################

# hcat! for 2 arguments
function hcat!(df1::DataFrame, df2::AbstractDataFrame)
    u = add_names(index(df1), index(df2))
    for i in 1:length(u)
        df1[u[i]] = df2[i]
    end

    return df1
end
hcat!(df::DataFrame, x::CategoricalArray) = hcat!(df, DataFrame(Any[x]))
hcat!(df::DataFrame, x::NullableCategoricalArray) = hcat!(df, DataFrame(Any[x]))
hcat!(df::DataFrame, x::NullableVector) = hcat!(df, DataFrame(Any[x]))
hcat!(df::DataFrame, x::Vector) = hcat!(df, DataFrame(Any[NullableArray(x)]))
hcat!(df::DataFrame, x) = hcat!(df, DataFrame(Any[NullableArray([x])]))

# hcat! for 1-n arguments
hcat!(df::DataFrame) = df
hcat!(a::DataFrame, b, c...) = hcat!(hcat!(a, b), c...)

# hcat
Base.hcat(df::DataFrame, x) = hcat!(copy(df), x)

##############################################################################
##
## Nullability
##
##############################################################################

function nullable!(df::DataFrame, col::ColumnIndex)
    df[col] = NullableArray(df[col])
    df
end
function nullable!{T <: ColumnIndex}(df::DataFrame, cols::Vector{T})
    for col in cols
        nullable!(df, col)
    end
    df
end

##############################################################################
##
## Pooling
##
##############################################################################

function categorical!(df::DataFrame, cname::@compat(Union{Integer, Symbol}), compact::Bool=true)
    df[cname] = categorical(df[cname], compact)
    return
end

function categorical!{T <: @compat(Union{Integer, Symbol})}(df::DataFrame, cnames::Vector{T},
                                                            compact::Bool=true)
    for cname in cnames
        df[cname] = categorical(df[cname], compact)
    end
    return
end

function categorical!(df::DataFrame, compact::Bool=true)
    for i in 1:size(df, 2)
        if eltype(df[i]) <: AbstractString
            df[i] = categorical(df[i], compact)
        end
    end
    return
end

function Base.append!(df1::DataFrame, df2::AbstractDataFrame)
   _names(df1) == _names(df2) || error("Column names do not match")
   eltypes(df1) == eltypes(df2) || error("Column eltypes do not match")
   ncols = size(df1, 2)
   # TODO: This needs to be a sort of transaction to be 100% safe
   for j in 1:ncols
       append!(df1[j], df2[j])
   end
   return df1
end

function Base.convert(::Type{DataFrame}, A::Matrix)
    n = size(A, 2)
    cols = Array(Any, n)
    for i in 1:n
        cols[i] = A[:, i]
    end
    return DataFrame(cols, Index(gennames(n)))
end

function _dataframe_from_associative(dnames, d::Associative)
    p = length(dnames)
    p == 0 && return DataFrame()
    columns  = Array(Any, p)
    colnames = Array(Symbol, p)
    n = length(d[dnames[1]])
    for j in 1:p
        name = dnames[j]
        col = d[name]
        if length(col) != n
            throw(ArgumentError("All columns in Dict must have the same length"))
        end
        columns[j] = NullableArray(col)
        colnames[j] = Symbol(name)
    end
    return DataFrame(columns, Index(colnames))
end

function Base.convert(::Type{DataFrame}, d::Associative)
    dnames = collect(keys(d))
    return _dataframe_from_associative(dnames, d)
end

# A Dict is not sorted or otherwise ordered, and it's nicer to return a
# DataFrame which is ordered in some way
function Base.convert(::Type{DataFrame}, d::Dict)
    dnames = collect(keys(d))
    sort!(dnames)
    return _dataframe_from_associative(dnames, d)
end

##############################################################################
##
## push! a row onto a DataFrame
##
##############################################################################

function Base.push!(df::DataFrame, associative::Associative{Symbol,Any})
    i = 1
    for nm in _names(df)
        try
            push!(df[nm], associative[nm])
        catch
            #clean up partial row
            for j in 1:(i - 1)
                pop!(df[_names(df)[j]])
            end
            msg = "Error adding value to column :$nm."
            throw(ArgumentError(msg))
        end
        i += 1
    end
end

function Base.push!(df::DataFrame, associative::Associative)
    i = 1
    for nm in _names(df)
        try
            val = get(() -> associative[string(nm)], associative, nm)
            push!(df[nm], val)
        catch
            #clean up partial row
            for j in 1:(i - 1)
                pop!(df[_names(df)[j]])
            end
            msg = "Error adding value to column :$nm."
            throw(ArgumentError(msg))
        end
        i += 1
    end
end

# array and tuple like collections
function Base.push!(df::DataFrame, iterable::Any)
    if length(iterable) != length(df.columns)
        msg = "Length of iterable does not match DataFrame column count."
        throw(ArgumentError(msg))
    end
    i = 1
    for t in iterable
        try
            push!(df.columns[i], t)
        catch
            #clean up partial row
            for j in 1:(i - 1)
                pop!(df.columns[j])
            end
            msg = "Error adding $t to column :$(_names(df)[i]). Possible type mis-match."
            throw(ArgumentError(msg))
        end
        i += 1
    end
end


text:
using JuMP

m = Model()

# ----- Variables ----- #
@variable(m, objvar)
x_Idx = Any[1, 2]
@variable(m, x[x_Idx])
set_lower_bound(x[1], -10.0)
set_upper_bound(x[1], 10.0)
set_lower_bound(x[2], -10.0)
set_upper_bound(x[2], 10.0)

# ----- Constraints ----- #
@NLconstraint(m, e1, -((-5+sqrt( (x[1])^2))*(-5+sqrt( (x[1])^2))+(-5+sqrt( (x[2])^2))*(-5+sqrt( (x[2])^2)))+objvar == 0.0)

# ----- Objective ----- #
@objective(m, Min, objvar)

m = m 		 # model get returned when including this script.

ext: jl

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: MatthewStuber/EAGO


text:
using Documenter, VortexLattice

makedocs(;
    modules = [VortexLattice],
    pages = [
        "Home" => "index.md",
        "Getting Started" => "guide.md",
        "Examples" => "examples.md",
        "Library" => "library.md",
        "Theory" => "theory.md"
    ],
    sitename = "VortexLattice.jl",
    authors = "Taylor McDonnell <taylormcd@byu.edu> and Andrew Ning <aning@byu.edu>",
)

deploydocs(
    repo = "github.com/byuflowlab/VortexLattice.jl.git",
)


text:
module GPLinearODEMaker

using LinearAlgebra
using SpecialFunctions
using Random
using Distributed
using SharedArrays
using SymEngine
using IterativeSolvers
using Statistics

# using Pkg
# Pkg.activate(".")
# Pkg.add("LinearAlgebra")
# Pkg.add("SpecialFunctions")
# Pkg.add("Random")
# Pkg.add("Distributed")
# Pkg.add("SharedArrays")
# Pkg.add("SymEngine")
# Pkg.pin(PackageSpec(name="SymEngine", version="0.6.0"))
# Pkg.add("IterativeSolvers")
# Pkg.add("Documenter")
# Pkg.add("Statistics")

import Base.ndims

include("general_functions.jl")
include("problem_definition_functions.jl")
include("gp_functions.jl")
include("prior_functions.jl")
include("kernel_base_functions.jl")
include("kernel_creation_functions.jl")
include("diagnostic_functions.jl")

end # module


text:
using Oceananigans.Grids: interior_parent_indices

struct PreconditionedConjugateGradientSolver{A, S}
    architecture :: A
        settings :: S
end

function PreconditionedConjugateGradientSolver(; arch = arch, grid = nothing, boundary_conditions = nothing, parameters = parameters)

    bcs = boundary_conditions
    if isnothing(boundary_conditions)
        bcs = parameters.Template_field.boundary_conditions
    end

    if isnothing(grid)
        grid = parameters.Template_field.grid
    end

    maxit = parameters.maxit
    if isnothing(parameters.maxit)
        maxit = grid.Nx * grid.Ny * grid.Nz
    end

    tol = parameters.tol
    if isnothing(parameters.tol)
        tol = 1e-13
    end

    tf = parameters.Template_field
    if isnothing(parameters.Template_field)
        tf = CenterField(arch, grid)
    end

    a_res = similar(tf.data)
    q     = similar(tf.data)
    p     = similar(tf.data)
    z     = similar(tf.data)
    r     = similar(tf.data)
    RHS   = similar(tf.data)
    x     = similar(tf.data)

    a_res.parent .= 0

    if isnothing(parameters.PCmatrix_function)
        # preconditioner not provided, use the Identity matrix
        PCmatrix_function(x; args...) = x
    else
        PCmatrix_function = parameters.PCmatrix_function
    end

    M(x; args...) = PCmatrix_function(x; args...)

    ii = interior_parent_indices(Center, topology(grid, 1), grid.Nx, grid.Hx)
    ji = interior_parent_indices(Center, topology(grid, 2), grid.Ny, grid.Hy)
    ki = interior_parent_indices(Center, topology(grid, 3), grid.Nz, grid.Hz)

    dotproduct(x,y)  = mapreduce((x,y)->x*y, + , view(x, ii, ji, ki), view(y, ii, ji, ki))
    norm(x)          = ( mapreduce((x)->x*x, + , view(x, ii, ji, ki)   ) )^0.5

    Amatrix_function = parameters.Amatrix_function
    A(x; args...) = ( Amatrix_function(a_res, x, arch, grid, bcs; args...); return  a_res )

    reference_pressure_solver = nothing
    if haskey(parameters, :reference_pressure_solver )
        reference_pressure_solver = parameters.reference_pressure_solver
    end

    settings = (
                                q=q,
                                p=p,
                                z=z,
                                r=r,
                                x=x,
                              RHS=RHS,
                              bcs=bcs,
                             grid=grid,
                                A=A,
                                M=M,
                            maxit=maxit,
                              tol=tol,
                          dotprod=dotproduct,
                             norm=norm,
                             arch=arch,
        reference_pressure_solver=reference_pressure_solver,
    )

   return PreconditionedConjugateGradientSolver(arch, settings)
end

@kernel function compute_residual!(r, RHS, A)
    i, j, k = @index(Global, NTuple)
    @inbounds r[i, j, k] = RHS[i, j, k] - A[i, j, k]
end

function quick_launch!(arch, grid, kernel!, args...)
    event = launch!(arch, grid, :xyz, kernel!, args...; dependencies=Event(device(arch)))
    wait(device(arch), event)
    return nothing
end

using Statistics

function solve_poisson_equation!(solver::PreconditionedConjugateGradientSolver, RHS, x; args...)
#
# Alg - see Fig 2.5 The Preconditioned Conjugate Gradient Method in
#                    "Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods"
#                    Barrett et. al, 2nd Edition.
#
#     given
#        linear Preconditioner operator M as a function M()
#        linear A matrix operator A as a function A()
#        a dot product function norm()
#        a right-hand side b
#        an initial guess x
#
#        local vectors: z, r, p, q
#
#     β  = 0
#     r .= b-A(x)
#     i  = 0
#
#     loop:
#      if i > MAXIT
#       break
#      end
#      z = M( r )
#      ρ    .= dotprod(r,z)
#      p = z+β*p
#      q = A(p)
#      α=ρ/dotprod(p,q)
#      x=x.+αp
#      r=r.-αq
#      if norm2(r) < tol
#       break
#      end
#      i=i+1
#      ρⁱᵐ1 .= ρ
#      β    .= ρⁱᵐ¹/ρ
#
    arch       = solver.architecture
    sset       = solver.settings
    grid       = sset.grid
    z, r, p, q = sset.z, sset.r, sset.p, sset.q
    A          = sset.A
    M          = sset.M
    maxit      = sset.maxit
    tol        = sset.tol
    dotprod    = sset.dotprod
    norm       = sset.norm

    β    = 0.0
    i    = 0
    ρ    = 0.0
    ρⁱᵐ¹ = 0.0

    r.parent .= 0
    # quick_launch!(arch, grid, compute_residual!, r, RHS, A(x))
    r.parent .= RHS.parent .- A(x; args...).parent
    ### println("PreconditionedConjugateGradientSolver ", i," RHS ", norm(RHS.parent) )
    ### println("PreconditionedConjugateGradientSolver ", i," A(x) ", norm(A(x; args...).parent) )

    while true
        ### println("PreconditionedConjugateGradientSolver ", i," norm(r.parent) ", norm(r.parent) )
        i >= maxit && break
        norm(r.parent) <= tol && break

        z.parent       .= M(r; args...).parent
        ### println("PreconditionedConjugateGradientSolver ", i," norm(z.parent) ", norm(z.parent) )
        ρ        = dotprod(z.parent, r.parent)
        ### println("PreconditionedConjugateGradientSolver ", i," ρ ", ρ )

        if i == 0
            p.parent   .= z.parent
        else
            β    = ρ/ρⁱᵐ¹
            p.parent   .= z.parent .+ β .* p.parent
        end

        q.parent       .= A(p; args...).parent
        ### println("PreconditionedConjugateGradientSolver ", i," norm(q.parent) ", norm(q.parent) )
        α        = ρ / dotprod(p.parent, q.parent)
        ### println("PreconditionedConjugateGradientSolver ", i," α ", α )
        x.parent       .= x.parent .+ α .* p.parent
        r.parent       .= r.parent .- α .* q.parent

        i     = i+1
        ρⁱᵐ¹  = ρ
    end

#==
#     No preconditioner verison
      i    = 0
      r   .= RHS .- A(x)
      p   .= r
      γ    = dotprod(r,r)
      while true
       if i > maxit
        break
       end
       q   .= A(p)
       α    = γ/dotprod(p,q)
       x   .= x .+ α .* p
       r   .= r .- α .* q
       println("Solver ", i," ", norm(r) )
       if norm(r) <= tol
        break
       end
       γ⁺   = dotprod(r,r)
       β    = γ⁺/γ
       p   .= r .+ β .* p
       γ    = γ⁺
       i    = i+1
      end
==#
      println("PreconditionedConjugateGradientSolver ", i," ", norm(r.parent) )

    fill_halo_regions!(x, sset.bcs, sset.arch, sset.grid)
    return x, norm(r.parent)
end

function Base.show(io::IO, solver::PreconditionedConjugateGradientSolver)
    print(io, "Oceanigans compatible preconditioned conjugate gradient solver.\n")
    print(io, " Problem size = "  , size(solver.settings.q) )
    print(io, "\n Boundary conditions = "  , solver.settings.bcs  )
    print(io, "\n Grid = "  , solver.settings.grid  )
    return nothing
end


text:
@symbol_func function phi_ind(cur_reactor::AbstractReactor)
  cur_phi = phi_RU(cur_reactor)

  cur_phi -= phi_res(cur_reactor)

  cur_phi
end


text:
################################################################################
abstract type BehaviorBelief end

struct DiscreteBehaviorBelief <: BehaviorBelief
    physical::MLPhysicalState
    models::AbstractVector
    weights::Vector{Vector{Float64}}
end
DiscreteBehaviorBelief(ps::MLPhysicalState, models::AbstractVector) = DiscreteBehaviorBelief(ps, models, Array(Vector{Float64}, length(ps.cars)))

function rand(rng::AbstractRNG,
              b::DiscreteBehaviorBelief,
              s::MLState=MLState(b.physical, Array(CarState, length(b.physical.cars))))
    s.crashed = b.physical.crashed
    resize!(s.cars, length(b.physical.cars))
    for i in 1:length(s.cars)
        s.cars[i] = CarState(b.physical.cars[i], sample(rng, b.models, WeightVec(b.weights[i])))
    end
    return s
end

@with_kw type WeightUpdateParams
    smoothing::Float64 = 0.02 # value between 0 and 1, adds this fraction of the max to each entry in the vector
    wrong_lane_factor::Float64 = 0.1
end

function weights_from_particles!(b::DiscreteBehaviorBelief, problem::NoCrashProblem, o::MLPhysicalState, particles, p::WeightUpdateParams)
    b.physical = o
    resize!(b.weights, length(o.cars))
    for i in 1:length(o.cars)
        b.weights[i] = zeros(length(b.models)) #XXX lots of allocation
    end
    for sp in particles
        isp = 1
        io = 1
        while io <= length(o.cars) && isp <= length(sp.cars)
            co = o.cars[io]
            csp = sp.cars[isp]
            if co.id == csp.id
                # sigma_acc = dmodel.vel_sigma/dt
                # dv = acc*dt
                # sigma_v = sigma_dv = dmodel.vel_sigma
                if abs(co.x-csp.x) < 0.3*problem.dmodel.phys_param.lane_length
                    proportional_likelihood = proportional_normal_pdf(csp.vel,
                                                                      co.vel,
                                                                      problem.dmodel.vel_sigma)
                    if co.y == csp.y
                        b.weights[io][csp.behavior.idx] += proportional_likelihood
                    elseif abs(co.y - csp.y) < 1.0
                        b.weights[io][csp.behavior.idx] += p.wrong_lane_factor*proportional_likelihood
                    end # if greater than one lane apart, do nothing
                end
                io += 1
                isp += 1
            elseif co.id < csp.id
                io += 1
            else
                @assert co.id > csp.id
                isp += 1
            end
        end
    end

    # smoothing
    for i in 1:length(b.weights)
        if sum(b.weights[i]) > 0.0
            b.weights[i] .+= p.smoothing * sum(b.weights[i])
        else
            b.weights[i] .+= 1.0
        end
    end

    return b
end

#abstract type Updater_DiscreteBehaviorBelief <: Updater{DiscreteBehaviorBelief} end
type BehaviorBeliefUpdater <: Updater  #{DiscreteBehaviorBelief}
    problem::NoCrashProblem
    smoothing::Float64
    nb_particles::Float64
end
function set_problem!(u::BehaviorBeliefUpdater, problem::Union{POMDP,MDP})
    u.problem = problem
end
#=
function update(up::BehaviorBeliefUpdater,
                b_old::DiscreteBehaviorBelief,
                a::MLAction,
                o::MLPhysicalState,
                b_new::DiscreteBehaviorBelief=DiscreteBehaviorBelief(o,
                                                       up.problem.dmodel.behaviors,
                                                       Array(Vector{Float64}, 0)))
    # particles = SharedArray(MLState, )
    # @parallel
    # # run simulationsabstract type

    # weights_from_particles!(b_new, )
    error("not implemented")
end
=#

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

#=
function update_cz(up::BehaviorBeliefUpdater,
                b_old::DiscreteBehaviorBelief,
                a::MLAction,
                o::MLPhysicalState,
                b_new::DiscreteBehaviorBelief=DiscreteBehaviorBelief(o,
                                                       up.problem.dmodel.behaviors,
                                                       Array(Vector{Float64}, 0)))

end
=#

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

type BehaviorRootUpdater <: Updater #{BeliefNode}
    problem::NoCrashProblem
    params::WeightUpdateParams
    min_particles::Int
    rng::AbstractRNG
end
BehaviorRootUpdater(p::NoCrashProblem, w::WeightUpdateParams, rng::AbstractRNG=Base.GLOBAL_RNG) = BehaviorRootUpdater(p,w,0,rng)
function set_problem!(u::BehaviorRootUpdater, problem::Union{POMDP,MDP})
    u.problem = problem
end

#initialize_belief(up::BehaviorRootUpdater, b) = BasicPOMCP.RootNode(b)
#=
function update(up::BehaviorRootUpdater,
                b_old::BeliefNode,
                a::MLAction,
                o::MLPhysicalState,
                b_new::RootNode=RootNode(DiscreteBehaviorBelief(o, up.problem.dmodel.behaviors.models)))

    b_new.B.models = up.problem.dmodel.behaviors.models   #### space for all behaviors
    #XXX hack
    particles = Iterators.chain([child.B.particles for child in values(b_old.children[a].children)]...)

    if up.min_particles > 0
        particles = collect(particles)
        sizehint!(particles, up.min_particles)
        for i in 1:up.min_particles-length(particles)
            s = rand(up.rng, b_old)
            push!(particles, generate_s(up.problem, s, a, up.rng))
        end
    end

    weights_from_particles!(b_new.B, up.problem, o, particles, up.params)
    return b_new
end
=#
proportional_normal_pdf(x, mu, sigma) = exp(-(x-mu)^2/(2.*sigma^2))
#=
function count_particles!(d::Dict{Int,Int}, n::POMCP.BeliefNode)
    for anode in values(n.children)
        for child in values(anode.children)
            count_particles!(d, child)
        end
    end
    if isa(n.B, POMCP.ParticleCollection)
        nb_particles = length(n.B.particles)
        if haskey(d, nb_particles)
            d[nb_particles] += 1
        else
            d[nb_particles] = 1
        end
    end
end
=#

# TODO clean up
#=
import DESPOT: fringe_upper_bound, lower_bound, init_lower_bound

function fringe_upper_bound(pomdp::NoCrashPOMDP, s::MLState)
  """
  Assume there are no cars in the environment: just go straight to target lane, and stay there
    for the rest of the time horizon
  """
  # NOTE: This function makes assumptions about the length of the time horizon!!!!
  T = 100 #time horizon
  pp = pomdp.dmodel.phys_param
  nb_lanes = pp.nb_lanes
  current_lane = s.cars[1].y
  t_single_lane_change = 1.0 / (pp.dt * pomdp.dmodel.lane_change_rate) #nb timesteps
  t_target_lane = (nb_lanes - current_lane) * t_single_lane_change
  return sum([pomdp.rmodel.reward_in_target_lane * discount(pomdp)^t for t = 0:(T-t_target_lane-1)])
end

type NoCrashLB <: DESPOTLowerBound end

init_lower_bound(lb::NoCrashLB, pomdp::NoCrashPOMDP, config::DESPOTConfig) = begin end

function lower_bound(lb::NoCrashLB, pomdp::NoCrashPOMDP, particles::Vector{DESPOTParticle{MLState}}, config:DESPOTConfig)
  """
  Just rollout the heuristic policy--it only really depends on observations, so
    bmodel (eg true behavior states) don't matter
  """

  s = particles[1].state

  sim = RolloutSimulator(initial_state=s, max_steps=100)

  policy = Simple(pomdp) # TODO check if this signature works with NoCrashPOMDP and stuff

  return simulate(sim, pomdp, policy)
end
=#

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
# #Inflación media ponderada
# InflationWeightedMean.jl - Función de inflación media ponderada interanual

"""
    InflationWeightedMean <: InflationFunction
Función de inflación para computar la media ponderada.
"""
struct InflationWeightedMean <: InflationFunction
end

# 2. Extender el método de nombre
measure_name(::InflationWeightedMean) = "Media ponderada interanual"

# Define cómo opera InflationSimpleMean sobre un objeto de tipo VarCPIBase.
function (inflfn::InflationWeightedMean)(base::VarCPIBase{T}) where T
    #Obtener el indice correspondiente a las variaciones intermensuales
    indmat = capitalize(base.v)
    #Cálculo de la variación interanual
    y2ymat = varinteran(indmat)
    #Cálculo de la media ponderada de las variaciones interanuales
    y2yweightedmean = y2ymat*base.w/100
    #Utilizar el índice de la media ponderada de las variaciones intermensuales de los primeros
    #11 meses para completar la serie del índice de media ponderada de variaciones interanuales
    m2mweightedmean = base.v*base.w/100
    weightedmeanind = zeros(T, periods(base)+1); weightedmeanind[1]=100;
    weightedmeanind[2:12] = capitalize(m2mweightedmean[1:11])
    for j = 13:periods(base)+1
        weightedmeanind[j] = (y2yweightedmean[j-12]/100 + 1)*weightedmeanind[j-12]
    end
    #Cálculo de la variación intermensual del índice de media simple
    weightedmean = varinterm(weightedmeanind[2:end])
    weightedmean
end


text:
import JuliaInterpreter.eval_code

# Simple evaling of function argument
function evalfoo1(x,y)
    x+y
end
frame = JuliaInterpreter.enter_call(evalfoo1, 1, 2)
@test eval_code(frame, "x") == 1
@test eval_code(frame, "y") == 2

# Evaling with sparams
evalsparams(x::T) where T = x
frame = JuliaInterpreter.enter_call(evalsparams, 1)
@test eval_code(frame, "x") == 1
eval_code(frame, "x = 3")
@test eval_code(frame, "x") == 3
@test eval_code(frame, "T") == Int
eval_code(frame, "T = Float32")
@test eval_code(frame, "T") == Float32

# Evaling with keywords
evalkw(x; bar=true) = x
frame = JuliaInterpreter.enter_call(evalkw, 2)
frame = JuliaInterpreter.maybe_step_through_wrapper!(frame)
@test eval_code(frame, "x") == 2
@test eval_code(frame, "bar") == true
eval_code(frame, "bar = false")
@test eval_code(frame, "bar") == false

# Evaling with symbols
evalsym() = (x = :foo)
frame = JuliaInterpreter.enter_call(evalsym)
# Step until the local actually end up getting defined
JuliaInterpreter.step_expr!(frame)
JuliaInterpreter.step_expr!(frame)
@test eval_code(frame, "x") == :foo

# Evaling multiple statements (https://github.com/JuliaDebug/Debugger.jl/issues/188)
frame = JuliaInterpreter.enter_call(evalfoo1, 1, 2)
@test eval_code(frame, "x = 1; y = 2") == 2
@test eval_code(frame, "x") == 1
@test eval_code(frame, "y") == 2

# https://github.com/JuliaDebug/Debugger.jl/issues/177
function f()
    x = 1
    f = ()->(x = 2)
    f()
    x
end
frame = JuliaInterpreter.enter_call(f)
JuliaInterpreter.step_expr!(frame)
JuliaInterpreter.step_expr!(frame)
@test eval_code(frame, "x") == 1
eval_code(frame, "x = 3")
@test eval_code(frame, "x") == 3
JuliaInterpreter.finish!(frame)
@test JuliaInterpreter.get_return(frame) == 2


text:
@testset "NewtonCorrector" begin
    @testset "square" begin
        F = equations(griewank_osborne())

        prob, starts = problem_startsolutions(F; seed = 12345)
        S = collect(starts)

        x₀ = S[1] .+ 1e-4
        t = 1.0
        H = HC.HomotopyWithCache(prob.homotopy, x₀, t)
        corr_cache = HC.cache(HC.NewtonCorrector(), H, x₀, t)

        JM = HC.JacobianMonitor(jacobian(H, x₀, t))
        x̄ = similar(x₀)

        R1 = HC.newton!(x̄, H, x₀, t, JM, HC.InfNorm(), corr_cache; tol = 1e-3)
        @test R1.iters == 1
        @test R1.accuracy ≤ 1e-3
        @test R1.norm_Δx₀ < 1e-1
        @test isnan(R1.ω)
        @test norm(x̄ - S[1]) < 1e-6

        R2 = HC.newton!(x̄, H, x₀, t, JM, HC.InfNorm(), corr_cache; tol = 1e-7)
        @test R2.return_code == HC.NEWT_CONVERGED
        @test R2.iters == 2
        @test R2.accuracy < 1e-7
        @test !isnan(R2.ω + R2.ω₀ + R2.θ + R2.θ₀)
        @test R2.norm_Δx₀ < 1.1e-4
        @test norm(x̄ - S[1]) < 1e-10

        R3 = HC.newton!(x̄, H, [1e-5, 1e-5], 0.0, JM, HC.InfNorm(), corr_cache)
        @test R3.return_code == HC.NEWT_TERMINATED
    end

    @testset "overdetermined" begin
        @polyvar x y
        F = SPSystem([x^2 + y^2 - 1, x - y, (x - y) * (x^2 + y^2 - 1)])
        t = rand()
        s = [√(0.5), √(0.5)]
        x₀ = s .+ 1e-4 .* rand(ComplexF64, 2)
        H = HomotopyWithCache(StraightLineHomotopy(F, F), x₀, t)

        corr_cache = HC.cache(HC.NewtonCorrector(), H, x₀, t)
        JM = HC.JacobianMonitor(jacobian(H, x₀, t))
        x̄ = similar(x₀)

        R1 = HC.newton!(x̄, H, x₀, t, JM, HC.InfNorm(), corr_cache; tol = 1e-3)
        @test R1.iters == 1
        @test R1.accuracy ≤ 1e-3
        @test R1.norm_Δx₀ ≤ 1e-3

        R2 = HC.newton!(x̄, H, x₀, t, JM, HC.InfNorm(), corr_cache; tol = 1e-7)
        @test R2.return_code == HC.NEWT_CONVERGED
        @test R2.iters == 2
        @test R2.accuracy < 1e-7
        @test !isnan(R2.ω + R2.ω₀ + R2.θ + R2.θ₀)
        @test R2.norm_Δx₀ < 1e-3
        @test norm(x̄ - s) < 1e-10
    end
end


text:
using KissSmoothing
using Test
using Random

@testset "Goodness of Smoothing" begin
    rng = Random.MersenneTwister(1337)
    for N in [100, 200, 400, 800]
        for α in [0.01, 0.05, 0.1]
            x = LinRange(0, 1, N)
            y = identity.(x)
            n = randn(rng, length(x)) .* α
            yr = y .+ n
            ys, yn = denoise(yr)
            c1 = sum(abs2, yr .- y)
            c2 = sum(abs2, ys .- y)
            @test c2 < c1
            s1 = sum(abs2, n)
            s2 = sum(abs2, yn)
            @test abs(log2(s2 / s1)) < 1
        end
    end
end

@testset "Too few points" begin
    rng = Random.MersenneTwister(1337)
    O = [1.0, 2.0, 10.0]
    S, N = denoise(O)
    @test all(O .== S)
    @test all(N .== 0)
    @test length(S) == 3
    @test length(N) == 3
end

@testset "Infinite smoothing" begin
    O = sign.(sin.(1:1000)) .+ 1
    S, N = denoise(O, factor = Inf)
    @test all(S .≈ 1)
    @test all(abs.(N) .≈ 1)
end

@testset "Verbose" begin
    O = sign.(sin.(1:1000)) .+ 1
    S, N = denoise(O, factor = 0.0, verbose = true)
    @test all(abs.(S .- O) .< 1e-10)
end

@testset "Fit 1D RBF" begin
    for μ in LinRange(-100,100,5)
        t = LinRange(0,2pi,150)
        y = sin.(t).+μ
        fn = fit_rbf(t,y,LinRange(0,2pi,50))
        pred_y = fn(t)
        error = sqrt(sum(abs2, pred_y .- y)/length(t))
        @test error < 0.0002
    end
end


text:
using Gadfly, Distributions

set_default_plot_size(6inch, 3inch)

plot(x=rand(Poisson(20), 1000), Scale.x_discrete, Geom.histogram)


text:
function initial_conditions_numerical(slob::SLOB, p_0, V₀)
    ud = (-V₀/(2.0*slob.Δx) + slob.D/(slob.Δx^2)) * ones(Float64, slob.M)
    md = ((-2.0*slob.D)/(slob.Δx^2) - slob.nu) * ones(Float64, slob.M+1)
    ld = (V₀/(2.0*slob.Δx) + slob.D/(slob.Δx^2)) * ones(Float64, slob.M)
    A = Tridiagonal(ld, md, ud)

    A[1,2] = 2*slob.D/(slob.Δx^2)
    A[end, end-1] = 2*slob.D/(slob.Δx^2)

    B = .-[slob.source_term(xᵢ, p_0) for xᵢ in slob.x]
    φ = A \ B
    return φ
end

function initial_conditions_numerical(slob::SLOB, p_0)
    ϵ = rand(Normal(0.0, 1.0))
    V₀ = sign(ϵ) * min(abs(slob.σ * ϵ), slob.Δx / slob.Δt)
    return initial_conditions_numerical(slob, p_0, V₀)
end

function extract_mid_price(slob, lob_density)
    mid_price_ind = 2
    while (lob_density[mid_price_ind] > 0) | (lob_density[mid_price_ind+1]>lob_density[mid_price_ind])
        mid_price_ind += 1
    end
    y1 = lob_density[mid_price_ind-1]
    y2 = lob_density[mid_price_ind]
    x1 = slob.x[mid_price_ind-1]

    mid_price = (-y1 * slob.Δx)/(y2 - y1) + x1
    return mid_price
end

function calculate_left_jump_probability(Z)
    return (exp(-Z))/(exp(Z) + exp(-Z) + 1.0)
end

function calculate_right_jump_probability(Z)
    return (exp(Z))/(exp(Z) + exp(-Z) + 1.0)
end

function calculate_self_jump_probability(Z)
    return (1.0)/(exp(Z) + exp(-Z) + 1.0)
end

function calculate_jump_probabilities(slob, V_t)
    Z = (3/4) * (V_t * slob.Δx) / (slob.D)
    # Z = (V_t * slob.Δx) / (slob.D)
    p⁻ = calculate_left_jump_probability(Z)
    p⁺ = calculate_right_jump_probability(Z)
    p = calculate_self_jump_probability(Z)
    return p⁺, p⁻, p
end

function get_sub_period_time(slob, t, time_steps)
    if slob.α <= 0.0
        return 0.0, time_steps - t + 1
    end
    τ = rand(Exponential(slob.α))
    remaining_time = time_steps - t + 1
    τ_periods = min(floor(Int, τ/slob.Δt), remaining_time)
    return τ, τ_periods
end

function intra_time_period_simulate(slob, φ, p)
    ϵ = rand(Normal(0.0, 1.0))
    V_t = sign(ϵ) * min(abs(slob.σ * ϵ), slob.Δx / slob.Δt)

    P⁺, P⁻, P = calculate_jump_probabilities(slob, V_t)

    φ₋₁ = φ[1]
    φ₊₁ = φ[end]
    φ_next = zeros(Float64, size(φ,1))

    φ_next[1] = P⁺ * φ₋₁ + P⁻ * φ[2] + P * φ[1] -
        slob.nu * slob.Δt * φ[1] + slob.Δt * slob.source_term(slob.x[1], p)

    φ_next[end] = P⁻ * φ₊₁ + P⁺ * φ[end-1] + P * φ[end] -
        slob.nu * slob.Δt * φ[end] + slob.Δt * slob.source_term(slob.x[end], p)

    φ_next[2:end-1] = P⁺ * φ[1:end-2] + P⁻ * φ[3:end] + P * φ[2:end-1] -
        slob.nu * slob.Δt * φ[2:end-1] +
        [slob.Δt * slob.source_term(xᵢ, p) for xᵢ in slob.x[2:end-1]]

    return φ_next, P⁺, P⁻, P
end

function dtrw_solver(slob::SLOB)
    time_steps = get_time_steps(slob.T, slob.Δt)
    p = ones(Float64, time_steps + 1)
    p[1] = slob.p₀

    t = 1
    φ = initial_conditions_numerical(slob, p[t], 0.0)

    while t <= time_steps
        τ, τ_periods = get_sub_period_time(slob, t, time_steps)

        for τₖ = 1:τ_periods
            t += 1
            φ, _, _, _  = intra_time_period_simulate(slob, φ, p[t-1])
            p[t] = extract_mid_price(slob, φ)

        end
        if t > time_steps
            mid_prices = sample_mid_price_path(slob, p)
            return mid_prices
        end
        t += 1
        if slob.α > 0.0
            φ = initial_conditions_numerical(slob, p[t-1])
        end
        p[t] = extract_mid_price(slob, φ)
    end
    mid_prices = sample_mid_price_path(slob, p)
    return mid_prices
end

function dtrw_solver(debug::Bool, slob::SLOB)
    time_steps = get_time_steps(slob.T, slob.Δt)
    φ = ones(Float64, slob.M + 1, time_steps + 1)

    p = ones(Float64, time_steps + 1)
    mid_prices = ones(Float64, slob.T + 1)

    p[1] = slob.p₀
    mid_prices[1] = slob.p₀

    P⁺s = fill(1/3, time_steps)
    P⁻s = fill(1/3, time_steps)
    Ps = fill(1/3, time_steps)
    t = 1
    φ[:, t] = initial_conditions_numerical(slob, p[t], 0.0)

    while t <= time_steps
        τ, τ_periods = get_sub_period_time(slob, t, time_steps)

        for τₖ = 1:τ_periods
            t += 1
            φ[:, t], P⁺s[t-1], P⁻s[t-1], Ps[t-1]  = intra_time_period_simulate(slob,
                φ[:, t-1], p[t-1])
            try
                p[t] = extract_mid_price(slob, φ[:, t])
            catch e
                println("Bounds Error at t=$t")
                mid_prices = sample_mid_price_path(slob, p)
                return φ, p, mid_prices, P⁺s, P⁻s, Ps
            end

            @info "Intra-period simulation. tick price = R$(p[t]) @t=$t"
        end
        if t > time_steps
            mid_prices = sample_mid_price_path(slob, p)
            return φ, p, mid_prices, P⁺s, P⁻s, Ps
        end
        t += 1
        if slob.α > 0.0
            φ[:, t] = initial_conditions_numerical(slob, p[t-1])
        end

        try
            p[t] = extract_mid_price(slob, φ[:, t])
        catch e
            println("Initial Conditions Bounds Error at t=$t")
            mid_prices = sample_mid_price_path(slob, p)
            return φ, p, mid_prices, P⁺s, P⁻s, Ps
        end

        @info "LOB Density recalculated. tick price = R$(p[t]) @t=$t"
    end

    mid_prices = sample_mid_price_path(slob, p)
    return φ, p, mid_prices, P⁺s, P⁻s, Ps
end


text:
# this file is not used currently, left here just in case

function select(sockets::Vector{TCPSocket}, readfn::Function;
                sleep_time::Float64=0.1, channel_size::Int64=1024)
    result_channel = Channel(channel_size)
    for s in sockets
        @async begin
            Base.start_reading(s)
            while isopen(s)
                while nb_available(s) == 0 sleep(sleep_time) end
                result = readfn(s)
                put!(result_channel, result)
            end
        end
    end
    return result_channel
end

# test code for select()

function echo_server()
    server = listen(8080)
    while true
        conn = accept(server)
        @async begin
            try
                while true
                    line = readline(conn)
                    write(conn,line)
                end
            catch err
                print("connection ended with error $err")
            end
        end
    end
end

function example()
    @async echo_server()
    s1 = connect("localhost", 8080)
    s2 = connect("localhost", 8080)
    ch = select([s1, s2], sock -> readline(sock))
    write(s1, "hello\n")
    take!(ch)

end


text:
"""
    ransac(pc, params, setenabled; reset_rand = false)

Run the RANSAC algorithm on a pointcloud with the given parameters.

Return the extracted primitives and the time it took to run the algorithm (in seconds).

# Arguments
- `pc::RANSACCloud`: the point cloud.
- `params::NamedTuple`: parameters.
- `setenabled::Bool`: if `true`: set every point to enabled.
- `reset_rand::Bool=false`: if `true`, resets the random seed with `Random.seed!(1234)`
"""
function ransac(pc, params, setenabled; reset_rand = false)
    if setenabled
        for i in eachindex(pc.isenabled)
            pc.isenabled[i] = true
        end
    end
    ransac(pc, params, reset_rand=reset_rand)
end

"""
    ransac(pc, params; reset_rand = false)

Run the RANSAC algorithm on a pointcloud with the given parameters.

Return the extracted primitives and the time it took to run the algorithm (in seconds).

# Arguments
- `pc::RANSACCloud`: the point cloud.
- `params::NamedTuple`: parameters.
- `reset_rand::Bool=false`: if `true`, resets the random seed with `Random.seed!(1234)`
"""
function ransac(pc, params; reset_rand = false)
    reset_rand && Random.seed!(1234)

    @extract params : iter_params=iteration
    @extract iter_params : drawN minsubsetN prob_det τ
    @extract iter_params : itermax shape_types
    @extract iter_params : extract_s terminate_s

    # @unpack drawN, minsubsetN, prob_det, τ = params
    # @unpack itermax, shape_types = params
    # @unpack extract_s, terminate_s = params
    start_time = time_ns()

    # FittedShape[] is abstract array, which is bad
    candidates = FittedShape[]
    scoredshapes = IterationCandidates()
    extracted = ExtractedShape[]

    # save octree level of a FittedShape
    shape_octree_level = Int[]

    # smallest distance in the pointcloud
    #lsd = smallestdistance(pc.vertices)

    # allocate for the random selected points
    sd = Vector{Int}(undef, drawN)
    @logmsg IterInf "Iteration begins."

    # for logging
    notifit = itermax > 10 ? div(itermax,10) : 1

    # track the number of candidates for the probabilities
    # 1. :lengthC - number of candidates in a given iteration
    # 2. :allcand - number of candidates that have been ever scored
    # 3. :nofminset - number of minimal sets that have been drawn so far
    countcandidates = [0, 0, 0]

    # iterate begin
    for k in 1:itermax
        if count(pc.isenabled) < τ
            @logmsg IterInf "Break at $k iteration, because left less points, then τ"
            break
        end
        # generate minsubsetN candidate
        for i in 1:minsubsetN
            sampling_res = samplepointcloud4!(sd, pc, params)
            sampling_res[1] || continue

            # sd: indexes of the actually selected points
            #TODO: this should be something more general
            # fit plane to the selected points
            f_v = @view pc.vertices[sd]
            f_n = @view pc.normals[sd]

            forcefitshapes!(f_v, f_n, params, candidates, shape_octree_level, sampling_res[2], pc)
        end # for t

        # evaluate the compatible points, currently used as score
        # TODO: do something with octree levels and scores
        which_ = 1

        # update number of candidates
        countcandidates[2] += size(candidates,1)
        scorecandidates!(pc, scoredshapes, candidates, which_, params, shape_octree_level)
        # by now every candidate is scored into scoredshapes

        # considered so far k*minsubsetN pcs. minimal sets
        countcandidates[3] = k*minsubsetN
        # length of the candidate struct

        countcandidates[1] = length(scoredshapes)
        if ! (length(scoredshapes) < 1)
            best = findhighestscore(scoredshapes)
            bestshape = scoredshapes.shapes[best.index]
            # TODO: refine if best.overlap
            scr = E(scoredshapes.scores[best.index])
            #best_length = size(bestshape.inpoints, 1)

            s = chooseS(countcandidates, extract_s)
            ppp = prob(scr, s, pc.size, drawN)

            #info printing: there's a best
            if k%notifit == 0
                @logmsg IterInf "$k. it, best: $scr score, prob: $ppp, scored shapes: $(length(scoredshapes)) pcs."
            end

            # if the probability is large enough, extract the shape
            if ppp > prob_det

                # TODO: proper refit, not only getting the points that fit to that shape
                # what do you mean by refit?
                # refit on the whole pointcloud
                extr_shape = refit(bestshape, pc, params)

                if ! (extr_shape === nothing)
                    @logmsg IterInf "Extracting best: $(strt(bestshape)) score: $scr, size: $(size(extr_shape.inpoints,1)) ps."
                    # invalidate indexes
                    invalidate_indexes!(pc, extr_shape.inpoints)
                    # extract the shape and delete from scoredshapes
                    push!(extracted, extr_shape)
                    deleteat!(scoredshapes, best.index)
                    # remove invalid shapes
                    removeinvalidshapes!(pc, scoredshapes)
                end
            end # if extract shape
        else
            #info printing: not best
            if k%notifit == 0
                @logmsg IterInf "$k. it, no candidates."
            end
        end
        # update octree levels
        updatelevelweight(pc)

        # check exit condition
        s = chooseS(countcandidates, terminate_s)
        #if prob(τ/subsetN, s, pc.size, drawN) > prob_det
        if prob(τ, s, pc.size, drawN) > prob_det
            @logmsg IterInf "Break, at this point all shapes should be extracted: $k. iteration."
            break
        end
    end # iterate end
    fint = trunc((time_ns() - start_time)/1_000_000_000, digits=2)
    @logmsg IterInf "Iteration finished in $fint seconds with $(length(extracted)) extracted and $(length(scoredshapes)) scored shapes."
    return extracted, fint
end # ransac function


text:
"""
    mpifit_hsvgp(get_data, n_parts, n_dims, bounds_low, bounds_high;
    ni_low     = 5,
    ni_high    = 50,
    batch_size = 20,
    n_iters    = 200,
    n_subcycle_ELBO_local  = 20,
    n_subcycle_ELBO_global = 20)

Fit HSVGP to data for distributed computation

    get_data:                function to define how a local partition should read in its data
    ni_low:                  number of local inducing points
    ni_high:                 number of global inducing points
    batch_size:              size of batches for stochastic optimizer
    n_iters:                 number of optimization steps
    n_subcycle_ELBO_local:   subcycles for local  ELBO before MPI communication
    n_subcycle_ELBO_global:  subcycles for global ELBO before MPI communication

WARNING: Currently writes out several jld2 files containing local and global GPss
"""
function mpifit_hsvgp(get_data, n_parts, n_dims, bounds_low, bounds_high;
    ni_low     = 5,               # Number of input dimensions
    ni_high    = 50,              # Number of local inducing points
    batch_size = 20,              # Number of points used in updating local ELBO
    n_iters    = 200,             # Total updates in whole cycle
    n_subcycle_ELBO_local  = 20,  # subcycles for local ELBO before comm.
    n_subcycle_ELBO_global = 20)  # subcycles for global ELBO before comm.

    #TODO:
    #     Need to clean up saving out of local and global objects at the end.
    #     It's messy now and will get worse with many MPI processes.
    #
    #     Would be nice to write a function for prediction to identify which
    #     partition a new value of X would be in, then direct the prediction to
    #     that local svgp
    #
    #     Need to generalize the interface here. Should have a way to pass in
    #     the structure of data to be worked with. Maybe pass in a get_data
    #     function that returns the x,y data for each partition.
    #
    #     Would like to clean up some of those arguments also. I dont think
    #     we should need to have to pass all that stuff in, but not sure the
    #     best way to give global bounds information without simply specifying
    #     it for cases where all the data cannot be stored on one node.
    #

    MPI.Init()
    comm = MPI.COMM_WORLD
    myrank = MPI.Comm_rank(comm)
    mysize = MPI.Comm_size(comm)

    root = 0 # global node

    part_split = [[ii] for ii in 1:(mysize-1)]

    if n_parts > mysize - 1
        for ii = mysize:n_parts
            append!(part_split[1 + ((ii-1) % (mysize - 1))], ii)
        end
    end

    ######################## DATA SPLIT ########################################
    # This part need not be done separately. Data will already be stored on nodes

    # Root process is the global node and rest of the processes are local processes
    # datarow goes from 2 to 48601
    # numrowscsv = 48600
    # numrowseachrank = Int(numrowscsv/(mysize-1))

    if myrank != root

        data_array = [get_data(part) for part in part_split[myrank]]
        len        = length(part_split[myrank])

    end # rank != root

    # Initializing "trace" array to store ELBO values
    if myrank == root
        trace_h       = zeros(n_iters)          # To collect the global ELBO value
        store_trace_l = zeros(n_iters,n_parts) # Store trace_l from all proesses
    end
    if myrank != root
        trace_l = [zeros(n_iters) for ii in 1:len] # To collect the local ELBO value
    end

    # Now that each local node has seen the data, it will send xi_low, cov_low,
    # mean_low to global node.
    # Once the global node gets data from local nodes, it will evaluate xi_high,
    # cov_high, mean_high and theta_high
    # Following that, global node will calculate prior_mean and prior_cov and
    # send it back to local nodes.

    ###################### INITIALIZATION ######################################

    if myrank != root
        # Step 1
        # Build local SVGPs ###################
        local_svgps = [SVGP_obj(data_array[ii][1], float(data_array[ii][2]), ni_low) for ii in 1:len]

        for ii in 1:len
            # Step 2
            # Sending mean_factor, xi_low, cov_low, mean_low and N to global #######
            send_inducing_points_to_global(local_svgps[ii].params.inducing_locs, local_svgps[ii].params.inducing_mean, root, comm)
        end

        # Step 3
        # Initialize local optimizers
        lopt_lrho = [ADAM(1. * 0.0075) for ii in 1:len]     #     opt  = ADAM(0.0075, (0.5, 0.8))
        lopt_lkap = [ADAM(1. * 0.0075) for ii in 1:len]      #     opt  = ADAM(0.0075, (0.5, 0.8))
        lopt_lsig = [ADAM(1. * 0.0075) for ii in 1:len]      #     opt  = ADAM(0.0075, (0.5, 0.8))
        loptm     = [ADAM(1. * 0.01)   for ii in 1:len]        #     optm = ADAM(0.02, (0.5, 0.8))
        loptS     = [ADAM(1. * 0.002)  for ii in 1:len]       #     optS = ADAM(0.002, (0.5, 0.8))
        loptx     = [ADAM(1. * 0.01)   for ii in 1:len]        #     optx = ADAM(0.01, (0.5, 0.8))

    end # myrank != root

    if myrank == root

        # Step 3
            # Initialize global svgp and optimizers
            # Root process will receive inducing point data
            # Root process will update global parameters
            # Root process will calculate prior_mean prior_cov for local nodes
            # Root process will send prior_mean prior_cov back to local processes

        # Initialize global optimizers
        opt_lrho = ADAM(1. * 0.0075)     #     opt  = ADAM(0.0075, (0.5, 0.8))
        opt_lkap = ADAM(1. * 0.0075)     #     opt  = ADAM(0.0075, (0.5, 0.8))
        opt_lsig = ADAM(1. * 0.0075)     #     opt  = ADAM(0.0075, (0.5, 0.8))
        optm     = ADAM(1. * 0.01)       #     optm = ADAM(0.02, (0.5, 0.8))
        optS     = ADAM(1. * 0.002)      #     optS = ADAM(0.002, (0.5, 0.8))
        optx     = ADAM(1. * 0.01)       #     optx = ADAM(0.01, (0.5, 0.8))

        # Initialize global svgp
        # TODO:
        #     Need bounds on input space and some idea of scale for output
        #         to set initial GP variance parameter and
        #         Latin Hypercube for inducing point locations
        #
        #     Assuming we will work in the case that we can't send all inducing point
        #         data to global. Instead should probably get bounds as arguments to original
        #         function or otherwise.
        fake_data_to_init = SVGP_data(
            vcat(bounds_low', bounds_high'),
            randn(4)*10 # sampling 4 points with something like the variance I was to initialize at
        )
        global_svgp = SVGP_obj(ni_high, SVGP_data("global"), SVGP_params(fake_data_to_init, ni_high), "gaussian")

        for part = (root+1):n_parts
            recv_id = 1 + (part - 1) % (mysize - 1)
            # TODO: Would prefer to not do this sequentially by rank number but as communications arrive
            #

            # println("In Process(1): ",recv_id,"\n")

            # Receiving xi_low and yi_low #
            xi_low, yi_low = receive_inducing_points_in_global(ni_low, n_dims, recv_id, comm) ### TODO: IMPORTANT! Need to communicate which partition is being recieved

            # Update global parameters from local data #
            grads  = gradient(gp -> svgp_elbo_global(xi_low, yi_low, gp, global_svgp), global_svgp.params)[1][]

            Flux.Optimise.update!(opt_lrho, global_svgp.params.log_rho,       -grads.log_rho)
            Flux.Optimise.update!(opt_lkap, global_svgp.params.log_kappa,     -grads.log_kappa)
            Flux.Optimise.update!(opt_lsig, global_svgp.params.log_sigma,     -grads.log_sigma)
            Flux.Optimise.update!(optm,     global_svgp.params.inducing_mean, -grads.inducing_mean)
            Flux.Optimise.update!(optS,     global_svgp.params.inducing_L,    -LowerTriangular(grads.inducing_L) )
            Flux.Optimise.update!(optx,     global_svgp.params.inducing_locs, -grads.inducing_locs)
            # println("global ELBO optimized")

            # Sending prior_mean, prior_cov to local process ###################
            send_prior_to_local( pred_vgp(xi_low, global_svgp), recv_id, comm) ### TODO: IMPORTANT! Need to communicate which partition is being recieved
            # println("prior_mean, prior_cov sent to local")

        end # for loop

    end # myrank == root

    ###################### END OF INITIALIZATION ###############################

    ##################### OUTER LOOP ###########################################

    for iter= 1:n_iters
        if mod(iter, 10) == 0
            print(iter, " ")
        end

        if myrank != root

            for ii in 1:len ## TODO: Clean this up with Kelin to process non-sequentially
                            ## TODO: ENSURE THE MPI PROPERLY ASSOCIATED DATA COMMUNICATED TO CORRECT PARTITION
                # Step 4
                # Each local process will receive prior_mean and prior_cov

                # Receiving prior_mean, prior_cov to local process #################
                prior_mean, prior_sd = receive_prior_in_local(
                                                        ni_low,
                                                        root,
                                                        comm)
                # println("prior_mean, prior_cov received in local")

                # Step 5
                # Update local svgp using prior_mean prior_cov from global
                for sub_iter = 1:n_subcycle_ELBO_local
                    inds   = rand(1:size(local_svgps[ii].data.y)[1], batch_size)

                    grads  = gradient(gp -> svgp_elbo_local(
                            local_svgps[ii].data.x[inds,:],
                            local_svgps[ii].data.y[inds],
                            gp,
                            local_svgps[ii],
                            (prior_mean, prior_sd)
                            ), local_svgps[ii].params)[1][]

                    Flux.Optimise.update!(lopt_lrho[ii], local_svgps[ii].params.log_rho,       -grads.log_rho)
                    Flux.Optimise.update!(lopt_lkap[ii], local_svgps[ii].params.log_kappa,     -grads.log_kappa)
                    Flux.Optimise.update!(lopt_lsig[ii], local_svgps[ii].params.log_sigma,     -grads.log_sigma)
                    Flux.Optimise.update!(loptm[ii],     local_svgps[ii].params.inducing_mean, -grads.inducing_mean)
                    Flux.Optimise.update!(loptS[ii],     local_svgps[ii].params.inducing_L,    -LowerTriangular(grads.inducing_L) )
                    Flux.Optimise.update!(loptx[ii],     local_svgps[ii].params.inducing_locs, -grads.inducing_locs)
                    if sub_iter == n_subcycle_ELBO_local
                        # Collecting local ELBO in an array ###########################
                        trace_l[ii][iter] = svgp_elbo_local(
                            local_svgps[ii].data.x[inds,:],
                            local_svgps[ii].data.y[inds],
                            local_svgps[ii].params,
                            local_svgps[ii],
                            (prior_mean, prior_sd)
                        )
                    end
                end
                # println("local ELBO optimized")

                # Step 6
                # Sending inducing point data to global
                send_inducing_points_to_global(local_svgps[ii].params.inducing_locs, local_svgps[ii].params.inducing_mean, root, comm)
            end

        end # myrank != root

        if myrank == root
            step_elbo = 0.
            for part = (root+1):n_parts
                status  = MPI.Probe(MPI.MPI_ANY_SOURCE, MPI.MPI_ANY_TAG, comm)
                recv_id = status.source

                xi_low  = zeros(ni_low, n_dims)
                rmsg    = MPI.Recv!(xi_low, recv_id, 0, comm)
                yi_low  = zeros(ni_low)
                rmsg    = MPI.Recv!(yi_low, recv_id, 1, comm)s

                # Update global parameters from local data #
                for sub_iter = 1:n_subcycle_ELBO_global
                    grads  = gradient(gp -> svgp_elbo_global(xi_low, yi_low, gp, global_svgp), global_svgp.params)[1][]

                    Flux.Optimise.update!(opt_lrho, global_svgp.params.log_rho,       -grads.log_rho)
                    Flux.Optimise.update!(opt_lkap, global_svgp.params.log_kappa,     -grads.log_kappa)
                    Flux.Optimise.update!(opt_lsig, global_svgp.params.log_sigma,     -grads.log_sigma)
                    Flux.Optimise.update!(optm,     global_svgp.params.inducing_mean, -grads.inducing_mean)
                    Flux.Optimise.update!(optS,     global_svgp.params.inducing_L,    -LowerTriangular(grads.inducing_L) )
                    Flux.Optimise.update!(optx,     global_svgp.params.inducing_locs, -grads.inducing_locs)
                end
                # println("global ELBO optimized")

                # Update step_elbo
                step_elbo += svgp_elbo(xi_low, yi_low, global_svgp.params, global_svgp)

                # Sending prior_mean, prior_cov to local process ###################
                send_prior_to_local( pred_vgp(xi_low, global_svgp), recv_id, comm) # MAKE SURE COMMUNICATION INDICATES WHICH PARTITION IN A PROCESS THIS COMM IS FOR
                # println("prior_mean, prior_cov sent to local")

            end # for loop

            trace_h[iter] = step_elbo
        end # myrank == root

        MPI.Barrier(comm)
    end # outeriter

    ############################################################################
    # Returning trace_h and trace_l to main function
    if myrank != root

        for ii in 1:len
            send_trace(trace_l[ii], root ,comm)
        end

    end #myrank != root

    if myrank == root

        for part = (root+1):n_parts
            recv_id = 1 + (part - 1) % (mysize - 1)
            trace_l = receive_trace(n_iters, recv_id, comm)
            store_trace_l[:,part] = trace_l
        end # recv_id

        return_trace = hcat(trace_h, store_trace_l)
        # println(return_trace)

        # Write trace to a csv file
        df_return_trace = DataFrame(return_trace)
        CSV.write("./return_trace.csv", df_return_trace)

    end # myrank == root
    ############################################################################

    ############################################################################
    # Save svgp_params
    if myrank != root

        # send_svgp_params(local_svgp.params, root ,comm)
        for ii in 1:len
            save(string("local_fitted_parameters",part_split[myrank][ii],"_",mysize-1, ".jld2"), string("local_svgp",part_split[myrank][ii]), local_svgps[ii])
        end

    end #myrank != root

    if myrank == root

        # local_svgp_params = [receive_svgp_params(recv_id, comm) for recv_id = 1:mysize-1]

        save(string("global_fitted_svgp_",mysize-1,".jld2"), "global_svgp", global_svgp)

    end # myrank == root
    ############################################################################

    MPI.Finalize()

end     # mpifit_HSVGP

max_stars_repo_name: ayanbiswas/PRISM

max_issues_repo_name: ayanbiswas/PRISM

max_forks_repo_name: ayanbiswas/PRISM


text:
if !isdefined(Main, :Groebner)
    import Groebner
end

import AbstractAlgebra
using BenchmarkTools
using Logging
global_logger(ConsoleLogger(stderr, Logging.Error))

# BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10
BenchmarkTools.DEFAULT_PARAMETERS.samples = 2

function benchmark_system_my(system, tablesize)
    system = Groebner.change_ordering(system, :degrevlex)
    Groebner.groebner([system[1]], tablesize)

    println("length = $(length(system)), nvars = $(AbstractAlgebra.nvars(AbstractAlgebra.parent(system[1])))")
    println("maxdeg = $(maximum(AbstractAlgebra.total_degree ∘ AbstractAlgebra.leading_monomial, system))")

    @btime Groebner.groebner($system, $tablesize, reduced=true)
end

function run_f4_ff_degrevlex_benchmarks(ground)
    systems = [
    ("cyclic 10", Groebner.rootn(10, ground=ground)),
    ("cyclic 11", Groebner.rootn(11, ground=ground)),
    ("cyclic 12", Groebner.rootn(12, ground=ground)),
    ("cyclic 13", Groebner.rootn(13, ground=ground)),
        ("katsura 6", Groebner.katsura6(ground=ground)),
        ("katsura 9", Groebner.katsura9(ground=ground)),
        ("katsura 10", Groebner.katsura10(ground=ground)),
        ("noon 4"    ,Groebner.noonn(4, ground=ground)),
        ("noon 5"    ,Groebner.noonn(5, ground=ground)),
        ("noon 6"    ,Groebner.noonn(6, ground=ground)),
        ("noon 7"    ,Groebner.noonn(7, ground=ground)),
        ("eco 7"    ,Groebner.eco7(ground=ground)),
        ("eco 10"    ,Groebner.eco10(ground=ground)),
        ("eco 11"    ,Groebner.eco11(ground=ground)),
        ("ku 10"    ,Groebner.ku10(ground=ground)),
        ("s9_1"    ,Groebner.s9_1(ground=ground)),
        ("kinema"    ,Groebner.kinema(ground=ground)),
        ("ojika4_d1R2_d2R5"    ,Groebner.ojika4_d1R2_d2R5(ground=ground))
    ]

    println()

    for (name, system) in systems
        for tablesize in (2^13, 2^14, 2^15, 2^16)
            println("$name - $tablesize")
            benchmark_system_my(system, tablesize)
        end
        println("########\n")
    end
end

function big_run_f4_ff_degrevlex_benchmarks(ground)
    systems = [
        #("cyclic 14", Groebner.rootn(14, ground=ground)),
        #("cyclic 15", Groebner.rootn(15, ground=ground)),
        #("katsura 11", Groebner.katsura11(ground=ground)),
        #("katsura 12", Groebner.katsura12(ground=ground)),
        #("katsura 13", Groebner.katsura13(ground=ground)),
        ("noon 7"    ,Groebner.noonn(7, ground=ground)),
        ("noon 8"    ,Groebner.noonn(8, ground=ground)),
        ("eco 11"    ,Groebner.eco11(ground=ground)),
        ("eco 12"    ,Groebner.eco12(ground=ground)),
        #("eco 13"    ,Groebner.eco13(ground=ground)),
    ]

    println()

    for (name, system) in systems
        for tablesize in (2^12, 2^14, 2^16)
            println("$name - $tablesize")
            benchmark_system_my(system, tablesize)
        end
        println("########\n")
    end
end

println()
ground = AbstractAlgebra.GF(2^31 - 1)
run_f4_ff_degrevlex_benchmarks(ground)

#=

=#


text:
"""get_h5attr(f::AbstractString, attr::AbstractString)
Get an attribute from H5 file.
"""
function get_h5attr(f::AbstractString, attr::AbstractString)
    HDF5.h5open(f, "r") do h5io
        return HDF5.attributes(h5io)[attr] |> read
    end
end

function get_h5attr(f::AbstractString, symbol::Symbol)
    return get_h5attr(f, String(symbol))
end

"""get_h5data(f::AbstractString, path::AbstractString)
Get a dataset from H5 file.
"""
function get_h5data(f::AbstractString, path::AbstractString)
    HDF5.h5open(f, "r") do h5io
        return h5io[path] |> read
    end
end

function get_h5data(f::AbstractString, symbol::Symbol)
    return get_h5data(f, String(symbol))
end

"""
    readcsv2dataframe(f::AbstractString; kwargs...)
Read a CSV file on disk to a DataFrame object. kwargs are forwarded to `CSV.read`.
"""
function readcsv2dataframe(f::AbstractString; kwargs...)
    CSV.read(f, DataFrame; kwargs...)
end