
"""

    simplex_split(;k::Int = 2, dim::Int = 3, orientations = false)

Instructions for splitting a ´dim´-dimensional simplex with a size 
reducing factor `k`.
"""
function simplex_split(;k::Int = 2, dim::Int = 3, orientations = false)

    sequences::Array{Int, 2} = tensordecomp(k, dim)
    n_seq = size(sequences, 1)

    χ1 = sequences .* (dim + 1)
    χ2 = repeat(transpose(collect(1:dim)), n_seq, 1)
    χ::Array{Int, 2} = sort(χ1 .+ χ2, dims = 2)

    matrices_simplicial_subdivision::Array{Int, 2} = zeros(Int, (dim + 1) * k^dim, k)

    for a = 1:n_seq
        tmp = ones(Int, k * (dim + 1))
        for b = 1:(dim - 1)
            tmp[(χ[a, b] + 1):(χ[a, b+1])] = (b + 1) * ones(χ[a, b + 1] - χ[a, b], 1)
        end

        tmp[(χ[a, dim] + 1):size(tmp, 1)] = (dim + 1) * ones(size(tmp, 1) - χ[a, dim], 1)

        indices = ((a - 1) * (dim + 1) + 1):(a*(dim + 1))
        matrices_simplicial_subdivision[indices, :] = reshape(tmp, dim + 1, k)
    end


    simplex_orientations = zeros(Float64, k^dim, 1)
    Χ_values = zeros(Int, dim, 1)

    for i = 1:size(simplex_orientations, 1)
        indices::Vector{Int} = ((dim + 1) * (i - 1) + 1):((dim + 1) * i)
        Χi = matrices_simplicial_subdivision[indices, :]

        ΔΧi = (Χi[2:(dim + 1), :] - Χi[1:dim, :]) .* repeat(transpose(collect(1:k)), dim, 1)

        transition_ind = round.(Int, ΔΧi * ones(k, 1))

        M = zeros(dim, dim)

        for j = 1:dim
            Χ_values[j] =  Χi[j, transition_ind[j]]
            M[j, Χ_values[j]] = 1
        end

        simplex_orientations[i] = det(M) / k^dim
    end


    if orientations # Should the orientations of the simplices also be returned?
        return matrices_simplicial_subdivision, simplex_orientations
    else
        return matrices_simplicial_subdivision
    end

end

"""
    heaviside0(x)

Returns 0 if `x < 0` and 1 if `x >= 0` . If `x` is a vector, do the same, but element-wise.
"""
function heaviside0 end

heaviside0(x::Real) = ifelse(x < 0, zero(x), ifelse(x >= 0, one(x), oftype(x, 0)))
heaviside0(x::AbstractArray{T}) where {T<:Real} = heaviside0.(x)

"""
    simplex_splitting_rules(k, dim)

Generate generic rules for generating vertices when splitting a `dim`-dimensional
simplex with a size reducing factor of k.

## Returns 
1. An array of dimension `r * k` (``R = (k + dim)!/(k! dim!) - dim - 1``), where each 
    row contains the indices of the original vertices (to be linearly combined to 
    form new vertices?); these vertices defining the new vertices of the simplicial 
    subdivision. The number of columns in this array defines how many new vertices 
    are generated each split.
2. An array of integers with dimension `(dim + 1)*k^dim x 1`. It consists of a
    stack of k^E columns of dimension E+1, each stack containing
    the indices of the vertices generating the corresponding simplex
    of the subdivision in the order given by 
    `[rows(original_vertices); rows(new_vertices)]`
"""
function simplex_splitting_rules(k::Int, dim::Int)

    # Get simplex splitting rules, ignore orientations for now.
    # What are these instructions?
    splits = simplex_split(;k = k, dim = dim, orientations = false)
    
    # Retrieve integer labels from the tensor decomposition. A way of uniquely identifying
    # the rows of splits. Each integer corresponds to a row in
    # splits (but are not necessarily in order).
    integer_labels = (splits .- 1) * ((dim + 1) .^ collect(0:(k - 1)))
    uniquelabels = unique(integer_labels)
    n_rows = size(uniquelabels, 1)

    J = indexin(integer_labels, uniquelabels)
    
    # What does I contain?
    # IndicesOfUnqiueIntegersinOriginal
    I = zeros(Int, length(uniquelabels))
    for i = 1:length(uniquelabels)
        I[i] = findfirst(integer_labels .== uniquelabels[i])
    end

    
    # What is all this aux stuff? Is there a simpler, more readable way to do this?
    Aux = splits[I, :]
    V = repeat(uniquelabels, 1, dim + 1)
    aux = repeat((((dim + 1)^k - 1) ./ dim ) * transpose(collect(0:dim)), size(uniquelabels, 1), 1)

    aux = heaviside0(-abs.(V - aux)) .* repeat(collect(1:n_rows), 1, dim + 1)
    aux = aux[findall(aux .> 0)]

    aux = round.(Int, aux)
    Caux = round.(Int, setdiff(1:n_rows, aux))
    Ipp = [aux; Caux]

    Ip = zeros(Int, size(I))

    # for any i\in Ipp, Ip(i) is the number j\in {1,...,size(Ipp,1)} such that Ipp(j)=i
    Ip[Ipp] = transpose(collect(1:size(I, 1)))
    
    # Rules for forming the strictly new vertices of the subtriangulation. 
    # The i-th row in this array contains what? 
    # Indices of the non-repeated vertices that are not part of the original simplex.
    # These are the truly new vertices generated by the splitting.
    rules_new_vertices = Aux[Caux, :]

    allverts = vcat(repeat(collect(1:(dim + 1)), 1, k), rules_new_vertices)

    # Subarrays of dimension E+1 x 1 vertically concatenated. Each subarray contains the
    # indices in the rows of NewVert corresponding to the vertices furnishing the corresponding
    # simplex (analogous to the output of SimplexSplitting, check notes).
    SubTriangulation = round.(Int, Ip[J])
    SubTriangulation = reshape(SubTriangulation, dim + 1, k^dim)

    return allverts, SubTriangulation
end

export simplex_splitting_rules