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

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

**Constructors**

```julia
DataTable(columns::Vector{Any}, names::Vector{Symbol})
DataTable(kwargs...)
DataTable() # an empty DataTable
DataTable(t::Type, nrows::Integer, ncols::Integer) # an empty DataTable of arbitrary size
DataTable(column_eltypes::Vector, names::Vector, nrows::Integer)
DataTable(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, `DataTable(columns::Vector{Any},
names::Vector{Symbol})` does not convert to `NullableArray`.

A `DataTable` is a lightweight object. As long as columns are not
manipulated, creation of a DataTable 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 DataTable is not type stable. For
performance-critical code, do not index into a DataTable inside of
loops.

**Examples**

```julia
dt = DataTable()
v = ["x","y","z"][rand(1:3, 10)]
dt1 = DataTable(Any[collect(1:10), v, rand(10)], [:A, :B, :C])  # columns are Arrays
dt2 = DataTable(A = 1:10, B = v, C = rand(10))           # columns are NullableArrays
dump(dt1)
dump(dt2)
describe(dt2)
DataTables.head(dt1)
dt1[:A] + dt2[:C]
dt1[1:4, 1:2]
dt1[[:A,:C]]
dt1[1:2, [:A,:C]]
dt1[:, [:A,:C]]
dt1[:, [1,3]]
dt1[1:4, :]
dt1[1:4, :C]
dt1[1:4, :C] = 40. * dt1[1:4, :C]
[dt1; dt2]  # vcat
[dt1  dt2]  # hcat
size(dt1)
```

"""
type DataTable <: AbstractDataTable
    columns::Vector{Any}
    colindex::Index

    function DataTable(columns::Vector{Any}, colindex::Index)
        if length(columns) == length(colindex) == 0
            return new(Vector{Any}(0), Index())
        elseif length(columns) != length(colindex)
            throw(DimensionMismatch("Number of columns ($(length(columns))) and number of column names ($(length(colindex))) are not equal"))
        end
        lengths = [isa(col, AbstractArray) ? length(col) : 1 for col in columns]
        minlen, maxlen = extrema(lengths)
        if minlen == 0 && maxlen == 0
            return new(columns, colindex)
        elseif minlen != maxlen || minlen == maxlen == 1
            # recycle scalars
            for i in 1:length(columns)
                isa(columns[i], AbstractArray) && continue
                columns[i] = fill(columns[i], maxlen)
                lengths[i] = maxlen
            end
            uls = unique(lengths)
            if length(uls) != 1
                strnames = string.(names(colindex))
                estrings = ["column length $u for column(s) " *
                            join(strnames[lengths .== u], ", ", " and ") for (i, u) in enumerate(uls)]
                throw(DimensionMismatch(join(estrings, " is incompatible with ", ", and is incompatible with ")))
            end
        end
        for (i, c) in enumerate(columns)
            if isa(c, Range)
                columns[i] = collect(c)
            elseif !isa(c, AbstractVector)
                throw(DimensionMismatch("columns must be 1-dimensional"))
            end
        end
        new(columns, colindex)
    end
end

function DataTable(; kwargs...)
    colnames = Symbol[k for (k,v) in kwargs]
    columns = Any[v for (k,v) in kwargs]
    DataTable(columns, Index(colnames))
end

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

# Initialize an empty DataTable with specific eltypes and names
function DataTable{T<:Type}(column_eltypes::AbstractVector{T}, cnames::AbstractVector{Symbol}, nrows::Integer)
    numcols = length(column_eltypes)
    columns = Vector{Any}(numcols)
    for j in 1:numcols
        elty = column_eltypes[j]
        if elty <: Nullable
            if eltype(elty) <: CategoricalValue
                columns[j] = NullableCategoricalArray{eltype(elty)}(nrows)
            else
                columns[j] = NullableVector{eltype(elty)}(nrows)
            end
        else
            if elty <: CategoricalValue
                columns[j] = CategoricalVector{elty}(nrows)
            else
                columns[j] = Vector{elty}(nrows)
            end
        end
    end
    return DataTable(columns, Index(convert(Vector{Symbol}, cnames)))
end

# Initialize an empty DataTable with specific eltypes and names
# and whether a nominal array should be created
function DataTable{T<:Type}(column_eltypes::AbstractVector{T}, cnames::AbstractVector{Symbol},
                            nominal::Vector{Bool}, nrows::Integer)
    # upcast Vector{DataType} -> Vector{Type} which can hold CategoricalValues
    updated_types = convert(Vector{Type}, column_eltypes)
    for i in eachindex(nominal)
        nominal[i] || continue
        if updated_types[i] <: Nullable
            updated_types[i] = Nullable{CategoricalValue{eltype(updated_types[i])}}
        else
            updated_types[i] = CategoricalValue{updated_types[i]}
        end
    end
    return DataTable(updated_types, cnames, nrows)
end

# Initialize empty DataTable objects of arbitrary size
function DataTable(t::Type, nrows::Integer, ncols::Integer)
    return DataTable(fill(t, ncols), nrows)
end

# Initialize an empty DataTable with specific eltypes
function DataTable{T<:Type}(column_eltypes::AbstractVector{T}, nrows::Integer)
    return DataTable(column_eltypes, gennames(length(column_eltypes)), nrows)
end

##############################################################################
##
## AbstractDataTable interface
##
##############################################################################

index(dt::DataTable) = dt.colindex
columns(dt::DataTable) = dt.columns

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

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

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

@compat const ColumnIndex = Union{Real, Symbol}

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

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

# dt[:] => DataTable
Base.getindex(dt::DataTable, col_inds::Colon) = copy(dt)

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

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

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

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

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

# dt[SingleRowIndex, :] => DataTable
Base.getindex(dt::DataTable, row_ind::Real, col_inds::Colon) = dt[[row_ind], col_inds]

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

# dt[:, :] => DataTable
Base.getindex(dt::DataTable, ::Colon, ::Colon) = copy(dt)

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

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

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

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

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

function insert_single_entry!(dt::DataTable, v::Any, row_ind::Real, col_ind::ColumnIndex)
    if haskey(index(dt), col_ind)
        dt.columns[index(dt)[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}(dt::DataTable,
                                            v::Any,
                                            row_inds::AbstractVector{T},
                                            col_ind::ColumnIndex)
    if haskey(index(dt), col_ind)
        dt.columns[index(dt)[col_ind]][row_inds] = v
        return v
    else
        error("Cannot assign to non-existent column: $col_ind")
    end
end

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# dt[:] = DataTable, dt[:, :] = DataTable
function Base.setindex!(dt::DataTable,
                                  new_dt::DataTable,
                                  row_inds::Colon,
                                  col_inds::Colon=Colon())
    dt.columns = copy(new_dt.columns)
    dt.colindex = copy(new_dt.colindex)
    dt
end

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

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

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

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

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

Base.empty!(dt::DataTable) = (empty!(dt.columns); empty!(index(dt)); dt)

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

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

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

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

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

function Base.merge!(dt::DataTable, others::AbstractDataTable...)
    for other in others
        for n in _names(other)
            dt[n] = other[n]
        end
    end
    return dt
end

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

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

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

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

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

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

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

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

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

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

# hcat! for 2 arguments
function hcat!(dt1::DataTable, dt2::AbstractDataTable)
    u = add_names(index(dt1), index(dt2))
    for i in 1:length(u)
        dt1[u[i]] = dt2[i]
    end
    return dt1
end
hcat!(dt::DataTable, x::AbstractVector) = hcat!(dt, DataTable(Any[x]))

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

# hcat
Base.hcat(dt::DataTable, x) = hcat!(copy(dt), x)
Base.hcat(dt1::DataTable, dt2::AbstractDataTable) = hcat!(copy(dt1), dt2)
Base.hcat(dt1::DataTable, dt2::AbstractDataTable, dtn::AbstractDataTable...) = hcat!(hcat(dt1, dt2), dtn...)

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

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

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

function categorical!(dt::DataTable, cname::Union{Integer, Symbol}, compact::Bool=true)
    dt[cname] = categorical(dt[cname], compact)
    dt
end

function categorical!{T <: Union{Integer, Symbol}}(dt::DataTable, cnames::Vector{T},
                                                   compact::Bool=true)
    for cname in cnames
        dt[cname] = categorical(dt[cname], compact)
    end
    dt
end

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

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

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

function Base.convert(::Type{DataTable}, d::Associative)
    colnames = keys(d)
    if isa(d, Dict)
        colnames = sort!(collect(keys(d)))
    else
        colnames = keys(d)
    end
    colindex = Index(Symbol[k for k in colnames])
    columns = Any[d[c] for c in colnames]
    DataTable(columns, colindex)
end

##############################################################################
##
## push! a row onto a DataTable
##
##############################################################################

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

function Base.push!(dt::DataTable, associative::Associative)
    i = 1
    for nm in _names(dt)
        try
            val = get(() -> associative[string(nm)], associative, nm)
            push!(dt[nm], val)
        catch
            #clean up partial row
            for j in 1:(i - 1)
                pop!(dt[_names(dt)[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!(dt::DataTable, iterable::Any)
    if length(iterable) != length(dt.columns)
        msg = "Length of iterable does not match DataTable column count."
        throw(ArgumentError(msg))
    end
    i = 1
    for t in iterable
        try
            push!(dt.columns[i], t)
        catch
            #clean up partial row
            for j in 1:(i - 1)
                pop!(dt.columns[j])
            end
            msg = "Error adding $t to column :$(_names(dt)[i]). Possible type mis-match."
            throw(ArgumentError(msg))
        end
        i += 1
    end
end


text:
##==============================================================================
## LEGACY SUPPORT FOR ZMQ IN CAESAR
##==============================================================================

@deprecate _evalType(pt::String) DFG.getTypeFromSerializationModule(pt)
# # FIXME, replace with T = DFG.getTypeFromSerializationModule(obj._type)
# export _evalType
# # not sure if and where this is still being used
# function _evalType(pt::String)::Type
#   @error "_evalType has been deprecated, use DFG serialization methods instead."
#   try
#     getfield(Main, Symbol(pt))
#   catch ex
#     io = IOBuffer()
#     showerror(io, ex, catch_backtrace())
#     err = String(take!(io))
#     error("_evalType: Unable to locate factor/distribution type '$pt' in main context (e.g. do a using on all your factor libraries). Please check that this factor type is loaded into main. Stack trace = $err")
#   end
# end

##==============================================================================
## Deprecate code below before v0.27
##==============================================================================

# """
#     $SIGNATURES

# Sample points from a regular grid scalar field level set.

# Notes
# - modifies first argument roi to save on memory allocations
# - user should duplicate if with a deepcopy if needed
# """
# function sampleLevelSetGaussian!( roi::AbstractMatrix{<:Real},
#                                   sigma::Real,
#                                   x_grid::AbstractVector{<:Real},
#                                   y_grid::AbstractVector{<:Real};
#                                   sigma_scale::Real=3  )
#   #
#   # make Gaussian
#   roi .^= 2
#   roi .*= 0.5/(sigma^2)
#     # truncate at sigma_scale*sigma
#     _roi = thres .- roi

#   sampleHeatmap(roi, x_grid, y_grid, sigma_scale^2)
# end

@deprecate HeatmapDensityRegular(w...;kw...) LevelSetGridNormal(w...;kw...)

# @deprecate getLevelSetSigma(data::AbstractMatrix{<:Real}, level::Real, w...; kw...) sampleLevelSetGaussian!(data.-level, w...; kw...)

getOutNeighbors(w...;kw...) = error("Obsolete, use DFG.getNeighbors instead.")

# # TODO -- there should be a better way, without retrieving full vertex
# # TODO -- Deprecated for DFG -- must update
# function getOutNeighbors(dfg::T, v::V; needdata::Bool=false, ready::Union{Nothing, Int}=nothing, backendset::Union{Nothing, Int}=nothing)::Vector{Symbol} where {T <: AbstractDFG, V <: DFGNode}
#   @warn "TODO: needdata is currently ignored. Symbols are returned."
#   nodes = getNeighbors(dfg, v, ready=ready, backendset=backendset)
#   return nodes
# end
# function getOutNeighbors(dfg::T, vertSym::Symbol; needdata::Bool=false, ready::Int=1, backendset::Int=1 )::Vector{Symbol} where {T <: AbstractDFG, V <: DFGNode}
#   @warn "TODO: needdata is currently ignored. Symbols are returned."
#   nodes = getNeighbors(dfg, vertSym, ready=ready, backendset=backendset)
#   return nodes
# end

# SolverParams(;dimID::Int=0,
#               registeredModuleFunctions=nothing,
#               reference=nothing,
#               stateless::Bool=false,
#               qfl::Int=99999999999,
#               isfixedlag::Bool=false,
#               limitfixeddown::Bool=false,
#               incremental::Bool=true,
#               useMsgLikelihoods::Bool=false,
#               upsolve::Bool=true,
#               downsolve::Bool=true,
#               drawtree::Bool=false,
#               drawCSMIters::Bool=true,
#               showtree::Bool=false,
#               drawtreerate::Float64=0.5,
#               dbg::Bool=false,
#               async::Bool=false,
#               limititers::Int=500,
#               N::Int=100,
#               multiproc::Bool=1 < nprocs(),
#               logpath::String="/tmp/caesar/$(now())",
#               graphinit::Bool=true,
#               treeinit::Bool=false,
#               limittreeinit_iters::Int=10,
#               algorithms::Vector{Symbol}=[:default],
#               spreadNH::Real=3.0,
#               inflation::Real=5.0,
#               inflateCycles::Int=3,
#               gibbsIters::Int=3,
#               maxincidence::Int=500,
#               alwaysFreshMeasurements::Bool=true,
#               attemptGradients::Bool=true,
#               devParams::Dict{Symbol,String}=Dict{Symbol,String}()
#             ) = begin useMsgLikelihoods==true && @warn "useMsgLikelihoods is under development, use with care, see #1010"
#                 SolverParams( dimID,
#                               registeredModuleFunctions,
#                               reference,
#                               stateless,
#                               qfl,
#                               isfixedlag,
#                               limitfixeddown,
#                               incremental,
#                               useMsgLikelihoods,
#                               upsolve,
#                               downsolve,
#                               drawtree,
#                               drawCSMIters,
#                               showtree,
#                               drawtreerate,
#                               dbg,
#                               async,
#                               limititers,
#                               N,
#                               multiproc,
#                               logpath,
#                               graphinit,
#                               treeinit,
#                               limittreeinit_iters,
#                               algorithms,
#                               spreadNH,
#                               inflation,
#                               inflateCycles,
#                               gibbsIters,
#                               maxincidence,
#                               alwaysFreshMeasurements,
#                               attemptGradients,
#                               devParams )
#             end
#

@deprecate setVariableInferDim!(w...;kw...) setIPC!(w...;kw...)

# moved to IncrementalInference/attic
# include("CSMOccuranceUtils.jl")

@deprecate prepgenericconvolution(w...;kw...) _prepCCW(w...;kw...)

function getSample(cf::CalcFactor, N::Int)
  # Base.depwarn("`getSample(cf, N)` is deprecated, use `getSample(cf)`", :getSample)
  error("`getSample(cf, N)` is deprecated, use `getSample(cf)`")
end

# """
#     $SIGNATURES

# Build an approximate density `[Y|X,DX,.]=[X|Y,DX][DX|.]` as proposed by the conditional convolution.

# Notes
# - Assume both are on circular manifold, `manikde!(pts, (:Circular,))`
# """
# function approxConvCircular(pX::ManifoldKernelDensity,
#                             pDX::ManifoldKernelDensity; N::Int=100)
#   #

#   # building basic factor graph
#   tfg = initfg()
#   addVariable!(tfg, :s1, Sphere1)
#   addVariable!(tfg, :s2, Sphere1)
#   addFactor!(tfg, [:s1;:s2], Sphere1Sphere1(pDX), graphinit=false)
#   initManual!(tfg,:s1, pX)

#   # solve for outgoing proposal value
#   approxConv(tfg,:s1s2f1,:s2)
# end

# function approxConvCircular(pX::ManifoldKernelDensity,
#                             pDX::SamplableBelief; N::Int=100)
#   #
#   pts = reshape(rand(pDX, N), 1, :)
#   pC = manikde!(pts, Sphere1)
#   approxConvCircular(pX, pC)
# end

# function approxConvCircular(pX::SamplableBelief,
#                             pDX::ManifoldKernelDensity; N::Int=100)
#   #
#     pts = reshape(rand(pX, N), 1, :)
#   pC = manikde!(pts, Sphere1)
#   approxConvCircular(pC, pDX)
# end

@deprecate printGraphSummary(dfg::AbstractDFG, logger=ConsoleLogger()) show(logger.stream, MIME("text/plain"), dfg)
@deprecate printSummary(dfg::AbstractDFG, logger=ConsoleLogger()) show(logger.stream, MIME("text/plain"), dfg)

# """
#     $SIGNATURES

# Print basic summary of graph to `logger=ConsoleLogger()`.
# """
# function printGraphSummary(dfg::G, logger=ConsoleLogger())::Nothing where {G <: AbstractDFG}
#     vars = ls(dfg)
#     fcts = lsf(dfg)

#     prio = lsfPriors(dfg)

#     isinit = map(x->isInitialized(dfg,x), vars)
#     infdim = map(x->getVariableInferredDim(dfg, x), vars)
#     numedges = map(v->length(ls(dfg, v)), vars)
#     numfed = map(fc->length(ls(dfg, fc)), fcts)
#     vardims = map(v->getDimension(getVariable(dfg, v)), vars)
#     fctdims = map(v->getDimension(getFactor(dfg, v)), fcts)
#     priodims = map(v->getDimension(getFactor(dfg, v)), prio)

#     with_logger(logger) do
#       @info "Distributed Factor Graph summary:"
#       @info "  num variables:    $(length(vars))"
#       @info "  num factors:      $(length(fcts)), w/ $(length(prio)) priors"
#       @info "  var initialized:  $(sum(isinit))"
#       @info ""
#       @info "  var num edges: min. $(minimum(numedges)) | mean $(round(Statistics.mean(numedges),digits=2)) | 90% $(round(quantile(numedges,0.9),digits=2)) | max. $(maximum(numedges))"
#       @info "  fct num edges: min. $(minimum(numfed)) | mean $(round(Statistics.mean(numfed),digits=2)) | 90% $(round(quantile(numfed,0.9),digits=2)) | max. $(maximum(numfed))"
#       @info "  Variable dims: min. $(minimum(vardims)) | mean $(round(Statistics.mean(vardims),digits=2)) | 90% $(round(quantile(vardims,0.9),digits=2)) | max. $(maximum(vardims))"
#       @info "  Factor dims:   min. $(minimum(fctdims)) | mean $(round(Statistics.mean(fctdims),digits=2)) | 90% $(round(quantile(fctdims,0.9),digits=2)) | max. $(maximum(fctdims))"
#       @info "  Prior dimens:  min. $(minimum(priodims)) | mean $(round(Statistics.mean(priodims),digits=2)) | 90% $(round(quantile(priodims,0.9),digits=2)) | max. $(maximum(priodims))"
#       @info "  var infr'dims: min. $(minimum(infdim)) | mean $(round(Statistics.mean(infdim),digits=2)) | 90% $(round(quantile(infdim,0.9),digits=2)) | max. $(maximum(infdim))"
#     end
#     nothing
# end

# """
#     $SIGNATURES

# Print basic summary of graph to `logger=ConsoleLogger()`.
# """
# function printSummary(dfg::G, logger=ConsoleLogger()) where G <: AbstractDFG
#     printGraphSummary(dfg, logger)
# end

# """
#     $SIGNATURES

# Calculate both measured and predicted relative variable values, starting with `from` at zeros up to `to::Symbol`.

# Notes
# - assume single variable separators only.

# DevNotes
# - TODO better consolidate with [`approxConvBelief`](@ref) which can now also work with factor chains.
# """
# function accumulateFactorChain( dfg::AbstractDFG,
#                                 from::Symbol,
#                                 to::Symbol,
#                                 fsyms::Vector{Symbol}=findFactorsBetweenNaive(dfg, from, to);
#                                 initval=zeros(size(getVal(dfg, from))))

#   # get associated variables
#   svars = union(ls.(dfg, fsyms)...)

#   # use subgraph copys to do calculations
#   tfg_meas = buildSubgraph(dfg, [svars;fsyms])
#   tfg_pred = buildSubgraph(dfg, [svars;fsyms])

#   # drive variable values manually to ensure no additional stochastics are introduced.
#   nextvar = from
#   initManual!(tfg_meas, nextvar, initval)
#   initManual!(tfg_pred, nextvar, initval)

#   # nextfct = fsyms[1] # for debugging
#   for nextfct in fsyms
#     nextvars = setdiff(ls(tfg_meas,nextfct),[nextvar])
#     @assert length(nextvars) == 1 "accumulateFactorChain requires each factor pair to separated by a single variable"
#     nextvar = nextvars[1]
#     meas, pred = approxDeconv(dfg, nextfct) # solveFactorMeasurements
#     pts_meas = approxConv(tfg_meas, nextfct, nextvar, (meas,ones(Int,100),collect(1:100)))
#     pts_pred = approxConv(tfg_pred, nextfct, nextvar, (pred,ones(Int,100),collect(1:100)))
#     initManual!(tfg_meas, nextvar, pts_meas)
#     initManual!(tfg_pred, nextvar, pts_pred)
#   end
#   return getVal(tfg_meas,nextvar), getVal(tfg_pred,nextvar)
# end

# # TODO should this be consolidated with regular approxConv?
# # TODO, perhaps pass Xi::Vector{DFGVariable} instead?
# function approxConvBinary(arr::Vector{Vector{Float64}},
#                           meas::AbstractFactor,
#                           outdims::Int,
#                           fmd::FactorMetadata,
#                           measurement::Tuple=(Vector{Vector{Float64}}(),);
#                           varidx::Int=2,
#                           N::Int=length(arr),
#                           vnds=DFGVariable[],
#                           _slack=nothing )
#   #
#   # N = N == 0 ? size(arr,2) : N
#   pts = [zeros(outdims) for _ in 1:N];
#   ptsArr = Vector{Vector{Vector{Float64}}}()
#   push!(ptsArr,arr)
#   push!(ptsArr,pts)

#   fmd.arrRef = ptsArr

#   # TODO consolidate with ccwl??
#   # FIXME do not divert Mixture for sampling
#   # cf = _buildCalcFactorMixture(ccwl, fmd, 1, ccwl.measurement, ARR) # TODO perhaps 0 is safer
#   # FIXME 0, 0, ()
#   cf = CalcFactor( meas, fmd, 0, 0, (), ptsArr)

#   measurement = length(measurement[1]) == 0 ? sampleFactor(cf, N) : measurement
#   # measurement = size(measurement[1],2) == 0 ? sampleFactor(meas, N, fmd, vnds) : measurement

#   zDim = length(measurement[1][1])
#   ccw = CommonConvWrapper(meas, ptsArr[varidx], zDim, ptsArr, fmd, varidx=varidx, measurement=measurement)  # N=> size(measurement[1],2)

#   for n in 1:N
#     ccw.cpt[Threads.threadid()].particleidx = n
#     _solveCCWNumeric!( ccw, _slack=_slack )
#   end
#   return pts
# end

@deprecate getParametricMeasurement(w...;kw...) getMeasurementParametric(w...;kw...)

# function prtslperr(s)
#   println(s)
#   sleep(0.1)
#   error(s)
# end

@deprecate getKDE(w...;kw...) getBelief(w...;kw...)

@deprecate findRelatedFromPotential(w...;kw...) (calcProposalBelief(w...;kw...),)

# function generateNullhypoEntropy( val::AbstractMatrix{<:Real},
#                                   maxlen::Int,
#                                   spreadfactor::Real=10  )
#   #
#   # covD = sqrt.(vec(Statistics.var(val,dims=2))) .+ 1e-3
#   # cVar = diagm((spreadfactor*covD).^2)
#   len = size(val, 1)
#   cVar = diagm((spreadfactor*ones(len)).^2)
#   mu = zeros( len )
#   MvNormal( mu, cVar )
# end

# @deprecate calcFactorResidualTemporary( fct::AbstractRelative, measurement::Tuple, T_param_args...; kw...) calcFactorResidualTemporary(fct, measurement, tuple(((x->x[1]).(T_param_args))...), tuple(((x->x[2]).(T_param_args))...); kw...)

# @deprecate _buildGraphByFactorAndTypes!(fct::AbstractFactor, TypeParams_vec...; kw...) _buildGraphByFactorAndTypes!(fct, (x->x[1]).(TypeParams_vec), (x->x[2]).(TypeParams_vec) ; kw...)

# @deprecate _evalFactorTemporary!( fct::AbstractFactor, sfidx::Int, meas_single::Tuple, TypeParams_args...; kw...) _evalFactorTemporary!( fct,sfidx,meas_single, (x->x[1]).(TypeParams_args), (x->x[2]).(TypeParams_args); kw...)

# MOVED TO AMP v0.4.6
# # NOTE SWITCHED TO ON-MANIFOLD VERSION, but this used to give deviation
# calcVariableCovarianceBasic(M::AbstractManifold, ptsArr::Vector{Vector{Float64}})
#   # cannot calculate the stdev from uninitialized state
#   # FIXME assume point type is only Vector{Float} at this time
#   @cast arr[i,j] := ptsArr[j][i]
#   msst = Statistics.std(arr, dims=2)
#   # FIXME use adaptive scale, see #802
#   msst_ = 0 < sum(1e-10 .< msst) ? maximum(msst) : 1.0
#   return msst_
# end

@deprecate testFactorResidualBinary(fct, T1::InstanceType{InferenceVariable}, T2::InstanceType{InferenceVariable}, param1, param2, meas) calcFactorResidualTemporary(fct, (T1,T2), meas, (param1,param2))

# """
#     $SIGNATURES
#
# Return the directly achievable dimensionality of solve for each variable in a clique.
#
# Related
#
# getFactorSolvableDim
# """
# function getCliqVarPossibleDim(dfg::G,
#                                     cliq::TreeClique,
#                                     saturate::Bool=true,
#                                     fraction::Bool=true  )::Dict{Symbol, Float64}
#   #
#   # variables and factors associated with this clique
#   vars = getCliqAllVarIds(cliq)
#   # fcts = getCliqAllFactIds(cliq)
#   # rows = length(fcts)
#   cols = length(vars)
#
#   # for output result
#   dict = Dict{Symbol,Float64}()
#
#   for j in 1:cols
#     dict[vars[j]] += getVariablePossibleDim(dfg, vars[j])
#   end
#
# end

# getManifold(::InstanceType{LinearRelative{N}}) where {N} = Euclidean(N)
# getManifolds(::T) where {T <: LinearRelative} = convert(Tuple, getManifold(T))
# getManifolds(::Type{<:T}) where {T <: LinearRelative} = convert(Tuple, getManifold(T))
# getManifolds(fctType::Type{LinearRelative}) = getManifolds(getDomain(fctType))

# # FIXME, why is Manifolds depdendent on the solveKey?? Should just be at DFGVariable level?

# getManifolds(vd::VariableNodeData) = getVariableType(vd) |> getManifolds
# getManifolds(::DFGVariable{T}) where T <: InferenceVariable = getManifolds(T)
# getManifolds(dfg::AbstractDFG, sym::Symbol) = getManifolds(getVariable(dfg, sym))
#
# getManifolds(vartype::InferenceVariable) = vartype.manifolds
# getManifolds(vartype::Type{<: InferenceVariable}) = getManifolds(vartype())

# import DistributedFactorGraphs: getfnctype
# # TODO: Refactor - was is das?
# function getfnctype(data::GenericFunctionNodeData)
#   if typeof(data).name.name == :VariableNodeData
#     return VariableNodeData
#   end
#   return data.fnc.usrfnc!
# end
#
# function getfnctype(fact::DFGFactor; solveKey::Symbol=:default)
#   data = getData(fact) # TODO , solveKey=solveKey)
#   return getfnctype(data)
# end
#
# function getfnctype(dfg::T, lbl::Symbol; solveKey::Symbol=:default) where T <: AbstractDFG
#   getfnctype(getFactor(dfg, exvertid))
# end

# function manikde!(ptsArr::Vector{Vector{Float64}}, manis::Tuple)
#   arr = Matrix{Float64}(undef, length(ptsArr[1]), length(ptsArr))
#   @cast arr[i,j] = ptsArr[j][i]
#   manikde!(arr, manis)
# end

##==============================================================================
## Deprecate code below before v0.26
##==============================================================================

export listSolvekeys

# wow, that was quite far off -- needs testing
# function factorCanInitFromOtherVars(dfg::T,
#                                     fct::Symbol,
#                                     loovar::Symbol)::Tuple{Bool, Vector{Symbol}, Vector{Symbol}} where T <: AbstractDFG
#   #
#   # all variables attached to this factor
#   varsyms = getNeighbors(dfg, fct)
#
#   # list of factors to use in init operation
#   useinitfct = Symbol[]
#   faillist = Symbol[]
#   for vsym in varsyms
#     xi = DFG.getVariable(dfg, vsym)
#     if (isInitialized(xi) && sum(useinitfct .== fct) == 0 ) || length(varsyms) == 1
#       push!(useinitfct, fct)
#     end
#   end
#
#   return (length(useinitfct)==length(varsyms)&&length(faillist)==0,
#           useinitfct,
#           faillist   )
# end

@deprecate calcPerturbationFromVariable( fgc::FactorGradientsCached!, fromVar::Int, infoPerCoord::AbstractVector;tol::Real=0.02*fgc._h ) calcPerturbationFromVariable(fgc, [fromVar => infoPerCoord;], tol=tol )

function getVal(vA::Vector{<:DFGVariable}, solveKey::Symbol=:default)
  @error "getVal(::Vector{DFGVariable}) is obsolete, use getVal.(DFGVariable) instead."
  # len = length(vA)
  # vals = Array{Array{Float64,2},1}()
  # cols = Array{Int,1}()
  # push!(cols,0)
  # rows = Array{Int,1}()
  # for v in vA
  #     push!(vals, getVal(v, solveKey=solveKey))
  #     c = size(vals[end],2)
  #     r = size(vals[end],1)
  #     push!(cols, floor(Int,c))
  #     push!(rows, floor(Int,r))
  # end
  # cols = cumsum(cols)
  # sc = cols[end]
  # rw = floor(Int,rows[1])
  # val = Array{Float64,2}(undef,rw, sc)
  # for i in 1:(len-1)
  #     val[:,(cols[i]+1):cols[i+1]] = vals[i]
  # end
  # val[:,(cols[len]+1):cols[len+1]] = vals[len] # and the last one
  # return val
end

# """
#     $SIGNATURES

# Get graph node (variable or factor) dimension.
# """
# DFG.getDimension(vartype::InferenceVariable) = vartype.dims #TODO Deprecate
# DFG.getDimension(vartype::Type{<:InferenceVariable}) = getDimension(vartype())

@deprecate ensureAllInitialized!(w...;kw...) initAll!(w...;kw...)

@deprecate getFactorMean(w...) IIF.getMeasurementParametric(w...)[1]

# """
#     $SIGNATURES

# Recover the mean (Gaussian) or estimate stochastic mean (non-Gaussian) value stored in a factor measurement.

# Related

# accumulateFactorMeans, solveBinaryFactorParameteric
# """
# function getFactorMean(fct::FunctorInferenceType)
#   fctt = typeof(fct)
#   error("no getFactorMean defined for $(fctt.name), has fields $(fieldnames(fctt))")
# end

# getFactorMean(fct::Normal) = [fct.μ]
# getFactorMean(fct::MvNormal) = fct.μ
# getFactorMean(fct::Union{<:BallTreeDensity,<:ManifoldKernelDensity}) = getKDEMean(fct)
# getFactorMean(fct::AliasingScalarSampler) = Statistics.mean(rand(fct,1000))

# getFactorMean(fct::DFGFactor) = getFactorMean(getFactorType(fct))

# getFactorMean(dfg::AbstractDFG, fctsym::Symbol) = getFactorMean(getFactor(dfg, fctsym))

# AMP.getManifolds(::T) where {T <: InferenceVariable} = getManifolds(getManifold(T))
# AMP.getManifolds(::Type{T}) where {T <: InferenceVariable} = getManifolds(getManifold(T))

# getManifolds(getManifold(T))
getManifolds(::InstanceType{T}) where {T <: Union{InferenceVariable, AbstractFactor}} = error("getManifolds is obsolete, use getManifold(...)::MB.AbstactManifold instead.")

#

ext: jl


text:
# TODO: replace lnLambda with exact version. Use CODEs
# TODO: make tests
# TODO: Verify nus, nupar, etc. against H&S
# TODO: GEt rid of nuhatvta3 so that some temperature is defined
# TODO: populate solveCollisionalEquation function

# Contains various routines related to collision operators
module collisions
using ..constants
using ..species
using ..grids

export slowingdowntime, nu_s, nu_par, collisionOperator, lnLambda, nuhat, gfunc, nu_D

"""
nu_s(v::Real, a::SpeciesData, b::SpeciesData)

Returns νₛᵃᵇ(v) from Helander and Sigmar (2002), equation 3.46 (without the mass ratio factor that is cancelled in 3.40), where a and b are the test particle and target particle respectively
"""
function nu_s(v::Float64,a::SpeciesData,b::SpeciesData)
  vtb = sqrt(2.0 * b.temperature / b.mass)

  return nuhatvta3(a,b) * (a.mass/b.temperature) * Goverx(v/vtb)/vtb
end

function nu_s(v::Vector,a::SpeciesData,b::SpeciesData)
  vtb = sqrt(2.0 * b.temperature / b.mass)
  return nuhatvta3(a,b)* (a.mass/b.temperature) * Goverx(v/vtb)/vtb
end

"""
nu_par(v::Real, a::SpeciesData, b::SpeciesData)

Returns ν∥ᵃᵇ(v) from Helander and Sigmar (2002), equation 3.47, where a and b are the test particle and target particle respectively
"""
function nu_par(v::Float64,a::SpeciesData,b::SpeciesData)
  vtb = sqrt(2.0 * b.temperature / b.mass)
  if abs(v) .> eps(Float64)
     return 0.0
  else
     return nuhatvta3(a,b)* 2.0 * gfunc(v/vtb) / v^3
  end
end

function nu_par(v::Vector,a::SpeciesData,b::SpeciesData)
  vtb = sqrt(2.0 * b.temperature / b.mass)
  y=zeros(length(v))
  idx=abs(v) .> eps(Float64)
  y[idx] = nuhatvta3(a,b)* 2.0 * gfunc(v[idx]/vtb) ./ v[idx].^3
  return y
end

function nu_D(v::Vector,a::SpeciesData,b::SpeciesData)
  vtb = sqrt(2.0 * b.temperature / b.mass)
  vta = sqrt(2.0 * a.temperature / a.mass)
  y=zeros(length(v))
  idx=abs(v) .> eps(Float64)
  y[idx] = nuhatvta3(a,b)* (erf(v[idx]/vtb) - gfunc(v[idx]/vtb)) ./ (v[idx]/vta).^3
  return y
end

"""
Chandresekhar G function
"""
function gfunc(x::Float64)
   if abs(x) > eps(Float64)
      return (erf(x)-(2.0/sqrt(pi))*x.*exp(-x.^2))./(2.0*x.^2)
   else
      return 0.0
   end
end
function gfunc(x::Vector)
   y=zeros(length(x))
   idx=abs(x) .> eps(Float64)
   y[idx] = (erf(x[idx])-(2.0/sqrt(pi))*x[idx].*exp(-x[idx].^2))./(2.0*x[idx].^2)
   return y
end

function Goverx(x::Float64)
   if abs(x) > eps(Float64)
      return (erf(x)-(2.0/sqrt(pi))*x.*exp(-x.^2))./(2.0*x.^3)
   else
      return 2.0/(3.0*sqrt(pi))
   end
end
function Goverx(x::Vector)
   y=zeros(length(x))
   idx=abs(x) .> eps(Float64)
   y[idx] = (erf(x[idx])-(2.0/sqrt(pi))*x[idx].*exp(-x[idx].^2))./(2.0*x[idx].^3)
   y[abs(x) .<= eps(Float64)] = 2.0/(3.0*sqrt(pi))
   return y
end

"""
collisionOperator(vgrid::Grid, a::SpeciesData, b::SpeciesData)

Constructs a finite volume form of the collision operator
"""
function collisionOperator(vgrid::Grid, a::SpeciesData, b::SpeciesData)
   Nv = length(vgrid.ctr)

   if false
   # Replace a zero element with epsilon*next largest element
   vgrid[abs(vgrid).<eps(Float64)] = eps(Float64)*sort(abs(vgrid))[2]

   if ddv_in == [0.0]
      size(ddv_in) == (Nv,Nv) || error("ddv_in must be a square matrix the same size as v")
      ddv = ddv_in
   else
      ddv = createDifferentiationMatrix(vgrid)
   end

   return vgrid.^(-2) * ddv * ( nu_s(vgrid,a,b).*vgrid.^3  + 0.5*nu_par(vgrid,a,b).*vgrid.^4*ddv)
   end

   # Define functions that are required by operator constructor
   jacobian(v) = v.^2
   nu_s_term(v) = nu_s(v,a,b).*v
   nu_par_term(v) = 0.5*nu_par(v,a,b).*v.^2
   zerofunc(v) = 0.0
   # Since boundary conditions for collision operator are zero, we can worry about the source term later
   source = zeros(Nv)

   collop = createAdvecDiffOperator(vgrid,jacobian,nu_s_term,nu_par_term,[0.0,0.0],source)

   return collop
end
function collisionOperator(vgrid::Grid, a::SpeciesData, b::Array{SpeciesData,1})
   Nv = length(vgrid.ctr)
   collop = spzeros(Nv,Nv)
   for is in 1:length(b)
      collop = collop + collisionOperator(vgrid,a,b[is])
   end
   return collop
end

"""
solveCollisionalEquation(a::SpeciesData,b::Array{SpeciesData,1};rhs::Array{Float64,1}=[0.0];extraDrag=[0.0];extraDiff=[0.0])

Solves equations of the form C[F] = S with additional terms modifying collision operator accordingly. Subject to constraint that the solution integrates to a.dens.
"""
function solveCollisionalEquation(a::SpeciesData,b::Array{SpeciesData,1};rhs::Array{Float64,1}=[0.0],extraDrag::Array{Float64,1}=[0.0],extraDiff::Array{Float64,1}=[0.0])

end

"""
nuhat*v_ta^3. Useful for non-Maxwellian species, where v_ta is not defined.
"""
function nuhatvta3(a::SpeciesData,b::SpeciesData)
  return b.density * b.charge^2 * a.charge^2 * lnLambda(a,b) / (4.0*pi*ep0^2 * a.mass^2)
end

"""
nuhat from Helander and Sigmar.
"""
function nuhat(a::SpeciesData,b::SpeciesData)
   return nuhatvta3(a,b) / (2.0*a.temperature/a.mass)^1.5
end

"""
lnLambda(a::SpeciesData, b::SpeciesData)

Calculates the Coulomb logarithm ln(Λᵃᵇ ) using approximations in NRL plasma formulary"
"""
function lnLambda(s1::SpeciesData, s2::SpeciesData)
   local me_threshold = 0.05*mp	# Maximum mass that can be considered an "electron", in case one needs reduced mass ratio

   # Check for electron-electron collisions
   if (s1.mass < me_threshold) && (s2.mass < me_threshold)
      # Convert to eV and cm^3
      Te = s1.temperature/el
      ne = s1.density*1.0e-6
      return 23.5 - log(sqrt(ne)/Te^1.25) - sqrt(1.0e-5 + (log(Te)-2.0)^2/16.0)
   # Check for electron-ion collisions
   elseif (s1.mass < me_threshold) || (s2.mass < me_threshold)
      if (s1.mass < me_threshold)
         me = s1.mass; Te = s1.temperature/el; ne = 1.e-6*s1.density
         mi = s2.mass; Ti = s2.temperature/el; ni = 1.e-6*s2.density; Z=s2.charge/el
      else
         mi = s1.mass; Ti = s1.temperature/el; ni = 1.e-6*s1.density; Z=s1.charge/el
         me = s2.mass; Te = s2.temperature/el; ne = 1.e-6*s2.density
      end
      if Te < Ti*Z*me/mi
         return 30.0 - log(sqrt(ni)*Z^2*(mp/mi)/Ti^1.5)
      elseif Te > 10.0*Z^2
         return 24.0 - log(sqrt(ne)/Te)
      else
         return 23.0 - log(sqrt(ne)*Z/Te^1.5)
      end
  # Ion-ion collision
  else
    Z1 = s1.charge/el; Z2 = s2.charge/el
    m1 = s1.mass/mp; m2 = s2.mass/mp
    n1 = s1.density*1.0e-6; n2 = s2.density*1.0e-6
    T1 = s1.temperature/el; T2 = s2.temperature/el
    return 23.0 - log( abs(Z1*Z2)* (m1+m2)*sqrt( (n1*Z1^2/T1) + (n2*Z2^2/T2)) /(m1*T2 + m2*T1 ))
  end
end

"""
slowingDownTime(a::SpeciesData, el::SpeciesData)

Calculates the slowing-down time of species a against electrons where v_ta << v_te.
"""
function slowingdowntime(a::SpeciesData, electron::SpeciesData)
   num = (3.0/(16.0*sqrt(pi)))*electron.mass*a.mass*(2.0*electron.temperature/electron.mass)^1.5
   return num/(a.charge*el^2*electron.density*lnLambda(a,electron))
end

end

max_stars_repo_name: gjwilkie/PlasmaTools

max_issues_repo_name: gjwilkie/PlasmaTools

max_forks_repo_name: gjwilkie/PlasmaTools


text:
module UnderscoreOh

export _o, _nt

# Use README as the docstring of the module:
@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), r"^```julia"m => "```jldoctest README")
end UnderscoreOh

using Base.Broadcast: Broadcasted

# --- Call graph

abstract type Graph <: Function end
struct Hole <: Graph end
struct Call{F,A,K} <: Graph
    f::F
    args::A
    kwargs::K
end
struct GetProperty{name,T} <: Graph
    object::T
end
GetProperty{name}(x::T) where {name,T} = GetProperty{name,T}(x)
call(f, args...; kwargs...) = Call(f, args, kwargs.data)

Base.getproperty(g::Graph, name::Symbol) = GetProperty{name}(g)
Base.getproperty(g::Graph, prop) = call(getproperty, g, prop)
Base.getindex(g::Graph, idx...) = call(getindex, g, idx...)

_f(x) = getfield(x, :f)
_args(x) = getfield(x, :args)
_kwargs(x) = getfield(x, :kwargs)
_object(x) = getfield(x, :object)

# --- Broadcasting

struct UnderscoreOhStyle <: Broadcast.BroadcastStyle end
Base.BroadcastStyle(::Type{<:Graph}) = UnderscoreOhStyle()
Base.BroadcastStyle(::UnderscoreOhStyle, ::Broadcast.BroadcastStyle) = UnderscoreOhStyle()

Base.broadcastable(x::Graph) = x
Broadcast.instantiate(bc::Broadcasted{UnderscoreOhStyle}) = bc

# Used only for `KWBroadcastedInner` below:
function air end
struct Aired{T}
    value::T
end
Broadcast.broadcasted(::typeof(air), x) = Aired(x)
Broadcast.materialize(x::Aired) = x.value

# A hack to support broadcasting with keyword arguments:
const KWBroadcastedInner = let T = typeof(air.(identity.(; dummy = 1)).f)
    getfield(parentmodule(T), nameof(T))
end

# Convert broadcasted to a callable
Base.copy(bc::Broadcasted{UnderscoreOhStyle}) = asgraph(bc)

asgraph(x) = x
asgraph(x::Ref) = x[]
function asgraph(bc::Broadcasted)
    if bc.f isa KWBroadcastedInner
        args = map(asgraph, bc.args)
        kwargs = _map(asgraph, (; bc.f.kwargs...))
        return call(bc.f.f, args...; kwargs...)
    else
        return call(bc.f, map(asgraph, bc.args)...)
    end
end

struct Guarded{F} <: Function
    f::F
end
(g::Guarded)(x) = _f(g)(x)

Base.:~(g::Graph) = Guarded(g)
Base.:~(::Hole, g::Graph) = Guarded(g)

function _nt(; kwargs...)
    if any(x -> x isa Graph, values(kwargs))
        return call(_nt; kwargs...)
    else
        return (; kwargs...)
    end
end

# --- Evaluation

@inline (g::Hole)(x) = materialize(g, x)
@inline (g::Call)(x) = materialize(g, x)
@inline (g::GetProperty)(x) = materialize(g, x)

@inline materialize(y, _) = y
@inline materialize(g::Hole, x) = x
@inline materialize(g::Call, x) = _f(g)(feed(x, _args(g))...; feed(x, _kwargs(g))...)
@inline materialize(g::GetProperty{name}, x) where {name} =
    getproperty(materialize(_object(g), x), name)

@inline feed(x, args::Tuple) = map(fix2(materialize, x), args)
@inline feed(x, kwargs::NamedTuple) = _map(fix2(materialize, x), kwargs)

@inline _map(f, xs::NamedTuple{names}) where {names} = NamedTuple{names}(map(f, Tuple(xs)))

@inline function fix2(f, a2)
    @inline function fixed(a1)
        return f(a1, a2)
    end
    return fixed
end

# --- Indexing

# Base.getindex(x::AbstractArray, g::Graph) = g.(x)

# --- Show it like you build it

maybe_print_dot(io) = get(io, :_secrete_print_dot_key, true) ? print(io, '.') : nothing
unset_print_dot(io) = IOContext(io, :_secrete_print_dot_key => false)

@nospecialize

Base.show(io::IO, ::Hole) = printstyled(io, "_o"; color = :light_black)
Base.show(io::IO, g::Call) = show_impl(io, _f(g), _args(g), _kwargs(g))
Base.show(io::IO, g::GetProperty{name}) where {name} =
    show_impl(io, getproperty, (_object(g), name), NamedTuple())

is_binop(f) = Base.isbinaryoperator(nameof(f))

show_impl(io, f, args, kwargs) = show_impl_default(io, f, args, kwargs)

function show_impl_default(io, f, args, kwargs)
    if length(args) > 0 && isempty(kwargs) && is_binop(f)
        # Print as binary operator:
        show_term(io, f, args[1])
        for a in args[2:end]
            print(io, ' ')
            maybe_print_dot(io)
            print(io, nameof(f), ' ')
            show_term(io, f, a)
        end
        return
    end
    print(io, f)
    maybe_print_dot(io)
    print(io, '(')
    if length(args) > 0
        show(io, args[1])
        for a in args[2:end]
            print(io, ", ")
            show(io, a)
        end
    end
    if length(kwargs) > 0
        print(io, "; ")
        isfirst = true
        for (k, v) in pairs(kwargs)
            if isfirst
                isfirst = false
            else
                print(io, ", ")
            end
            print(io, k, "=")
            show(io, v)
        end
    end
    print(io, ')')
end

function show_term(io, f, a)
    if need_paren(f, a)
        print(io, '(')
        show(io, a)
        print(io, ')')
    else
        show(io, a)
    end
end

need_paren(f, g) = false
need_paren(f, g::Call) = is_binop(_f(g)) && _f(g) !== f

function show_impl(io, ::typeof(getproperty), args, kwargs)
    if length(args) == 2 && length(kwargs) == 0
        value, prop = args
        if prop isa String
            show(io, value)
            print(io, '.')
            show(io, prop)
            return
        end
        tmpstr = string(Expr(:., :_, QuoteNode(prop)))
        if startswith(tmpstr, "_.")
            show(io, value)
            print(io, tmpstr[2:end])
            return
        end
    end
    show_impl_default(io, getproperty, args, kwargs)
    return
end

function show_impl(io, ::typeof(getindex), args, kwargs)
    @assert length(kwargs) == 0
    show(io, args[1])
    print(io, '[')
    if length(args) > 1
        show(io, args[2])
        for a in args[3:end]
            print(io, ", ")
            show(io, a)
        end
    end
    print(io, ']')
    return
end

function Base.show(io::IO, g::Guarded)
    print(io, "~(")
    show(io, _f(g))
    print(io, ")")
end

function Base.show(io::IO, ::MIME"text/plain", g::Union{Graph,Guarded})
    if get(io, :compact, false) === true
        printstyled(io, '('; color = :light_black)
        show(io, g)
        printstyled(io, ')'; color = :light_black)
        return
    end
    print(io, _o, " -> ")
    show(unset_print_dot(io), g)
    n = length(methods(g))
    m = n == 1 ? "method" : "methods"
    print(io, " (generic function with $n $m)")
end

# `Base.Function` defines those methods to be something different
Base.print(io::IO, g::Graph) = show(io, g)
Base.print(io::IO, g::Guarded) = show(io, g)

@specialize

# --- User interface

const _o = Hole()

end # module