text:
##
## Join / merge
##

# Like similar, but returns a array that can have missings and is initialized with missings
similar_missing(dv::AbstractArray{T}, dims::Union{Int, Tuple{Vararg{Int}}}) where {T} =
    fill!(similar(dv, Union{T, Missing}, dims), missing)

similar_outer(leftcol::AbstractVector, rightcol::AbstractVector, n::Int) =
    Tables.allocatecolumn(promote_type(eltype(leftcol), eltype(rightcol)), n)

const OnType = Union{SymbolOrString, NTuple{2, Symbol}, Pair{Symbol, Symbol},
                     Pair{<:AbstractString, <:AbstractString}}

# helper structure for DataFrames joining
struct DataFrameJoiner
    dfl::AbstractDataFrame
    dfr::AbstractDataFrame
    dfl_on::AbstractDataFrame
    dfr_on::AbstractDataFrame
    left_on::Vector{Symbol}
    right_on::Vector{Symbol}

    function DataFrameJoiner(dfl::AbstractDataFrame, dfr::AbstractDataFrame,
                             on::Union{<:OnType, AbstractVector},
                             matchmissing::Symbol)
        on_cols = isa(on, AbstractVector) ? on : [on]
        left_on = Symbol[]
        right_on = Symbol[]
        for v in on_cols
            if v isa SymbolOrString
                push!(left_on, Symbol(v))
                push!(right_on, Symbol(v))
            elseif v isa Union{Pair{Symbol, Symbol},
                               Pair{<:AbstractString, <:AbstractString}}
                push!(left_on, Symbol(first(v)))
                push!(right_on, Symbol(last(v)))
            elseif v isa NTuple{2, Symbol}
                # an explicit error is thrown as Tuple{Symbol, Symbol} was supported in the past
                throw(ArgumentError("Using a `Tuple{Symbol, Symbol}` or a vector containing " *
                                    "such tuples as a value of `on` keyword argument is " *
                                    "not supported: use `Pair{Symbol, Symbol}` instead."))
            else
                throw(ArgumentError("All elements of `on` argument to `join` must be " *
                                    "Symbol or Pair{Symbol, Symbol}."))
            end
        end
        dfl_on = dfl[!, left_on]
        dfr_on = dfr[!, right_on]

        if matchmissing === :error
            for df in (dfl_on, dfr_on), col in eachcol(df)
                if any(ismissing, col)
                    throw(ArgumentError("missing values in key columns are not allowed " *
                                        "when matchmissing == :error"))
                end
            end
        elseif matchmissing !== :equal
            throw(ArgumentError("matchmissing allows only :error or :equal"))
        end

        for df in (dfl_on, dfr_on), col in eachcol(df)
            if any(x -> (x isa Union{Complex, Real}) &&
                        (isnan(x) || isequal(real(x), -0.0) || isequal(imag(x), -0.0)), col)
                throw(ArgumentError("currently for numeric values NaN and `-0.0` " *
                                    "in their real or imaginary components are not " *
                                    "allowed. Use CategoricalArrays.jl to wrap " *
                                    "these values in a CategoricalVector to perform " *
                                    "the requested join."))
            end
        end

        new(dfl, dfr, dfl_on, dfr_on, left_on, right_on)
    end
end

# helper map between the row indices in original and joined table
struct RowIndexMap
    "row indices in the original table"
    orig::Vector{Int}
    "row indices in the resulting joined table"
    join::Vector{Int}
end

Base.length(x::RowIndexMap) = length(x.orig)

_rename_cols(old_names::AbstractVector{Symbol},
             renamecols::Union{Function, Symbol, AbstractString},
             exclude::AbstractVector{Symbol} = Symbol[]) =
    Symbol[n in exclude ? n :
           (renamecols isa Function ? Symbol(renamecols(string(n))) : Symbol(n, renamecols))
           for n in old_names]

function compose_inner_table(joiner::DataFrameJoiner,
                             makeunique::Bool,
                             left_rename::Union{Function, AbstractString, Symbol},
                             right_rename::Union{Function, AbstractString, Symbol})
    left_ixs, right_ixs = find_inner_rows(joiner)
    dfl = joiner.dfl[left_ixs, :]
    dfr_noon = joiner.dfr[right_ixs, Not(joiner.right_on)]

    ncleft = ncol(dfl)
    cols = Vector{AbstractVector}(undef, ncleft + ncol(dfr_noon))

    for (i, col) in enumerate(eachcol(dfl))
        cols[i] = col
    end
    for (i, col) in enumerate(eachcol(dfr_noon))
        cols[i+ncleft] = col
    end

    new_names = vcat(_rename_cols(_names(joiner.dfl), left_rename, joiner.left_on),
                     _rename_cols(_names(dfr_noon), right_rename))
    res = DataFrame(cols, new_names, makeunique=makeunique, copycols=false)

    return res
end

function find_missing_idxs(present::Vector{Int}, target_len::Int)
    not_seen = trues(target_len)
    @inbounds for v in present
        not_seen[v] = false
    end
    return findall(not_seen)
end

function compose_joined_table(joiner::DataFrameJoiner, kind::Symbol, makeunique::Bool,
                              left_rename::Union{Function, AbstractString, Symbol},
                              right_rename::Union{Function, AbstractString, Symbol},
                              indicator::Union{Nothing, Symbol, AbstractString})
    @assert kind == :left || kind == :right || kind == :outer
    left_ixs, right_ixs = find_inner_rows(joiner)

    if kind == :left || kind == :outer
        leftonly_ixs = find_missing_idxs(left_ixs, nrow(joiner.dfl))
    else
        leftonly_ixs = Int[]
    end

    if kind == :right || kind == :outer
        rightonly_ixs = find_missing_idxs(right_ixs, nrow(joiner.dfr))
    else
        rightonly_ixs = Int[]
    end

    lil = length(left_ixs)
    ril = length(right_ixs)
    loil = length(leftonly_ixs)
    roil = length(rightonly_ixs)

    @assert lil == ril

    dfl_noon = select(joiner.dfl, Not(joiner.left_on), copycols=false)
    dfr_noon = select(joiner.dfr, Not(joiner.right_on), copycols=false)

    target_nrow = lil + loil + roil

    _similar_left = kind == :left ? similar : similar_missing
    _similar_right = kind == :right ? similar : similar_missing

    if isnothing(indicator)
        src_indicator = nothing
    else
        src_indicator = Vector{UInt32}(undef, target_nrow)
        src_indicator[1:lil] .= 3
        src_indicator[lil + 1:lil + loil] .= 1
        src_indicator[lil + loil + 1:target_nrow] .= 2
    end

    cols = Vector{AbstractVector}(undef, ncol(joiner.dfl) + ncol(dfr_noon))

    col_idx = 1

    left_idxs = [columnindex(joiner.dfl, n) for n in joiner.left_on]
    append!(left_idxs, setdiff(1:ncol(joiner.dfl), left_idxs))

    if kind == :left
        @assert loil + lil == target_nrow
        for col in eachcol(joiner.dfl_on)
            cols_i = left_idxs[col_idx]
            cols[cols_i] = similar(col, target_nrow)
            copyto!(cols[cols_i], view(col, left_ixs))
            copyto!(cols[cols_i], lil + 1, view(col, leftonly_ixs), 1, loil)
            col_idx += 1
        end
    elseif kind == :right
        @assert roil + ril == target_nrow
        for col in eachcol(joiner.dfr_on)
            cols_i = left_idxs[col_idx]
            cols[cols_i] = similar(col, target_nrow)
            copyto!(cols[cols_i], view(col, right_ixs))
            copyto!(cols[cols_i], lil + 1, view(col, rightonly_ixs), 1, roil)
            col_idx += 1
        end
    else
        @assert kind == :outer
        @assert loil + roil + ril == target_nrow
        for (lcol, rcol) in zip(eachcol(joiner.dfl_on), eachcol(joiner.dfr_on))
            cols_i = left_idxs[col_idx]
            cols[cols_i] = similar_outer(lcol, rcol, target_nrow)
            copyto!(cols[cols_i], view(lcol, left_ixs))
            copyto!(cols[cols_i], lil + 1, view(lcol, leftonly_ixs), 1, loil)
            copyto!(cols[cols_i], lil + loil + 1, view(rcol, rightonly_ixs), 1, roil)
            col_idx += 1
        end
    end

    @assert col_idx == ncol(joiner.dfl_on) + 1

    for col in eachcol(dfl_noon)
        cols_i = left_idxs[col_idx]
        cols[cols_i] = _similar_left(col, target_nrow)
        copyto!(cols[cols_i], view(col, left_ixs))
        copyto!(cols[cols_i], lil + 1, view(col, leftonly_ixs), 1, loil)
        col_idx += 1
    end

    @assert col_idx == ncol(joiner.dfl) + 1

    for col in eachcol(dfr_noon)
        cols[col_idx] = _similar_right(col, target_nrow)
        copyto!(cols[col_idx], view(col, right_ixs))
        copyto!(cols[col_idx], lil + loil + 1, view(col, rightonly_ixs), 1, roil)
        col_idx += 1
    end

    @assert col_idx == length(cols) + 1

    new_names = vcat(_rename_cols(_names(joiner.dfl), left_rename, joiner.left_on),
                     _rename_cols(_names(dfr_noon), right_rename))
    res = DataFrame(cols, new_names, makeunique=makeunique, copycols=false)

    return res, src_indicator
end

function _join(df1::AbstractDataFrame, df2::AbstractDataFrame;
               on::Union{<:OnType, AbstractVector}, kind::Symbol, makeunique::Bool,
               indicator::Union{Nothing, Symbol, AbstractString},
               validate::Union{Pair{Bool, Bool}, Tuple{Bool, Bool}},
               left_rename::Union{Function, AbstractString, Symbol},
               right_rename::Union{Function, AbstractString, Symbol},
               matchmissing::Symbol)
    _check_consistency(df1)
    _check_consistency(df2)

    if on == []
        throw(ArgumentError("Missing join argument 'on'."))
    end

    joiner = DataFrameJoiner(df1, df2, on, matchmissing)

    # Check merge key validity
    left_invalid = validate[1] ? any(nonunique(joiner.dfl, joiner.left_on)) : false
    right_invalid = validate[2] ? any(nonunique(joiner.dfr, joiner.right_on)) : false

    if validate[1]
        non_unique_left = nonunique(joiner.dfl, joiner.left_on)
        if any(non_unique_left)
            left_invalid = true
            non_unique_left_df = unique(joiner.dfl[non_unique_left, joiner.left_on])
            nrow_nul_df = nrow(non_unique_left_df)
            @assert nrow_nul_df > 0
            if nrow_nul_df == 1
                left_invalid_msg = "df1 contains 1 duplicate key: " *
                                   "$(NamedTuple(non_unique_left_df[1, :])). "
            elseif nrow_nul_df == 2
                left_invalid_msg = "df1 contains 2 duplicate keys: " *
                                   "$(NamedTuple(non_unique_left_df[1, :])) and " *
                                   "$(NamedTuple(non_unique_left_df[2, :])). "
            else
                left_invalid_msg = "df1 contains $(nrow_nul_df) duplicate keys: " *
                                   "$(NamedTuple(non_unique_left_df[1, :])), ..., " *
                                   "$(NamedTuple(non_unique_left_df[end, :])). "
            end

        else
            left_invalid = false
            left_invalid_msg = ""
        end
    else
        left_invalid = false
        left_invalid_msg = ""
    end

    if validate[2]
        non_unique_right = nonunique(joiner.dfr, joiner.right_on)
        if any(non_unique_right)
            right_invalid = true
            non_unique_right_df = unique(joiner.dfr[non_unique_right, joiner.right_on])
            nrow_nur_df = nrow(non_unique_right_df)
            @assert nrow_nur_df > 0
            if nrow_nur_df == 1
                right_invalid_msg = "df2 contains 1 duplicate key: " *
                                    "$(NamedTuple(non_unique_right_df[1, :]))."
            elseif nrow_nur_df == 2
                right_invalid_msg = "df2 contains 2 duplicate keys: " *
                                    "$(NamedTuple(non_unique_right_df[1, :])) and " *
                                    "$(NamedTuple(non_unique_right_df[2, :]))."
            else
                right_invalid_msg = "df2 contains $(nrow_nur_df) duplicate keys: " *
                                    "$(NamedTuple(non_unique_right_df[1, :])), ..., " *
                                    "$(NamedTuple(non_unique_right_df[end, :]))."
            end

        else
            right_invalid = false
            right_invalid_msg = ""
        end
    else
        right_invalid = false
        right_invalid_msg = ""
    end

    if left_invalid && right_invalid
        first_error_df1 = findfirst(nonunique(joiner.dfl, joiner.left_on))
        first_error_df2 = findfirst(nonunique(joiner.dfr, joiner.right_on))
        throw(ArgumentError("Merge key(s) are not unique in both df1 and df2. " *
                            left_invalid_msg * right_invalid_msg))
    elseif left_invalid
        first_error = findfirst(nonunique(joiner.dfl, joiner.left_on))
        throw(ArgumentError("Merge key(s) in df1 are not unique. " *
                            left_invalid_msg))
    elseif right_invalid
        first_error = findfirst(nonunique(joiner.dfr, joiner.right_on))
        throw(ArgumentError("Merge key(s) in df2 are not unique. " *
                            right_invalid_msg))
    end

    src_indicator = nothing
    if kind == :inner
        joined = compose_inner_table(joiner, makeunique, left_rename, right_rename)
    elseif kind == :left
        joined, src_indicator =
            compose_joined_table(joiner, kind, makeunique, left_rename, right_rename, indicator)
    elseif kind == :right
        joined, src_indicator =
            compose_joined_table(joiner, kind, makeunique, left_rename, right_rename, indicator)
    elseif kind == :outer
        joined, src_indicator =
            compose_joined_table(joiner, kind, makeunique, left_rename, right_rename, indicator)
    elseif kind == :semi
        joined = joiner.dfl[find_semi_rows(joiner), :]
    elseif kind == :anti
        joined = joiner.dfl[.!find_semi_rows(joiner), :]
    else
        throw(ArgumentError("Unknown kind of join requested: $kind"))
    end

    if indicator !== nothing
        pool = ["left_only", "right_only", "both"]
        invpool = Dict{String, UInt32}("left_only" => 1,
                                       "right_only" => 2,
                                       "both" => 3)
        indicatorcol = PooledArray(PooledArrays.RefArray(src_indicator),
                                   invpool, pool)

        unique_indicator = indicator
        if makeunique
            try_idx = 0
            while hasproperty(joined, unique_indicator)
                try_idx += 1
                unique_indicator = Symbol(string(indicator, "_", try_idx))
            end
        end

        if hasproperty(joined, unique_indicator)
            throw(ArgumentError("joined data frame already has column " *
                                ":$unique_indicator. Pass makeunique=true to " *
                                "make it unique using a suffix automatically."))
        end
        joined[!, unique_indicator] = indicatorcol
    else
        @assert isnothing(src_indicator)
    end

    return joined
end

"""
    innerjoin(df1, df2; on, makeunique=false, validate=(false, false),
              renamecols=(identity => identity), matchmissing=:error)
    innerjoin(df1, df2, dfs...; on, makeunique=false,
              validate=(false, false), matchmissing=:error)

Perform an inner join of two or more data frame objects and return a `DataFrame`
containing the result. An inner join includes rows with keys that match in all
passed data frames.

The order of rows in the result is undefined and may change in the future releases.

In the returned data frame the type of the columns on which the data frames are
joined is determined by the type of these columns in `df1`. This behavior may
change in future releases.

# Arguments
- `df1`, `df2`, `dfs...`: the `AbstractDataFrames` to be joined

# Keyword Arguments
- `on` : A column name to join `df1` and `df2` on. If the columns on which
  `df1` and `df2` will be joined have different names, then a `left=>right`
  pair can be passed. It is also allowed to perform a join on multiple columns,
  in which case a vector of column names or column name pairs can be passed
  (mixing names and pairs is allowed). If more than two data frames are joined
  then only a column name or a vector of column names are allowed.
  `on` is a required argument.
- `makeunique` : if `false` (the default), an error will be raised
  if duplicate names are found in columns not joined on;
  if `true`, duplicate names will be suffixed with `_i`
  (`i` starting at 1 for the first duplicate).
- `validate` : whether to check that columns passed as the `on` argument
  define unique keys in each input data frame (according to `isequal`).
  Can be a tuple or a pair, with the first element indicating whether to
  run check for `df1` and the second element for `df2`.
  By default no check is performed.
- `renamecols` : a `Pair` specifying how columns of left and right data frames should
  be renamed in the resulting data frame. Each element of the pair can be a
  string or a `Symbol` can be passed in which case it is appended to the original
  column name; alternatively a function can be passed in which case it is applied
  to each column name, which is passed to it as a `String`. Note that `renamecols`
  does not affect `on` columns, whose names are always taken from the left
  data frame and left unchanged.
- `matchmissing` : if equal to `:error` throw an error if `missing` is present
  in `on` columns; if equal to `:equal` then `missing` is allowed and missings are
  matched (`isequal` is used for comparisons of rows for equality)

It is not allowed to join on columns that contain `NaN` or `-0.0` in real or
imaginary part of the number. If you need to perform a join on such values use
CategoricalArrays.jl and transform a column containing such values into a
`CategoricalVector`.

When merging `on` categorical columns that differ in the ordering of their
levels, the ordering of the left data frame takes precedence over the ordering
of the right data frame.

If more than two data frames are passed, the join is performed recursively with
left associativity. In this case the `validate` keyword argument is applied
recursively with left associativity.

See also: [`leftjoin`](@ref), [`rightjoin`](@ref), [`outerjoin`](@ref),
          [`semijoin`](@ref), [`antijoin`](@ref), [`crossjoin`](@ref).

# Examples
```julia
julia> name = DataFrame(ID = [1, 2, 3], Name = ["John Doe", "Jane Doe", "Joe Blogs"])
3×2 DataFrame
 Row │ ID     Name
     │ Int64  String
─────┼──────────────────
   1 │     1  John Doe
   2 │     2  Jane Doe
   3 │     3  Joe Blogs

julia> job = DataFrame(ID = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])
3×2 DataFrame
 Row │ ID     Job
     │ Int64  String
─────┼───────────────
   1 │     1  Lawyer
   2 │     2  Doctor
   3 │     4  Farmer

julia> innerjoin(name, job, on = :ID)
2×3 DataFrame
 Row │ ID     Name      Job
     │ Int64  String    String
─────┼─────────────────────────
   1 │     1  John Doe  Lawyer
   2 │     2  Jane Doe  Doctor

julia> job2 = DataFrame(identifier = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])
3×2 DataFrame
 Row │ identifier  Job
     │ Int64       String
─────┼────────────────────
   1 │          1  Lawyer
   2 │          2  Doctor
   3 │          4  Farmer

julia> innerjoin(name, job2, on = :ID => :identifier, renamecols = "_left" => "_right")
2×3 DataFrame
 Row │ ID     Name_left  Job_right
     │ Int64  String     String
─────┼─────────────────────────────
   1 │     1  John Doe   Lawyer
   2 │     2  Jane Doe   Doctor

julia> innerjoin(name, job2, on = [:ID => :identifier], renamecols = uppercase => lowercase)
2×3 DataFrame
 Row │ ID     NAME      job
     │ Int64  String    String
─────┼─────────────────────────
   1 │     1  John Doe  Lawyer
   2 │     2  Jane Doe  Doctor
```
"""
function innerjoin(df1::AbstractDataFrame, df2::AbstractDataFrame;
                   on::Union{<:OnType, AbstractVector} = Symbol[],
                   makeunique::Bool=false,
                   validate::Union{Pair{Bool, Bool}, Tuple{Bool, Bool}}=(false, false),
                   renamecols::Pair=identity => identity,
                   matchmissing::Symbol=:error)
    if !all(x -> x isa Union{Function, AbstractString, Symbol}, renamecols)
        throw(ArgumentError("renamecols keyword argument must be a `Pair` " *
                            "containing functions, strings, or `Symbol`s"))
    end
    return _join(df1, df2, on=on, kind=:inner, makeunique=makeunique,
                 indicator=nothing, validate=validate,
                 left_rename=first(renamecols), right_rename=last(renamecols),
                 matchmissing=matchmissing)
end

innerjoin(df1::AbstractDataFrame, df2::AbstractDataFrame, dfs::AbstractDataFrame...;
          on::Union{<:OnType, AbstractVector} = Symbol[],
          makeunique::Bool=false,
          validate::Union{Pair{Bool, Bool}, Tuple{Bool, Bool}}=(false, false),
          matchmissing::Symbol=:error) =
    innerjoin(innerjoin(df1, df2, on=on, makeunique=makeunique, validate=validate),
              dfs..., on=on, makeunique=makeunique, validate=validate,
              matchmissing=matchmissing)

"""
    leftjoin(df1, df2; on, makeunique=false, indicator=nothing, validate=(false, false),
             renamecols=(identity => identity), matchmissing=:error)

Perform a left join of twodata frame objects and return a `DataFrame` containing
the result. A left join includes all rows from `df1`.

The order of rows in the result is undefined and may change in the future releases.

In the returned data frame the type of the columns on which the data frames are
joined is determined by the type of these columns in `df1`. This behavior may
change in future releases.

# Arguments
- `df1`, `df2`: the `AbstractDataFrames` to be joined

# Keyword Arguments
- `on` : A column name to join `df1` and `df2` on. If the columns on which
  `df1` and `df2` will be joined have different names, then a `left=>right`
  pair can be passed. It is also allowed to perform a join on multiple columns,
  in which case a vector of column names or column name pairs can be passed
  (mixing names and pairs is allowed).
- `makeunique` : if `false` (the default), an error will be raised
  if duplicate names are found in columns not joined on;
  if `true`, duplicate names will be suffixed with `_i`
  (`i` starting at 1 for the first duplicate).
- `indicator` : Default: `nothing`. If a `Symbol` or string, adds categorical indicator
  column with the given name, for whether a row appeared in only `df1` (`"left_only"`),
  only `df2` (`"right_only"`) or in both (`"both"`). If the name is already in use,
  the column name will be modified if `makeunique=true`.
- `validate` : whether to check that columns passed as the `on` argument
  define unique keys in each input data frame (according to `isequal`).
  Can be a tuple or a pair, with the first element indicating whether to
  run check for `df1` and the second element for `df2`.
  By default no check is performed.
- `renamecols` : a `Pair` specifying how columns of left and right data frames should
  be renamed in the resulting data frame. Each element of the pair can be a
  string or a `Symbol` can be passed in which case it is appended to the original
  column name; alternatively a function can be passed in which case it is applied
  to each column name, which is passed to it as a `String`. Note that `renamecols`
  does not affect `on` columns, whose names are always taken from the left
  data frame and left unchanged.
- `matchmissing` : if equal to `:error` throw an error if `missing` is present
  in `on` columns; if equal to `:equal` then `missing` is allowed and missings are
  matched (`isequal` is used for comparisons of rows for equality)

All columns of the returned data table will support missing values.

It is not allowed to join on columns that contain `NaN` or `-0.0` in real or
imaginary part of the number. If you need to perform a join on such values use
CategoricalArrays.jl and transform a column containing such values into a
`CategoricalVector`.

When merging `on` categorical columns that differ in the ordering of their
levels, the ordering of the left data frame takes precedence over the ordering
of the right data frame.

See also: [`innerjoin`](@ref), [`rightjoin`](@ref), [`outerjoin`](@ref),
          [`semijoin`](@ref), [`antijoin`](@ref), [`crossjoin`](@ref).

# Examples
```julia
julia> name = DataFrame(ID = [1, 2, 3], Name = ["John Doe", "Jane Doe", "Joe Blogs"])
3×2 DataFrame
 Row │ ID     Name
     │ Int64  String
─────┼──────────────────
   1 │     1  John Doe
   2 │     2  Jane Doe
   3 │     3  Joe Blogs

julia> job = DataFrame(ID = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])
3×2 DataFrame
 Row │ ID     Job
     │ Int64  String
─────┼───────────────
   1 │     1  Lawyer
   2 │     2  Doctor
   3 │     4  Farmer

julia> leftjoin(name, job, on = :ID)
3×3 DataFrame
 Row │ ID     Name       Job
     │ Int64  String     String?
─────┼───────────────────────────
   1 │     1  John Doe   Lawyer
   2 │     2  Jane Doe   Doctor
   3 │     3  Joe Blogs  missing

julia> job2 = DataFrame(identifier = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])
3×2 DataFrame
 Row │ identifier  Job
     │ Int64       String
─────┼────────────────────
   1 │          1  Lawyer
   2 │          2  Doctor
   3 │          4  Farmer

julia> leftjoin(name, job2, on = :ID => :identifier, renamecols = "_left" => "_right")
3×3 DataFrame
 Row │ ID     Name_left  Job_right
     │ Int64  String     String?
─────┼─────────────────────────────
   1 │     1  John Doe   Lawyer
   2 │     2  Jane Doe   Doctor
   3 │     3  Joe Blogs  missing

julia> leftjoin(name, job2, on = [:ID => :identifier], renamecols = uppercase => lowercase)
3×3 DataFrame
 Row │ ID     NAME       job
     │ Int64  String     String?
─────┼───────────────────────────
   1 │     1  John Doe   Lawyer
   2 │     2  Jane Doe   Doctor
   3 │     3  Joe Blogs  missing
```
"""
function leftjoin(df1::AbstractDataFrame, df2::AbstractDataFrame;
         on::Union{<:OnType, AbstractVector} = Symbol[],
         makeunique::Bool=false, indicator::Union{Nothing, Symbol, AbstractString} = nothing,
         validate::Union{Pair{Bool, Bool}, Tuple{Bool, Bool}}=(false, false),
         renamecols::Pair=identity => identity, matchmissing::Symbol=:error)
    if !all(x -> x isa Union{Function, AbstractString, Symbol}, renamecols)
        throw(ArgumentError("renamecols keyword argument must be a `Pair` " *
                            "containing functions, strings, or `Symbol`s"))
    end
    return _join(df1, df2, on=on, kind=:left, makeunique=makeunique,
                 indicator=indicator, validate=validate,
                 left_rename=first(renamecols), right_rename=last(renamecols),
                 matchmissing=matchmissing)
end

"""
    rightjoin(df1, df2; on, makeunique=false, indicator = nothing,
              validate=(false, false), renamecols=(identity => identity),
              matchmissing=:error)

Perform a right join on two data frame objects and return a `DataFrame` containing
the result. A right join includes all rows from `df2`.

The order of rows in the result is undefined and may change in the future releases.

In the returned data frame the type of the columns on which the data frames are
joined is determined by the type of these columns in `df2`. This behavior may
change in future releases.

# Arguments
- `df1`, `df2`: the `AbstractDataFrames` to be joined

# Keyword Arguments
- `on` : A column name to join `df1` and `df2` on. If the columns on which
  `df1` and `df2` will be joined have different names, then a `left=>right`
  pair can be passed. It is also allowed to perform a join on multiple columns,
  in which case a vector of column names or column name pairs can be passed
  (mixing names and pairs is allowed).
- `makeunique` : if `false` (the default), an error will be raised
  if duplicate names are found in columns not joined on;
  if `true`, duplicate names will be suffixed with `_i`
  (`i` starting at 1 for the first duplicate).
- `indicator` : Default: `nothing`. If a `Symbol` or string, adds categorical indicator
  column with the given name for whether a row appeared in only `df1` (`"left_only"`),
  only `df2` (`"right_only"`) or in both (`"both"`). If the name is already in use,
  the column name will be modified if `makeunique=true`.
- `validate` : whether to check that columns passed as the `on` argument
  define unique keys in each input data frame (according to `isequal`).
  Can be a tuple or a pair, with the first element indicating whether to
  run check for `df1` and the second element for `df2`.
  By default no check is performed.
- `renamecols` : a `Pair` specifying how columns of left and right data frames should
  be renamed in the resulting data frame. Each element of the pair can be a
  string or a `Symbol` can be passed in which case it is appended to the original
  column name; alternatively a function can be passed in which case it is applied
  to each column name, which is passed to it as a `String`. Note that `renamecols`
  does not affect `on` columns, whose names are always taken from the left
  data frame and left unchanged.
- `matchmissing` : if equal to `:error` throw an error if `missing` is present
  in `on` columns; if equal to `:equal` then `missing` is allowed and missings are
  matched (`isequal` is used for comparisons of rows for equality)

All columns of the returned data table will support missing values.

It is not allowed to join on columns that contain `NaN` or `-0.0` in real or
imaginary part of the number. If you need to perform a join on such values use
CategoricalArrays.jl and transform a column containing such values into a
`CategoricalVector`.

When merging `on` categorical columns that differ in the ordering of their
levels, the ordering of the left data frame takes precedence over the ordering
of the right data frame.

See also: [`innerjoin`](@ref), [`leftjoin`](@ref), [`outerjoin`](@ref),
          [`semijoin`](@ref), [`antijoin`](@ref), [`crossjoin`](@ref).

# Examples
```julia
julia> name = DataFrame(ID = [1, 2, 3], Name = ["John Doe", "Jane Doe", "Joe Blogs"])
3×2 DataFrame
 Row │ ID     Name
     │ Int64  String
─────┼──────────────────
   1 │     1  John Doe
   2 │     2  Jane Doe
   3 │     3  Joe Blogs

julia> job = DataFrame(ID = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])
3×2 DataFrame
 Row │ ID     Job
     │ Int64  String
─────┼───────────────
   1 │     1  Lawyer
   2 │     2  Doctor
   3 │     4  Farmer

julia> rightjoin(name, job, on = :ID)
3×3 DataFrame
 Row │ ID     Name      Job
     │ Int64  String?   String
─────┼─────────────────────────
   1 │     1  John Doe  Lawyer
   2 │     2  Jane Doe  Doctor
   3 │     4  missing   Farmer

julia> job2 = DataFrame(identifier = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])
3×2 DataFrame
 Row │ identifier  Job
     │ Int64       String
─────┼────────────────────
   1 │          1  Lawyer
   2 │          2  Doctor
   3 │          4  Farmer

julia> rightjoin(name, job2, on = :ID => :identifier, renamecols = "_left" => "_right")
3×3 DataFrame
 Row │ ID     Name_left  Job_right
     │ Int64  String?    String
─────┼─────────────────────────────
   1 │     1  John Doe   Lawyer
   2 │     2  Jane Doe   Doctor
   3 │     4  missing    Farmer

julia> rightjoin(name, job2, on = [:ID => :identifier], renamecols = uppercase => lowercase)
3×3 DataFrame
 Row │ ID     NAME      job
     │ Int64  String?   String
─────┼─────────────────────────
   1 │     1  John Doe  Lawyer
   2 │     2  Jane Doe  Doctor
   3 │     4  missing   Farmer
```
"""
function rightjoin(df1::AbstractDataFrame, df2::AbstractDataFrame;
          on::Union{<:OnType, AbstractVector} = Symbol[], makeunique::Bool=false,
          indicator::Union{Nothing, Symbol, AbstractString} = nothing,
          validate::Union{Pair{Bool, Bool}, Tuple{Bool, Bool}}=(false, false),
          renamecols::Pair=identity => identity, matchmissing::Symbol=:error)
    if !all(x -> x isa Union{Function, AbstractString, Symbol}, renamecols)
        throw(ArgumentError("renamecols keyword argument must be a `Pair` " *
                            "containing functions, strings, or `Symbol`s"))
    end
    return _join(df1, df2, on=on, kind=:right, makeunique=makeunique,
                 indicator=indicator, validate=validate,
                 left_rename=first(renamecols), right_rename=last(renamecols),
                 matchmissing=matchmissing)
end

"""
    outerjoin(df1, df2; on, makeunique=false, indicator=nothing, validate=(false, false),
              renamecols=(identity => identity), matchmissing=:error)
    outerjoin(df1, df2, dfs...; on, makeunique = false,
              validate = (false, false), matchmissing=:error)

Perform an outer join of two or more data frame objects and return a `DataFrame`
containing the result. An outer join includes rows with keys that appear in any
of the passed data frames.

The order of rows in the result is undefined and may change in the future releases.

In the returned data frame the type of the columns on which the data frames are
joined is determined by the element type of these columns both `df1` and `df2`.
This behavior may change in future releases.

# Arguments
- `df1`, `df2`, `dfs...` : the `AbstractDataFrames` to be joined

# Keyword Arguments
- `on` : A column name to join `df1` and `df2` on. If the columns on which
  `df1` and `df2` will be joined have different names, then a `left=>right`
  pair can be passed. It is also allowed to perform a join on multiple columns,
  in which case a vector of column names or column name pairs can be passed
  (mixing names and pairs is allowed). If more than two data frames are joined
  then only a column name or a vector of column names are allowed.
  `on` is a required argument.
- `makeunique` : if `false` (the default), an error will be raised
  if duplicate names are found in columns not joined on;
  if `true`, duplicate names will be suffixed with `_i`
  (`i` starting at 1 for the first duplicate).
- `indicator` : Default: `nothing`. If a `Symbol` or string, adds categorical indicator
  column with the given name for whether a row appeared in only `df1` (`"left_only"`),
  only `df2` (`"right_only"`) or in both (`"both"`). If the name is already in use,
  the column name will be modified if `makeunique=true`.
  This argument is only supported when joining exactly two data frames.
- `validate` : whether to check that columns passed as the `on` argument
  define unique keys in each input data frame (according to `isequal`).
  Can be a tuple or a pair, with the first element indicating whether to
  run check for `df1` and the second element for `df2`.
  By default no check is performed.
- `renamecols` : a `Pair` specifying how columns of left and right data frames should
  be renamed in the resulting data frame. Each element of the pair can be a
  string or a `Symbol` can be passed in which case it is appended to the original
  column name; alternatively a function can be passed in which case it is applied
  to each column name, which is passed to it as a `String`. Note that `renamecols`
  does not affect `on` columns, whose names are always taken from the left
  data frame and left unchanged.
- `matchmissing` : if equal to `:error` throw an error if `missing` is present
  in `on` columns; if equal to `:equal` then `missing` is allowed and missings are
  matched (`isequal` is used for comparisons of rows for equality)

All columns of the returned data table will support missing values.

It is not allowed to join on columns that contain `NaN` or `-0.0` in real or
imaginary part of the number. If you need to perform a join on such values use
CategoricalArrays.jl and transform a column containing such values into a
`CategoricalVector`.

When merging `on` categorical columns that differ in the ordering of their
levels, the ordering of the left data frame takes precedence over the ordering
of the right data frame.

If more than two data frames are passed, the join is performed
recursively with left associativity.
In this case the `indicator` keyword argument is not supported
and `validate` keyword argument is applied recursively with left associativity.

See also: [`innerjoin`](@ref), [`leftjoin`](@ref), [`rightjoin`](@ref),
          [`semijoin`](@ref), [`antijoin`](@ref), [`crossjoin`](@ref).

# Examples
```julia
julia> name = DataFrame(ID = [1, 2, 3], Name = ["John Doe", "Jane Doe", "Joe Blogs"])
3×2 DataFrame
 Row │ ID     Name
     │ Int64  String
─────┼──────────────────
   1 │     1  John Doe
   2 │     2  Jane Doe
   3 │     3  Joe Blogs

julia> job = DataFrame(ID = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])
3×2 DataFrame
 Row │ ID     Job
     │ Int64  String
─────┼───────────────
   1 │     1  Lawyer
   2 │     2  Doctor
   3 │     4  Farmer

julia> outerjoin(name, job, on = :ID)
4×3 DataFrame
 Row │ ID     Name       Job
     │ Int64  String?    String?
─────┼───────────────────────────
   1 │     1  John Doe   Lawyer
   2 │     2  Jane Doe   Doctor
   3 │     3  Joe Blogs  missing
   4 │     4  missing    Farmer

julia> job2 = DataFrame(identifier = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])
3×2 DataFrame
 Row │ identifier  Job
     │ Int64       String
─────┼────────────────────
   1 │          1  Lawyer
   2 │          2  Doctor
   3 │          4  Farmer

julia> rightjoin(name, job2, on = :ID => :identifier, renamecols = "_left" => "_right")
3×3 DataFrame
 Row │ ID     Name_left  Job_right
     │ Int64  String?    String
─────┼─────────────────────────────
   1 │     1  John Doe   Lawyer
   2 │     2  Jane Doe   Doctor
   3 │     4  missing    Farmer

julia> rightjoin(name, job2, on = [:ID => :identifier], renamecols = uppercase => lowercase)
3×3 DataFrame
 Row │ ID     NAME      job
     │ Int64  String?   String
─────┼─────────────────────────
   1 │     1  John Doe  Lawyer
   2 │     2  Jane Doe  Doctor
   3 │     4  missing   Farmer
```
"""
function outerjoin(df1::AbstractDataFrame, df2::AbstractDataFrame;
          on::Union{<:OnType, AbstractVector} = Symbol[],
          makeunique::Bool=false, indicator::Union{Nothing, Symbol, AbstractString} = nothing,
          validate::Union{Pair{Bool, Bool}, Tuple{Bool, Bool}}=(false, false),
          renamecols::Pair=identity => identity, matchmissing::Symbol=:error)
    if !all(x -> x isa Union{Function, AbstractString, Symbol}, renamecols)
        throw(ArgumentError("renamecols keyword argument must be a `Pair` " *
                            "containing functions, strings, or `Symbol`s"))
    end
    return _join(df1, df2, on=on, kind=:outer, makeunique=makeunique,
                 indicator=indicator, validate=validate,
                 left_rename=first(renamecols), right_rename=last(renamecols),
                 matchmissing=matchmissing)
end

outerjoin(df1::AbstractDataFrame, df2::AbstractDataFrame, dfs::AbstractDataFrame...;
          on::Union{<:OnType, AbstractVector} = Symbol[], makeunique::Bool=false,
          validate::Union{Pair{Bool, Bool}, Tuple{Bool, Bool}}=(false, false),
          matchmissing::Symbol=:error) =
    outerjoin(outerjoin(df1, df2, on=on, makeunique=makeunique, validate=validate),
              dfs..., on=on, makeunique=makeunique, validate=validate,
              matchmissing=matchmissing)

"""
    semijoin(df1, df2; on, makeunique=false, validate=(false, false), matchmissing=:error)

Perform a semi join of two data frame objects and return a `DataFrame`
containing the result. A semi join returns the subset of rows of `df1` that
match with the keys in `df2`.

The order of rows in the result is undefined and may change in the future releases.

# Arguments
- `df1`, `df2`: the `AbstractDataFrames` to be joined

# Keyword Arguments
- `on` : A column name to join `df1` and `df2` on. If the columns on which
  `df1` and `df2` will be joined have different names, then a `left=>right`
  pair can be passed. It is also allowed to perform a join on multiple columns,
  in which case a vector of column names or column name pairs can be passed
  (mixing names and pairs is allowed).
- `makeunique` : if `false` (the default), an error will be raised
  if duplicate names are found in columns not joined on;
  if `true`, duplicate names will be suffixed with `_i`
  (`i` starting at 1 for the first duplicate).
- `indicator` : Default: `nothing`. If a `Symbol` or string, adds categorical indicator
   column with the given name for whether a row appeared in only `df1` (`"left_only"`),
   only `df2` (`"right_only"`) or in both (`"both"`). If the name is already in use,
   the column name will be modified if `makeunique=true`.
- `validate` : whether to check that columns passed as the `on` argument
   define unique keys in each input data frame (according to `isequal`).
   Can be a tuple or a pair, with the first element indicating whether to
   run check for `df1` and the second element for `df2`.
   By default no check is performed.
- `matchmissing` : if equal to `:error` throw an error if `missing` is present
  in `on` columns; if equal to `:equal` then `missing` is allowed and missings are
  matched (`isequal` is used for comparisons of rows for equality)

It is not allowed to join on columns that contain `NaN` or `-0.0` in real or
imaginary part of the number. If you need to perform a join on such values use
CategoricalArrays.jl and transform a column containing such values into a
`CategoricalVector`.

When merging `on` categorical columns that differ in the ordering of their
levels, the ordering of the left data frame takes precedence over the ordering
of the right data frame.

See also: [`innerjoin`](@ref), [`leftjoin`](@ref), [`rightjoin`](@ref),
          [`outerjoin`](@ref), [`antijoin`](@ref), [`crossjoin`](@ref).

# Examples
```julia
julia> name = DataFrame(ID = [1, 2, 3], Name = ["John Doe", "Jane Doe", "Joe Blogs"])
3×2 DataFrame
 Row │ ID     Name
     │ Int64  String
─────┼──────────────────
   1 │     1  John Doe
   2 │     2  Jane Doe
   3 │     3  Joe Blogs

julia> job = DataFrame(ID = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])
3×2 DataFrame
 Row │ ID     Job
     │ Int64  String
─────┼───────────────
   1 │     1  Lawyer
   2 │     2  Doctor
   3 │     4  Farmer

julia> semijoin(name, job, on = :ID)
2×2 DataFrame
 Row │ ID     Name
     │ Int64  String
─────┼─────────────────
   1 │     1  John Doe
   2 │     2  Jane Doe

julia> job2 = DataFrame(identifier = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])
3×2 DataFrame
 Row │ identifier  Job
     │ Int64       String
─────┼────────────────────
   1 │          1  Lawyer
   2 │          2  Doctor
   3 │          4  Farmer

julia> semijoin(name, job2, on = :ID => :identifier)
2×2 DataFrame
 Row │ ID     Name
     │ Int64  String
─────┼─────────────────
   1 │     1  John Doe
   2 │     2  Jane Doe

julia> semijoin(name, job2, on = [:ID => :identifier])
2×2 DataFrame
 Row │ ID     Name
     │ Int64  String
─────┼─────────────────
   1 │     1  John Doe
   2 │     2  Jane Doe
```
"""
semijoin(df1::AbstractDataFrame, df2::AbstractDataFrame;
         on::Union{<:OnType, AbstractVector} = Symbol[], makeunique::Bool=false,
         validate::Union{Pair{Bool, Bool}, Tuple{Bool, Bool}}=(false, false),
         matchmissing::Symbol=:error) =
    _join(df1, df2, on=on, kind=:semi, makeunique=makeunique,
          indicator=nothing, validate=validate,
          left_rename=identity, right_rename=identity, matchmissing=matchmissing)

"""
    antijoin(df1, df2; on, makeunique=false, validate=(false, false), matchmissing=:error)

Perform an anti join of two data frame objects and return a `DataFrame`
containing the result. An anti join returns the subset of rows of `df1` that do
not match with the keys in `df2`.

The order of rows in the result is undefined and may change in the future releases.

# Arguments
- `df1`, `df2`: the `AbstractDataFrames` to be joined

# Keyword Arguments
- `on` : A column name to join `df1` and `df2` on. If the columns on which
  `df1` and `df2` will be joined have different names, then a `left=>right`
  pair can be passed. It is also allowed to perform a join on multiple columns,
  in which case a vector of column names or column name pairs can be passed
  (mixing names and pairs is allowed).
- `makeunique` : if `false` (the default), an error will be raised
  if duplicate names are found in columns not joined on;
  if `true`, duplicate names will be suffixed with `_i`
  (`i` starting at 1 for the first duplicate).
- `validate` : whether to check that columns passed as the `on` argument
   define unique keys in each input data frame (according to `isequal`).
   Can be a tuple or a pair, with the first element indicating whether to
   run check for `df1` and the second element for `df2`.
   By default no check is performed.
- `matchmissing` : if equal to `:error` throw an error if `missing` is present
  in `on` columns; if equal to `:equal` then `missing` is allowed and missings are
  matched (`isequal` is used for comparisons of rows for equality)

It is not allowed to join on columns that contain `NaN` or `-0.0` in real or
imaginary part of the number. If you need to perform a join on such values use
CategoricalArrays.jl and transform a column containing such values into a
`CategoricalVector`.

When merging `on` categorical columns that differ in the ordering of their
levels, the ordering of the left data frame takes precedence over the ordering
of the right data frame.

See also: [`innerjoin`](@ref), [`leftjoin`](@ref), [`rightjoin`](@ref),
          [`outerjoin`](@ref), [`semijoin`](@ref), [`crossjoin`](@ref).

# Examples
```julia
julia> name = DataFrame(ID = [1, 2, 3], Name = ["John Doe", "Jane Doe", "Joe Blogs"])
3×2 DataFrame
 Row │ ID     Name
     │ Int64  String
─────┼──────────────────
   1 │     1  John Doe
   2 │     2  Jane Doe
   3 │     3  Joe Blogs

julia> job = DataFrame(ID = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])
3×2 DataFrame
 Row │ ID     Job
     │ Int64  String
─────┼───────────────
   1 │     1  Lawyer
   2 │     2  Doctor
   3 │     4  Farmer

julia> antijoin(name, job, on = :ID)
1×2 DataFrame
 Row │ ID     Name
     │ Int64  String
─────┼──────────────────
   1 │     3  Joe Blogs

julia> job2 = DataFrame(identifier = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])
3×2 DataFrame
 Row │ identifier  Job
     │ Int64       String
─────┼────────────────────
   1 │          1  Lawyer
   2 │          2  Doctor
   3 │          4  Farmer

julia> antijoin(name, job2, on = :ID => :identifier)
1×2 DataFrame
 Row │ ID     Name
     │ Int64  String
─────┼──────────────────
   1 │     3  Joe Blogs

julia> antijoin(name, job2, on = [:ID => :identifier])
1×2 DataFrame
 Row │ ID     Name
     │ Int64  String
─────┼──────────────────
   1 │     3  Joe Blogs
```
"""
antijoin(df1::AbstractDataFrame, df2::AbstractDataFrame;
         on::Union{<:OnType, AbstractVector} = Symbol[], makeunique::Bool=false,
         validate::Union{Pair{Bool, Bool}, Tuple{Bool, Bool}}=(false, false),
         matchmissing::Symbol=:error) =
    _join(df1, df2, on=on, kind=:anti, makeunique=makeunique,
          indicator=nothing, validate=validate,
          left_rename=identity, right_rename=identity,
          matchmissing=matchmissing)

"""
    crossjoin(df1, df2, dfs...; makeunique = false)

Perform a cross join of two or more data frame objects and return a `DataFrame`
containing the result. A cross join returns the cartesian product of rows from
all passed data frames, where the first passed data frame is assigned to the
dimension that changes the slowest and the last data frame is assigned to the
dimension that changes the fastest.

# Arguments
- `df1`, `df2`, `dfs...` : the `AbstractDataFrames` to be joined

# Keyword Arguments
- `makeunique` : if `false` (the default), an error will be raised
  if duplicate names are found in columns not joined on;
  if `true`, duplicate names will be suffixed with `_i`
  (`i` starting at 1 for the first duplicate).

If more than two data frames are passed, the join is performed
recursively with left associativity.

See also: [`innerjoin`](@ref), [`leftjoin`](@ref), [`rightjoin`](@ref),
          [`outerjoin`](@ref), [`semijoin`](@ref), [`antijoin`](@ref).

# Examples
```julia
julia> df1 = DataFrame(X=1:3)
3×1 DataFrame
 Row │ X
     │ Int64
─────┼───────
   1 │     1
   2 │     2
   3 │     3

julia> df2 = DataFrame(Y=["a", "b"])
2×1 DataFrame
 Row │ Y
     │ String
─────┼────────
   1 │ a
   2 │ b

julia> crossjoin(df1, df2)
6×2 DataFrame
 Row │ X      Y
     │ Int64  String
─────┼───────────────
   1 │     1  a
   2 │     1  b
   3 │     2  a
   4 │     2  b
   5 │     3  a
   6 │     3  b
```
"""
function crossjoin(df1::AbstractDataFrame, df2::AbstractDataFrame; makeunique::Bool=false)
    _check_consistency(df1)
    _check_consistency(df2)
    r1, r2 = size(df1, 1), size(df2, 1)
    colindex = merge(index(df1), index(df2), makeunique=makeunique)
    cols = Any[[repeat(c, inner=r2) for c in eachcol(df1)];
               [repeat(c, outer=r1) for c in eachcol(df2)]]
    return DataFrame(cols, colindex, copycols=false)
end

crossjoin(df1::AbstractDataFrame, df2::AbstractDataFrame, dfs::AbstractDataFrame...;
          makeunique::Bool=false) =
    crossjoin(crossjoin(df1, df2, makeunique=makeunique), dfs..., makeunique=makeunique)

# an explicit error is thrown as join was supported in the past
Base.join(df1::AbstractDataFrame, df2::AbstractDataFrame, dfs::AbstractDataFrame...;
          on::Union{<:OnType, AbstractVector} = Symbol[],
          kind::Symbol = :inner, makeunique::Bool=false,
          indicator::Union{Nothing, Symbol} = nothing,
          validate::Union{Pair{Bool, Bool}, Tuple{Bool, Bool}}=(false, false)) =
    throw(ArgumentError("join function for data frames is not supported. Use innerjoin, " *
                        "leftjoin, rightjoin, outerjoin, semijoin, antijoin, or crossjoin"))


text:
"""
    hierarchical_fit(eventsTimesList, maxT, its)

# Arguments

* `eventTimesList` An array of arrays of events.
* `maxT` The boundary time over which the events were observed.
* `its` Number for of iterations to sample for.

"""
function hierarchical_fit(eventTimesList, maxT::Number, its::Int64)
    nSamples = length(eventTimesList)
    nEvents = map(length, eventTimesList)
    eventDifferences = map(event_difference_list, eventTimesList)
    maxTDifferences = map(x-> maxT .- x, eventTimesList)

    bgSamples = zeros(Float64, its)
    kappaSamples = zeros(Float64, its)
    kernSamples = zeros(Float64, its)

    bgSamples[1] = rand()
    kappaSamples[1] = rand()
    kernSamples[1] = rand()

    kernSample = kernSamples[1]
    kernDistribution::Exponential{Float64} = Distributions.Exponential(1/kernSample)
    kernFunc(x::Number) = Distributions.pdf(kernDistribution, x)

    sampleData = map(j-> sample_parents(eventTimesList[j], bgSamples[1], kappaSamples[1], kernFunc, eventDifferences[j]), 1:nSamples)

    for i in 2:its
        totalBG = mapreduce(x->x[2], +, sampleData)
        bgSamples[i] = Distributions.rand(Distributions.Gamma(0.01 + totalBG, 1/(0.01+nSamples))) / maxT

        H_tildes = map(x->sum(cdf.(kernDistribution, x)), maxTDifferences)
        totalChildEvents = mapreduce(x -> sum(x[3]), +, sampleData)
        kappaSamples[i] = Distributions.rand(Distributions.Gamma(0.01 + totalChildEvents, 1/(0.01 + sum(H_tildes))))

        totalShiftTimes = mapreduce(x -> sum(x[4]), +, sampleData)
        nShiftTimes = mapreduce(x -> length(x[4]), +, sampleData)
        kernSample = Distributions.rand(Distributions.Gamma(0.01 + nShiftTimes, 1/(0.01 + totalShiftTimes)))
        kernSamples[i] = kernSample
        kernDistribution = Distributions.Exponential(1/kernSamples[i])

        sampleData = map(j-> sample_parents(eventTimesList[j], bgSamples[i], kappaSamples[i], kernFunc, eventDifferences[j]), 1:nSamples)

    end

    (bgSamples, kappaSamples, kernSamples)
end

ext: jl


text:
function count_ops(x::SymbolicUtils.Mul)
    ops = x.coeff == 1 ? 0 : 1
    for (term, pow) in pairs(x.dict)
        ops += 1 + count_ops(term) + (pow == 1 ? 0 : 1)
    end
    return ops - 1 # The for loop overcounts operations by one
end
function count_ops(x::SymbolicUtils.Add)
    ops = x.coeff == 0 ? 0 : 1
    for (term, mult) in pairs(x.dict)
        ops += 1 + count_ops(term) + (mult == 1 ? 0 : 1)
    end
    return ops - 1 # The for loop overcounts operations by one
end
function count_ops(x::SymbolicUtils.Term)
    ops = 1
    for a in x.arguments
        ops += count_ops(a)
    end
    return ops
end
count_ops(x::SymbolicUtils.Pow) = 1 + count_ops(x.base) + count_ops(x.exp)
count_ops(::SymbolicUtils.Sym) = 0
count_ops(::Int) = 0

struct TrialRewriter
    trial_rewrites
    simplify_rewrite
    complexity_measure

    function TrialRewriter(trial_rewrites, simplify_rewrite,
            complexity_measure=count_ops)

        return new(trial_rewrites, simplify_rewrite, complexity_measure)
    end
end

function (tr::TrialRewriter)(x)
    complexity = tr.complexity_measure(x)
    for trial_rewrite in tr.trial_rewrites
        trial_walk = SymbolicUtils.Postwalk(SymbolicUtils.PassThrough(trial_rewrite))
        y = tr.simplify_rewrite(trial_walk(x))

        # y shouldn't ever be able to be nothing, but might as well check.
        if y !== nothing
            new_complexity = tr.complexity_measure(y)
            if new_complexity < complexity
                complexity = new_complexity
                x = y
            end
        end
    end
    return x
end

full_simplify_trial_rewrites = [
    # Double angle formulas
    @rule(sin(2 * ~x) => 2 * sin(~x) * cos(~x))
    @rule(cos(2 * ~x) => cos(~x)^2 + -1*sin(~x)^2)
    @rule(tan(2 * ~x) => 2*tan(~x) / (1 + -1*tan(x)^2))

    # TODO: Half angle formulas. Somewhat tricky because there are branches depending
    # on the quadrant that the angle occupies.

    # Trig sum formulas
    # Difference formulas should be derivable using these rules and other basic
    # trig rules.
    @rule(sin(~x + ~y) => sin(~x)*cos(~y) + cos(~x)*sin(~y))
    @rule(cos(~x + ~y) => cos(~x)*cos(~y) + -1*sin(~x)*sin(~y))
    @rule(tan(~x + ~y) => (tan(~x) + tan(~y)) / (1 - tan(~x)*tan(~y)))

    # Hyperbolic trig functions
    @acrule(exp(~x) +  1 => 2 * exp(1//2 * ~x) * cosh(1//2 * ~x))
    @acrule(exp(~x) + -1 => 2 * exp(1//2 * ~x) * sinh(1//2 * ~x))
]

extra_simplify_rules = [
    # Exponential rules
    @rule(exp(~x)^(~y) => exp(~x * ~y))
    @acrule(exp(~x)*exp(~y) => exp(~x + ~y))

    # Trig rules
    @rule(cos(-1*~x) => cos(~x))
    @rule(sin(-1*~x) => -1*sin(~x))
    @rule(tan(-1*~x) => -1*tan(~x))

]
extra_rewritter =
    SymbolicUtils.Postwalk(
    SymbolicUtils.PassThrough(
    SymbolicUtils.Chain(extra_simplify_rules)))

# This construction with mutiple Fixpoint rewritters seems terrible for
# performance. It should be rethought eventually.
simplify_rewritter =
    SymbolicUtils.Fixpoint(
    SymbolicUtils.Chain((extra_rewritter, SymbolicUtils.serial_simplifier)))

full_simplify_rewriter =
    TrialRewriter(
        full_simplify_trial_rewrites,
        simplify_rewritter)

fullsimplify(ex) = simplify(ex, rewriter=full_simplify_rewriter)


text:
module ExpressionIterator
"""Iterator of Substrings on which each contains a line"""
each_line_iterator(s::AbstractString) = Channel(c->begin
            for substring in split(s, "\n", limit=0, keepempty=false)
                put!(c, substring)
            end
            end, ctype=SubString{String})

struct EachExpressionAndOffset
    original::String
    next_line_offset_it # `lines_offsets(…)`
end

""" Get [begin, end] offsets for each line (end is inclusive)"""
lines_offsets(s::AbstractString) = map(line->begin
                                       real_offset=line.offset +1 # SubString.offset + 1 ↔ String.index
                                       (real_offset, real_offset + length(line) -1) # TODO(felipe) use `prevind` ?
                                       end, each_line_iterator(s))

"""points to info about "where should we get our *next* expression"""
struct _EachExpressionAndOffsetState
    maybe_next_line_offset # as in ` = iterate(…)`
    offset_where_last_expression_ends::Union{Nothing,Int}
end

_EachExpressionAndOffsetState(iter::EachExpressionAndOffset) = _EachExpressionAndOffsetState(
    Base.iterate(iter.next_line_offset_it),
    nothing)

function Base.iterate(iter::EachExpressionAndOffset, state=_EachExpressionAndOffsetState(iter)) #::Union{Nothing, Tuple{EachExpressionAndOffset, _EachExpressionAndOffsetState}}
    if state.maybe_next_line_offset == nothing
        # no more lines → we're done
        return nothing
    end

    # line ↔ [begin, end] offsets that describe a line

    # move next-line-wards …
    (current_line, current_line_state) = state.maybe_next_line_offset
    (line_begin, line_end) = current_line
    maybe_next_line = Base.iterate(iter.next_line_offset_it, current_line_state)

    # … at least we catch up with last parsed expression
    if state.offset_where_last_expression_ends !== nothing
        while line_begin < state.offset_where_last_expression_ends
            if maybe_next_line isa Nothing
                # ran out of lines with the last expression → nothing left to parse
                return nothing
            end
            # update current line info
            (current_line, current_line_state) = maybe_next_line
            (line_begin, line_end) = current_line

            # update next line info
            maybe_next_line = Base.iterate(iter.next_line_offset_it, current_line_state)
        end
    end

    (ex, i_for_end_of_expression) = Meta.parse(iter.original, line_begin)
    iter_value=(ex, line_begin, line_end)
    iter_state=_EachExpressionAndOffsetState(maybe_next_line,
                                     i_for_end_of_expression)
    return (iter_value, iter_state)
end

each_expression_and_offset(original::AbstractString) = EachExpressionAndOffset(original, lines_offsets(original))

function Base.length(iter::EachExpressionAndOffset)
    count=0
    for ex in iter
        count += 1
    end
    count
end

each_expression(original::AbstractString) = map(ex_off->begin
                                                (ex, line_begin, line_end) = ex_off
                                                ex
                                                end, each_expression_and_offset(original))

end # module ExpressionIterator


text:
#######################################################################
#
# An example of creating Excel Scatter charts with  XlsxWriter.
#
# Original Python Copyright 2013-2016, John McNamara, jmcnamara@cpan.org
# https://github.com/jmcnamara/XlsxWriter

using XlsxWriter

function test()

	wb = Workbook("chart_scatter.xlsx")
	ws = add_worksheet!(wb)
	bold = add_format!(wb, Dict("bold"=> 1))

	# Add the worksheet data that the charts will refer to.
	headings = ["Number", "Batch 1", "Batch 2"]
	data = [
		[2 3 4 5 6 7];
		[10 40 50 20 10 50];
		[30 60 70 50 40 30];
	]

	write_row!(ws, "A1", headings, fmt=bold)
	write_column!(ws, "A2", data[1,:])
	write_column!(ws, "B2", data[2,:])
	write_column!(ws, "C2", data[3,:])

	#######################################################################
	#
	# Create a new scatter chart.
	#
	chart1 = add_chart!(wb, Dict("type"=> "scatter"))

	# Configure the first series.
	add_series!(chart1, Dict(
		"name"=> "=Sheet1!\$B\$1",
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=> "=Sheet1!\$B\$2:\$B\$7",
	))

	# Configure second series. Note use of alternative syntax to define ranges.
	add_series!(chart1, Dict(
		"name"=>       ["Sheet1", 0, 2],
		"categories"=> ["Sheet1", 1, 0, 6, 0],
		"values"=>     ["Sheet1", 1, 2, 6, 2],
	))

	# Add a chart title and some axis labels.
	set_title!(chart1, Dict("name"=> "Results of sample analysis"))
	set_x_axis!(chart1, Dict("name"=> "Test number"))
	set_y_axis!(chart1, Dict("name"=> "Sample length (mm)"))

	# Set an Excel chart style.
	set_style!(chart1, 11)

	# Insert the chart into the worksheet (with an offset).
	insert_chart!(ws, "D2", chart1, Dict("x_offset"=> 25, "y_offset"=> 10))

	#######################################################################
	#
	# Create a scatter chart sub-type with straight lines and markers.
	#
	chart2 = add_chart!(wb, Dict("type"=> "scatter",
								 "subtype"=> "straight_with_markers"))

	# Configure the first series.
	add_series!(chart2, Dict(
		"name"=>       "=Sheet1!\$B\$1",
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$B\$2:\$B\$7",
	))

	# Configure second series.
	add_series!(chart2, Dict(
		"name"=>       "=Sheet1!\$C\$1",
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$C\$2:\$C\$7",
	))

	# Add a chart title and some axis labels.
	set_title!(chart2, Dict("name"=> "Straight line with markers"))
	set_x_axis!(chart2, Dict("name"=> "Test number"))
	set_y_axis!(chart2, Dict("name"=> "Sample length (mm)"))

	# Set an Excel chart style.
	set_style!(chart2, 12)

	# Insert the chart into the worksheet (with an offset).
	insert_chart!(ws, "D18", chart2, Dict("x_offset"=> 25, "y_offset"=> 10))

	#######################################################################
	#
	# Create a scatter chart sub-type with straight lines and no markers.
	#
	chart3 = add_chart!(wb, Dict("type"=> "scatter",
								 "subtype"=> "straight"))

	# Configure the first series.
	add_series!(chart3, Dict(
		"name"=>       "=Sheet1!\$B\$1",
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$B\$2:\$B\$7",
	))

	# Configure second series.
	add_series!(chart3, Dict(
		"name"=>       "=Sheet1!\$C\$1",
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$C\$2:\$C\$7",
	))

	# Add a chart title and some axis labels.
	set_title!(chart3, Dict("name"=> "Straight line"))
	set_x_axis!(chart3, Dict("name"=> "Test number"))
	set_y_axis!(chart3, Dict("name"=> "Sample length (mm)"))

	# Set an Excel chart style.
	set_style!(chart3, 13)

	# Insert the chart into the worksheet (with an offset).
	insert_chart!(ws, "D34", chart3, Dict("x_offset"=> 25, "y_offset"=> 10))

	#######################################################################
	#
	# Create a scatter chart sub-type with smooth lines and markers.
	#
	chart4 = add_chart!(wb, Dict("type"=> "scatter",
								 "subtype"=> "smooth_with_markers"))

	# Configure the first series.
	add_series!(chart4, Dict(
		"name"=>       "=Sheet1!\$B\$1",
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$B\$2:\$B\$7",
	))

	# Configure second series.
	add_series!(chart4, Dict(
		"name"=>       "=Sheet1!\$C\$1",
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$C\$2:\$C\$7",
	))

	# Add a chart title and some axis labels.
	set_title!(chart4, Dict("name"=> "Smooth line with markers"))
	set_x_axis!(chart4, Dict("name"=> "Test number"))
	set_y_axis!(chart4, Dict("name"=> "Sample length (mm)"))

	# Set an Excel chart style.
	set_style!(chart4, 14)

	# Insert the chart into the worksheet (with an offset).
	insert_chart!(ws, "D50", chart4, Dict("x_offset"=> 25, "y_offset"=> 10))

	#######################################################################
	#
	# Create a scatter chart sub-type with smooth lines and no markers.
	#
	chart5 = add_chart!(wb, Dict("type"=> "scatter",
								 "subtype"=> "smooth"))

	# Configure the first series.
	add_series!(chart5, Dict(
		"name"=>       "=Sheet1!\$B\$1",
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$B\$2:\$B\$7",
	))

	# Configure second series.
	add_series!(chart5, Dict(
		"name"=>       "=Sheet1!\$C\$1",
		"categories"=> "=Sheet1!\$A\$2:\$A\$7",
		"values"=>     "=Sheet1!\$C\$2:\$C\$7",
	))

	# Add a chart title and some axis labels.
	set_title!(chart5, Dict("name"=> "Smooth line"))
	set_x_axis!(chart5, Dict("name"=> "Test number"))
	set_y_axis!(chart5, Dict("name"=> "Sample length (mm)"))

	# Set an Excel chart style.
	set_style!(chart5, 15)

	# Insert the chart into the worksheet (with an offset).
	insert_chart!(ws, "D66", chart5, Dict("x_offset"=> 25, "y_offset"=> 10))

	close(wb)
	isfile("chart_scatter.xlsx")
end

test()


text:
module Pressure

using JuLIP
using JuLIP.Potentials: pairs
using MultiHZ
import JuLIP: energy, forces, hessian_pos, cutoff

export cal_Pressure

function cal_Pressure(at::AbstractAtoms)
    dim = get_data(at, "Dimension")
    Lx = get_data(at, "Length")
    Ly = get_data(at, "Height")

    α = get_data(at, "Alpha")
    α = [α α; α α]
    k = 1.0
    k = [k k; k k]

    r = at.M
    r1 = minimum(r)
    r2 = maximum(r)
    r0 = [2*r1 r1+r2; r1+r2 2*r2]

    z = [1, 2]
    V = MultiHZ.Multi_HZ(z, k, α, r0)

    p = zeros(dim, dim)

    for (i, j, r, R) in pairs(at, cutoff(V))
        a = V.Z2idx[at.Z[i]]
        b = V.Z2idx[at.Z[j]]
        f = 0.5 * grad(V.V[a,b], r, R)

        if dim == 2
            p[1,1] += f[1]*R[1]
            p[1,2] += f[1]*R[2]
            p[2,1] += f[2]*R[1]
            p[2,2] += f[2]*R[2]
        end

    end

    if dim == 2
        p = -p/Lx/Ly
    end

    set_data!(at, "Pressure", p)
    return p

end

end


text:
using StanSample, MCMCChains, CSV

df = CSV.read(joinpath(@__DIR__, "..",  "..", "data",  "Kline.csv"), DataFrame);

# New col log_pop, set log() for population data
df[!, :log_pop] = map((x) -> log(x), df[!, :population]);
df[!, :society] = 1:10;

m12_6s = "
data {
    int N;
    int N_societies;
    int total_tools[N];
    real logpop[N];
    int society[N];
}
parameters{
    real a;
    real bp;
    vector[N_societies] a_society;
    real<lower=0> sigma_society;
}
model{
    vector[N_societies] mu;
    sigma_society ~ cauchy( 0 , 1 );
    a_society ~ normal( 0 , sigma_society );
    bp ~ normal( 0 , 1 );
    a ~ normal( 0 , 10 );
    for ( i in 1:N ) {
        mu[i] = a + a_society[society[i]] + bp * logpop[i];
        mu[i] = exp(mu[i]);
    }
    total_tools ~ poisson( mu );
}
";

# Define the Stanmodel and set the output format to :mcmcchains.

m_12_6s = SampleModel("m12.6s",  m12_6s);

# Input data for cmdstan

m12_6_data = Dict("N" => size(df, 1), "N_societies" => 10,
"total_tools" => df[!, :total_tools], "logpop" => df[!, :log_pop],
"society" => df[!, :society]);

# Sample using cmdstan

rc = stan_sample(m_12_6s, data=m12_6_data);

# Describe the draws

if success(rc)
  chn = read_samples(m_12_6s; output_format=:mcmcchains)
  describe(chn)
end


text:
struct Ship
    direction::Char
    x::Int64
    y::Int64
end

function parse_instructions(lines)
    instructions = NamedTuple{(:action, :value), Tuple{Char, Int64}}[]
    for line in lines
        a = line[1]
        v = parse(Int64, line[2:end])
        push!(instructions, (action=a, value=v))
    end
    instructions
end

function direction_to_deg(direction)
    if direction == 'N'
        90
    elseif direction == 'W'
        180
    elseif direction == 'S'
        270
    elseif direction == 'E'
        0
    end
end

function deg_to_direction(degree)
    if degree == 0
        'E'
    elseif degree == 90
        'N'
    elseif degree == 180
        'W'
    elseif degree == 270
        'S'
    end
end

function turn(ship, value)
    deg = direction_to_deg(ship.direction)
    deg = mod(deg + value, 360)
    Ship(deg_to_direction(deg), ship.x, ship.y)
end

function move(ship, instruction)
    if instruction.action == 'N'
        Ship(ship.direction, ship.x, ship.y + instruction.value)
    elseif instruction.action == 'S'
        Ship(ship.direction, ship.x, ship.y - instruction.value)
    elseif instruction.action == 'E'
        Ship(ship.direction, ship.x + instruction.value, ship.y)
    elseif instruction.action == 'W'
        Ship(ship.direction, ship.x - instruction.value, ship.y)
    elseif instruction.action == 'L'
        turn(ship, instruction.value)
    elseif instruction.action == 'R'
        turn(ship, -instruction.value)
    elseif instruction.action == 'F'
        move(ship, (action=ship.direction, value=instruction.value))
    end
end

function navigate(instructions)
    initial = Ship('E', 0, 0)
    ship = foldl(move, instructions; init=initial)
    abs(ship.x - initial.x) + abs(ship.y - initial.y)
end

function main()
    lines = readlines("../resources/input_12.txt")
    instructions = parse_instructions(lines)
    result = navigate(instructions)
    println(result)
end

main()


text:
module testUtils_states
using Test
using LinearAlgebra

include(dirname(dirname(@__FILE__))*"\\src\\Examples\\Utils_states.jl")
using .Utils_states

@testset "take_Pᵀ" begin
    d = Dict()
    d[1] = (2,2)
    d[2] = (2,3)
    d[3] = (3,3)
    d[4] = (3,4)

    ρ_test = Dict()
    ρ_test[1] =  [1  2   9  10
                  5  6  13  14
                  3  4  11  12
                  7  8  15  16]

    ρ_test[2] =     [   1   2   3  19  20  21
                        7   8   9  25  26  27
                        13  14  15  31  32  33
                        4   5   6  22  23  24
                        10  11  12  28  29  30
                        16  17  18  34  35  36]
    ρ_test[3] =    [ 1   2   3  28  29  30  55  56  57
                    10  11  12  37  38  39  64  65  66
                    19  20  21  46  47  48  73  74  75
                    4   5   6  31  32  33  58  59  60
                    13  14  15  40  41  42  67  68  69
                    22  23  24  49  50  51  76  77  78
                    7   8   9  34  35  36  61  62  63
                    16  17  18  43  44  45  70  71  72
                    25  26  27  52  53  54  79  80  81]

    ρ_test[4] =      [  1   2   3   4  49  50  51  52   97   98   99  100
                        13  14  15  16  61  62  63  64  109  110  111  112
                        25  26  27  28  73  74  75  76  121  122  123  124
                        37  38  39  40  85  86  87  88  133  134  135  136
                        5   6   7   8  53  54  55  56  101  102  103  104
                        17  18  19  20  65  66  67  68  113  114  115  116
                        29  30  31  32  77  78  79  80  125  126  127  128
                        41  42  43  44  89  90  91  92  137  138  139  140
                        9  10  11  12  57  58  59  60  105  106  107  108
                        21  22  23  24  69  70  71  72  117  118  119  120
                        33  34  35  36  81  82  83  84  129  130  131  132
                        45  46  47  48  93  94  95  96  141  142  143  144]

    gen_ρ(d) =  reshape(1:(d[1]*d[2])^2,d[1]*d[2],d[1]*d[2])

    ρ = gen_ρ(d[1])
    Pᵀρ =  Utils_states.take_Pᵀ(ρ, 2, d[1])

    for i in 1:4
        @test Utils_states.take_Pᵀ(gen_ρ(d[i]), 1, d[i]) == ρ_test[i]
    end

    d = (3,2)
    ρ = reshape(1:36,6,6)
    @test ρ == Utils_states.take_Pᵀ(ρ,0,d)
    ρᵀ¹ = Utils_states.take_Pᵀ(ρ,1,d)
    ρᵀ² = Utils_states.take_Pᵀ(ρ,2,d)
    @test ρᵀ¹ ==   [1   2  13  14  25  26
                    7   8  19  20  31  32
                    3   4  15  16  27  28
                    9  10  21  22  33  34
                    5   6  17  18  29  30
                    11  12  23  24  35  36]
    @test ρᵀ² == [  1   7   3   9   5  11
                    2   8   4  10   6  12
                    13  19  15  21  17  23
                    14  20  16  22  18  24
                    25  31  27  33  29  35
                    26  32  28  34  30  36]

end

@testset "isPPᵀ" begin
    ρ_1=   [ 0.0833333  0.0        0.0        0.0        0.0833333  0.0        0.0       0.0        0.0833333
             0.0        0.0833333  0.0        0.0        0.0        0.0        0.0       0.0        0.0
             0.0        0.0        0.0833333  0.0        0.0        0.0        0.0       0.0        0.0
             0.0        0.0        0.0        0.0833333  0.0        0.0        0.0       0.0        0.0
             0.0833333  0.0        0.0        0.0        0.0833333  0.0        0.0       0.0        0.0833333
             0.0        0.0        0.0        0.0        0.0        0.0833333  0.0       0.0        0.0
             0.0        0.0        0.0        0.0        0.0        0.0        0.208333  0.0        0.161374
             0.0        0.0        0.0        0.0        0.0        0.0        0.0       0.0833333  0.0
             0.0833333  0.0        0.0        0.0        0.0833333  0.0        0.161374  0.0        0.208333]

    ρ_2=   [ 0.5  0.0  0.0  0.0
             0.0  0.0  0.0  0.0
             0.0  0.0  0.0  0.0
             0.0  0.0  0.0  0.5]

    @test !Utils_states.isPPᵀ(ρ_1,(3,3))
    @test Utils_states.isPPᵀ(ρ_2,(2,2))

    ρ_true = ones(4,4)
    ρ_false = reshape(1:16,4,4)
    d = (2,2)
    @test Utils_states.isPPᵀ(ρ_true,d) == true
    @test Utils_states.isPPᵀ(ρ_false,d) == false

end

@testset "maketraceone" begin
    @test tr(Utils_states.maketraceone(reshape(1:81,9,9))) == 1
    @test tr(Utils_states.maketraceone(reshape(1:36,6,6))) == 1
    @test tr(Utils_states.maketraceone(rand(4,4))) ≈ 1
    @test tr(Utils_states.maketraceone(rand(4,4))) ≈ 1
    @test tr(Utils_states.maketraceone(rand(5,5))) ≈ 1
end

end

max_stars_repo_name: JAndriesJ/sep-rank

max_issues_repo_name: JAndriesJ/sep-rank

max_forks_repo_name: JAndriesJ/sep-rank


text:
# The types defined below are used for specifying the options of the word embedding training
abstract type Option end

abstract type InitializationMethod <: Option end
struct RandomInited <: InitializationMethod
    # Initialize the embedding randomly
end
const random_inited = RandomInited()

abstract type NetworkType <: Option end
struct NaiveSoftmax <: NetworkType
    # |V| outputs softmax
end
struct HuffmanTree <: NetworkType
    # Predicate step by step on the huffman tree
end
const naive_softmax = NaiveSoftmax()
const huffman_tree = HuffmanTree()

mutable struct WordEmbedding
    vocabulary::Array{String}
    embedding::Dict{String, Array{Float64}}
    classification_tree::TreeNode
    distribution::Dict{String, Float64}
    codebook::Dict{String, Vector{Int}}

    init_type::InitializationMethod
    network_type::NetworkType
    dimension::Int
    lsize::Int    # left window size in training
    rsize::Int    # right window size
    trained_times::Dict{String,Int}
    trained_count::Int
    corpus_size::Int
    subsampling::Float64
    init_learning_rate::Float64
    iter::Int
    min_count::Int
end

function WordEmbedding(dim::Int, init_type::InitializationMethod, network_type::NetworkType; lsize=5, rsize=5, subsampling=1e-5, init_learning_rate=0.025, iter=5, min_count=5)
    if dim <= 0 || lsize <= 0 || rsize <= 0
        throw(ArgumentError("dimension should be a positive integer"))
    end
    WordEmbedding(String[],
                    Dict{String,Array{Float64}}(),
                    nullnode,
                    Dict{String,Array{Float64}}(),
                    Dict{String,Vector{Int}}(),
                    init_type, network_type,
                    dim,
                    lsize, rsize,
                    Dict{String,Int}(),
                    0, 0,
                    subsampling, init_learning_rate, iter, min_count)
end

function Base.show(io::IO, x::WordEmbedding)
    println(io, "Word embedding(dimension = $(x.dimension)) of $(length(x.vocabulary)) words, trained on $(x.trained_count) words")
    nothing
end

# strip embedding and retain only word vectors
function _strip(embed::WordEmbedding)
    embed.vocabulary = String[]
    embed.classification_tree = nullnode
    embed.distribution = Dict{String,Array{Float64}}()
    embed.codebook = Dict{String,Vector{Int}}()
    embed
end

# print the code book for debug
function _print_codebook(embed::WordEmbedding, N=10)
    for (word,code) in embed.codebook
        println("$word => $code")
        N -= 1
        (N > 0) || break
    end
    nothing
end

#===============================================================================
Step 1: Find Word Distribution
===============================================================================#

# function merge_distributions(dst_ref::RemoteRef, src_ref::RemoteRef)
#     word_count1,distribution1 = fetch(dst_ref)
#     word_count2,distribution2 = fetch(src_ref)

#     for (n,v) in distribution2
#         if haskey(distribution1, n)
#             distribution1[n] += v
#         else
#             distribution1[n] = v
#         end
#     end

#     ((word_count1 + word_count2), distribution1)
# end

# merge_distributions(refs) = merge_distributions(convert(Array{RemoteRef,1}, refs))
# function merge_distributions(refs::Array{RemoteRef,1})
#     result = refs
#     while length(result) > 1
#         iter_result = Array(RemoteRef, 0)
#         while length(result) > 1
#             refpair = splice!(result, 1:2)
#             push!(iter_result, remotecall(refpair[1].where, merge_distributions, refpair[1], refpair[2]))
#         end
#         isempty(result) || push!(iter_result, pop!(result))
#         result = iter_result
#     end
#     fetch(result[1])
# end

function get_distribution(corpus_file::IO)
    distribution = Dict{String,Float64}()
    word_count = 0

    for i in words_of(corpus_file)
        if !haskey(distribution, i)
            distribution[i] = 1
        else
            distribution[i] += 1
        end
        word_count += 1
    end

    (word_count, distribution)
end

function get_distribution(corpus_filename::String)
    open(corpus_filename, "r") do fs
        return get_distribution(fs)
    end
end

# function get_distribution(b::Block)
#     count_refs = pmap(get_distribution, b; fetch_results=false)
#     merge_distributions(count_refs)
# end

function strip_infrequent(distribution::Dict{String,Float64}, min_count::Int)
    stripped_distr = Dict{String,Float64}()
    word_count = 0

    for (k,v) in distribution
        if v >= min_count
            word_count += round(Int, v)
            stripped_distr[k] = v
        end
    end

    (word_count, stripped_distr)
end

function compute_frequency!(distribution::Dict{String,Float64}, word_count::Int)
    for (k, v) in distribution
        distribution[k] /= word_count
    end
    nothing
end

function word_distribution(source::String, min_count::Int=5)
    t1 = time()

    println("Finding word distribution...")
    word_count, distribution = get_distribution(source)
    println("Word Count: $word_count, Vocabulary Size: $(length(keys(distribution)))")

    println("Stripping infrequent words...")
    word_count, distribution = strip_infrequent(distribution, min_count)
    println("Word Count: $word_count, Vocabulary Size: $(length(keys(distribution)))")

    compute_frequency!(distribution, word_count)
    println("Compute time: $(time()-t1)")

    distribution
end

#===============================================================================
Step 2: Calculate Word Embedding
===============================================================================#

# function merge_embeds(ref1::RemoteRef, ref2::RemoteRef)
#     e1 = fetch(ref1)
#     e2 = fetch(ref2)
#     for word in e1.vocabulary
#         e1.embedding[word] .+= e2.embedding[word]
#     end
#     e1.trained_count += e2.trained_count
#     println("merged $ref1 and $ref2")
#     e1
# end

# merge_embeds(refs) = merge_embeds(convert(Array{RemoteRef,1}, refs))
# function merge_embeds(refs::Array{RemoteRef,1})
#     n = length(refs)
#     result = refs
#     while length(result) > 1
#         iter_result = Array(RemoteRef, 0)
#         while length(result) > 1
#             refpair = splice!(result, 1:2)
#             push!(iter_result, remotecall(refpair[1].where, merge_embeds, refpair[1], refpair[2]))
#         end
#         isempty(result) || push!(iter_result, pop!(result))
#         result = iter_result
#     end

#     emb = fetch(result[1])
#     for word in emb.vocabulary
#         emb.embedding[word] /= n
#     end

#     emb
# end

function work_process(ser_embed::String, words_stream::WordStream, strip::Bool=false)
    embed = restore(ser_embed)
    work_process(embed, words_stream, strip)
end

function work_process(embed::WordEmbedding, words_stream::WordStream, strip::Bool=false)
    t1 = time()
    middle = embed.lsize + 1
    input_gradient = zeros(Float64, embed.dimension)
    α = embed.init_learning_rate
    trained_count = 0
    trained_times = Dict{String, Int}()

    for current_iter = 1:embed.iter
        for window = sliding_window(words_stream, lsize=embed.lsize, rsize=embed.rsize)
            trained_word = window[middle]
            trained_times[trained_word] = get(trained_times, trained_word, 0) + 1
            trained_count += 1

            if trained_count % 10000 == 0
                flen = words_stream.endpoint - words_stream.startpoint
                iter_progress = (position(words_stream.fp) - words_stream.startpoint) / flen
                progress = ((current_iter - 1) + iter_progress) / embed.iter
                println("trained on $trained_count words, progress $progress, α = $α")
                α = embed.init_learning_rate * (1 - progress)
                if α < embed.init_learning_rate * 0.0001
                    α = embed.init_learning_rate * 0.0001
                end
            end

            local_lsize = round(Int, (rand(UInt) % embed.lsize))
            local_rsize = round(Int, (rand(UInt) % embed.rsize))

            for ind = (middle - local_lsize) : (middle + local_rsize)
                (ind == middle) && continue

                target_word = window[ind]
                # discard words not presenting in the classification tree
                (haskey(embed.codebook, target_word) && haskey(embed.codebook, trained_word)) || continue

                node = embed.classification_tree::TreeNode

                fill!(input_gradient, 0.0)
                input = embed.embedding[trained_word]

                for code = embed.codebook[target_word]
                    train_one(node.data, input, code, input_gradient, α)
                    node = node.children[code]
                end

                for i = 1:embed.dimension
                    input[i] -= input_gradient[i]
                end
            end
        end
    end

    embed.trained_count = trained_count
    embed.trained_times = trained_times

    t2 = time()
    println("Finished training. Trained on $(embed.trained_count) words in $(t2-t1) seconds.")

    # strip to remove unnecessary members and make serialization faster
    strip && _strip(embed)

    embed
end

#===============================================================================
Train a corpus
Step 1: Find Word Distribution
Step 2: Calculate Word Embedding
===============================================================================#

function initialize_embedding(embed::WordEmbedding, randomly::RandomInited)
    for i in embed.vocabulary
        embed.embedding[i] = rand(1, embed.dimension) * 2 .- 1
    end
    embed
end

function initialize_network(embed::WordEmbedding, huffman::HuffmanTree)
    heap = PriorityQueue()
    for (word, freq) in embed.distribution
        node = BranchNode(Vector{BranchNode}(), word, nothing)    # the data field of leaf node is its corresponding word.
        enqueue!(heap, node, freq)
    end
    while length(heap) > 1
        (node1, freq1) = peek(heap)
        dequeue!(heap)
        (node2, freq2) = peek(heap)
        dequeue!(heap)
        newnode = BranchNode([node1, node2], LinearClassifier(2, embed.dimension), nothing) # the data field of internal node is the classifier
        enqueue!(heap, newnode, freq1 + freq2)
    end
    embed.classification_tree = dequeue!(heap)
    embed
end

# function train(embed::WordEmbedding, corpus::Block, prepared_corpus_location::AbstractString)
#     corpus = corpus |> as_io |> as_wordio
#     embed.distribution = word_distribution(corpus)
#     embed.vocabulary = collect(keys(embed.distribution))

#     initialize_embedding(embed, embed.init_type)        # initialize by the specified method
#     initialize_network(embed, embed.network_type)

#     # determine the position in the tree for every word
#     for (w, code) in leaves_of(embed.classification_tree)
#         embed.codebook[w] = code
#     end

#     # serialize the word distribution for workers to pick up
#     save(embed, prepared_corpus_location)

#     t1 = time()
#     println("Starting parallel training...")

#     # Note: subsampling is not honored here, probably not required also?
#     corpus = corpus |> words_of

#     embs = pmap((x)->work_process(prepared_corpus_location, x, true), corpus; fetch_results=false)
#     t2 = time()
#     println("Partial training done at $(t2-t1) time")

#     println("Merging results...")
#     training_result = merge_embeds(embs)
#     embed.embedding = training_result.embedding
#     embed.trained_count = training_result.trained_count
#     t3 = time()

#     println("Training complete at $(t3-t1) time")
#     embed
# end

function train(embed::WordEmbedding, corpus_filename::String)
    embed.distribution = word_distribution(corpus_filename)
    embed.vocabulary = collect(keys(embed.distribution))

    initialize_embedding(embed, embed.init_type)        # initialize by the specified method
    initialize_network(embed, embed.network_type)

    # determine the position in the tree for every word
    for (w, code) = leaves_of(embed.classification_tree)
        embed.codebook[w] = code
    end

    t1 = time()
    println("Starting sequential training...")
    open(corpus_filename, "r") do fs
        words_stream = words_of(fs, subsampling = (embed.subsampling, true, embed.distribution))
        work_process(embed, words_stream, false)
    end

    t2 = time()
    println("Training complete at $(t2-t1) time")
    embed
end


text:
export pwc_density, pwc_densities, pwc_densities!

"""
    pwc_density(x::AbstractVector)

Returns the piecewise constant probability density from the quantile particle positions.

Let the quantile particles be at positions `x₀, x₁, …, xₙ` (remember that in Julia
they are actually `x[1], x[2], ..., x[n], x[n+1]`).
The returned array `R` is such that `R[1] = R[n+2] = 0` and `R[i] = 1 / (N * (x[i] - x[i-1]))`
for the intermediate indices (this formula holds also for the first and last entry
if we assume `x[0] = -∞` and `x[n+2] = ∞`).

# Examples

```jldoctest
julia> pwc_density([0, 1, 3])
4-element Vector{Float64}:
 0.0
 0.5
 0.25
 0.0
```
"""
function pwc_density(x::AbstractVector{<:Real})
    len = length(x)
    R = Array{float(eltype(x))}(undef, len+1)
    R[1] = R[end] = 0
    d = 1 / (len-1)
    for i in 2:len
        R[i] = d / (x[i] - x[i-1])
    end
    R
end

"""
    pwc_density(x::AbstractVector, y::AbstractVector)

The returned `x_dens` is indexed as `x_dens[x_or_y, left_or_right, i]`.
"""
function pwc_density(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})
    x_len, y_len = length(x), length(y)
    T = promote_type(float(eltype(x)), float(eltype(y)))
    x_dens = Array{T, 3}(undef, (2, 2, x_len))
    y_dens = Array{T, 3}(undef, (2, 2, y_len))
    # indices of current particles being examined
    x_i::Int, y_i::Int = 1, 1
    # current densities
    x_d::T, y_d::T = zero(T), zero(T)
    # normalization factors
    x_f::T, y_f::T = 1 / (x_len - 1), 1 / (y_len - 1)
    # @debug "Entering loop" x_len y_len
    @inbounds while x_i <= x_len || y_i <= y_len
        if y_i > y_len || x_i <= x_len && x[x_i] < y[y_i]
            # @debug "Iteration x < y" x_i y_i
            x_dens[1, 1, x_i] = x_d
            x_dens[2, 1, x_i] = y_d
            if x_i < x_len
                x_d = x_f / (x[x_i + 1] - x[x_i])
            else
                x_d = zero(T)
            end
            x_dens[1, 2, x_i] = x_d
            x_dens[2, 2, x_i] = y_d
            x_i += 1
        elseif x_i > x_len || y_i <= y_len && y[y_i] < x[x_i]
            # @debug "Iteration y < x" x_i y_i
            y_dens[1, 1, y_i] = x_d
            y_dens[2, 1, y_i] = y_d
            if y_i < y_len
                y_d = y_f / (y[y_i + 1] - y[y_i])
            else
                y_d = zero(T)
            end
            y_dens[1, 2, y_i] = x_d
            y_dens[2, 2, y_i] = y_d
            y_i += 1
        else
            # @debug "Iteration x == y" x_i y_i
            x_dens[1, 1, x_i] = x_d
            x_dens[2, 1, x_i] = y_d
            y_dens[1, 1, y_i] = x_d
            y_dens[2, 1, y_i] = y_d
            if x_i < x_len
                x_d = x_f / (x[x_i + 1] - x[x_i])
            else
                x_d = zero(T)
            end
            if y_i < y_len
                y_d = y_f / (y[y_i + 1] - y[y_i])
            else
                y_d = zero(T)
            end
            x_dens[1, 2, x_i] = x_d
            x_dens[2, 2, x_i] = y_d
            y_dens[1, 2, y_i] = x_d
            y_dens[2, 2, y_i] = y_d
            x_i += 1
            y_i += 1
        end
    end
    (x_dens, y_dens)
end

"""
    pwc_densities(xs::AbstractVector{<:Real}...)
    pwc_densities(xs::Tuple{Vararg{AbstractVector{<:Real}}})

Given a list of `n` families of particles of lengths `l₁, …, lₙ`,
returns an `n`-tuple whose `s`-th entry is an `n×2×lₛ` array `densₛ`.

This array contains the information about the densities of all the species
around the particles of the `s`-th species.
This array is indexed as `[o,side,i]`, where `o∈{1, …, n}` is the index of the other
species whose density we want to know, `side` is either `1` if we want to know
the left density or `2` for the right density, and `i∈{1, …, lₛ}` is the index
of the particle.

In other words, `densₛ[o,side,i]` is the density of the `o`-th species on the left (`side=1`)
or right (`side=2`) of the `i`-th particle of the `s`-th species.

See also [`pwc_densities!`](@ref) for an inplace version.

# Examples

```jldoctest
julia> pwc_densities([0,1,2], [0,2,3])
([0.0 0.5; 0.0 0.25]

[0.5 0.5; 0.25 0.25]

[0.5 0.0; 0.25 0.5], [0.0 0.5; 0.0 0.25]

[0.5 0.0; 0.25 0.5]

[0.0 0.0; 0.5 0.0])
```
"""
function pwc_densities end

pwc_densities(xs::AbstractVector{<:Real}...) = pwc_densities(xs)

function pwc_densities(xs::Tuple{Vararg{AbstractVector{<:Real}}})
    N = nfields(xs)
    T = promote_type(float.(eltype.(xs))...)
    len = length.(xs)
    dens = map(x -> new_undef_densities(T, N, length(x)), xs)::NTuple{N, Array{T, 3}}
    pwc_densities!(dens, xs)
end

"""
    pwc_densities!(dens, xs::AbstractVector{<:Real}...)
    pwc_densities!(dens, xs::Tuple{Vararg{AbstractVector{<:Real}}})

This is an inplace version of [`pwc_densities`](@ref) that writes the result in `dens`.

See the documentation of [`pwc_densities`](@ref) for an explanation.
"""
function pwc_densities! end

pwc_densities!(dens, xs::AbstractVector{<:Real}...) = pwc_densities!(dens, xs)

function pwc_densities!(dens, xs::Tuple{Vararg{AbstractVector{<:Real}}})
    N = nfields(xs)
    T = promote_type(float.(eltype.(xs))...)
    len = length.(xs)
    # indices of current particles being examined
    ind = ones(Int, N)
    # current densities
    ds = zeros(T, N)
    # normalization factors
    fs::NTuple{N, T} = 1 ./ (len .- 1)
    # @debug "Entering loop" T len ds fs
    i_min = Int[0]
    @inbounds while true
        n_mins::Int = 0
        for i in 1:N
            if ind[i] <= len[i]
                if n_mins == 0 || xs[i][ind[i]] < xs[i_min[1]][ind[i_min[1]]]
                    i_min[1] = i
                    n_mins = 1
                elseif xs[i][ind[i]] <= xs[i_min[1]][ind[i_min[1]]]
                    n_mins += 1
                    if n_mins <= length(i_min)
                        i_min[n_mins] = i
                    else
                        push!(i_min, i)
                    end
                end
            end
        end
        if n_mins == 0
            break
        end
        i_mins = @view i_min[1:n_mins]
        # let ind = (ind...,), i_min = (i_min...,), i_mins = (i_mins...,)
        #     @debug "Iteration" ind i_min n_mins i_mins
        # end
        for i in i_mins
            # @debug "Assigning left" i ind[i] ds
            # dens[i][:, 1, ind[i]] = ds
            for s in 1:N
                dens[i][s, 1, ind[i]] = ds[s]
            end
        end
        for i in i_mins
            if ind[i] < length(xs[i])
                ds[i] = fs[i] / (xs[i][ind[i]+1] - xs[i][ind[i]])
            else
                ds[i] = zero(T)
            end
            # ds[i] = right_local_density(xs[i], ind[i], fs[i])
        end
        for i in i_mins
            # @debug "Assigning right" i ind[i] ds
            # dens[i][:, 2, ind[i]] = ds
            for s in 1:N
                dens[i][s, 2, ind[i]] = ds[s]
            end
        end
        # ind[i_mins] .+= 1
        for i in i_mins
            ind[i] += 1
        end
    end
    dens
end

new_undef_densities(T::Type, N::Int, len::Int) = Array{T, 3}(undef, (N, 2, len))