Combining Matlab's `repmat` with R's `rep` for tensors of arbitrary order

repeat.jl

function repeat{T}(v::Array{T},
                   dims::Integer...;
                   repetitions::Integer = 1)
    n = length(v)
    n_res = n * repetitions * prod(dims)
    res = Array(T, tuple(repetitions * n * dims[1], dims[2:end]...)...)
    for i in 0:(n_res - 1)
        index = mod(fld(i, repetitions), n)
        res[i + 1] = v[index + 1]
    end
    return res
end

myrepmat{T}(v::Vector{T}, a::Integer, b::Integer) = repeat(v, a, b)
myrepmat{T}(m::Matrix{T}, a::Integer, b::Integer) = repeat(m, a, b)

@assert isequal(repeat(["A", "B"], 1, repetitions = 1),
                ["A", "B"])

@assert isequal(repeat(["A", "B"], 2, repetitions = 1),
                ["A", "B", "A", "B"])

@assert isequal(repeat(["A", "B"], 1, repetitions = 2),
                ["A", "A", "B", "B"])

@assert isequal(repeat(["A", "B"], 2, repetitions = 2),
                ["A", "A", "B", "B", "A", "A", "B", "B"])

@assert isequal(repeat(["A", "B"], 1, 1, repetitions = 1),
                ["A", "B"]'')

@assert isequal(repeat(["A", "B"], 1, 2, repetitions = 1),
                ["A" "A";
                 "B" "B"])

@assert isequal(repeat(["A", "B"], 2, 1, repetitions = 1),
                ["A";
                 "B";
                 "A";
                 "B";]'')

@assert isequal(repeat(["A", "B"], 2, 2, repetitions = 1),
                ["A" "A";
                 "B" "B";
                 "A" "A";
                 "B" "B";])

@assert isequal(repeat(["A", "B"], 2, 2, repetitions = 2),
                ["A" "A";
                 "A" "A";
                 "B" "B";
                 "B" "B";
                 "A" "A";
                 "A" "A";
                 "B" "B";
                 "B" "B";])

for i in 1:3
    for j in 1:3
        @assert isequal(myrepmat(["A", "B"], i, j),
                        repmat(["A", "B"], i, j))
    end
end

for i in 1:3
    for j in 1:3
        @assert isequal(myrepmat(["A", "B"]'', i, j),
                        repmat(["A", "B"]'', i, j))
    end
end

Another approach would be repeat{T}(v::Array{T}, dims::Integer...; repetitions::Integer = 1). This would look more like repmat, but would be much more general.

As I think about it, using keyword arguments for repeat is better since it clearly distinguishes the per-item repetitions from the output size. This means that something very close to R's rep and Matlab's repmat is exactly embedded in the proposed version of repeat.

The problem with this approach is that it doesn't give the right answer for matrices. The behavior it should have is:

julia> repmat([1, 2], 2, 1)
4x1 Int64 Array:
 1
 2
 1
 2

julia> repmat([1 2], 2, 1)
2x2 Int64 Array:
 1  2
 1  2

In contrast, it has this behavior:

julia> repeat([1, 2], 2, 1)
4x1 Int64 Array:
 1
 2
 1
 2

julia> repeat([1 2], 2, 1)
4x1 Int64 Array:
 1
 2
 1
 2

While it's easy to see how to hardcode the vector and matrix cases, it's not clear to me how to generalize repmat to produce sane behavior for tensors or under repetitions.

Next set of ideas:

repeat(A, inner = [a1, a2, a3], outer = [b1, b2, b3])

# order_res = max(length(inner), length(outer))
# Retain specified singleton dimensions
# Order of Result >= Order of Input

repeat([1, 2, 3], inner = [1], outer = [1])
 => [1, 2, 3]

repeat([1, 2, 3], inner = [1], outer = [1, 1])
 <=> repeat([1, 2, 3], inner = [1, 1], outer = [1, 1])
 => [1;
     2;
     3;]

repeat([1, 2, 3], inner = [2], outer = [1])
 => [1, 1, 2, 2, 3, 3]

repeat([1, 2, 3], inner = [1], outer = [2])
 => [1, 2, 3, 1, 2, 3]

repeat([1 2;
        3 4], inner = [1, 1], outer = [1, 1])
 => [1 2;
     3 4;]

repeat([1 2;
        3 4], inner = [2, 1], outer = [1, 1])
 => [1 2;
     1 2;
     3 4;
     3 4;]

repeat([1 2;
        3 4], inner = [1, 2], outer = [1, 1])
 => [1 1 2 2;
     3 3 4 4;]

# Index in original to indices in result
#1 -> 1, 3
#2 -> 2, 4
# ...
#4 -> 6, 8

repeat([1 2;
        3 4], inner = [1, 1], outer = [2, 1])
 => [1 2;
     3 4;
     1 2;
     3 4;]

# Index in original to indices in result
#1 -> 1, 3
#2 -> 2, 4
# ...
#4 -> 6, 8

repeat([1 2;
        3 4], inner = [1, 1], outer = [1, 2])
 => [1 2 1 2;
     3 4 3 4;]

# Index in original to indices in result
#1 -> 1, 5
#2 -> 2, 6
# ...
#4 -> 4, 8

repeat([1 2;
        3 4;
        5 6;], inner = [1, 1], outer = [2, 1])
 => [1 2;
     3 4;
     5 6;
     1 2;
     3 4;
     5 6;]

# Index in original to indices in result
#1 -> 1, 4
#2 -> 2, 5
# ...
#4 -> 7, 10
# mod(7, 6) = mod(7, outer[1] * sizes[0])
# Distance between successive slots is
#  outer[i] * sizes[i - 1]
#  outer[i - 2] * sizes[i - 1]
# Offset?
# Output Sizes:
#  output_sizes[i] = sizes[i] * inner[i] * outer[i]

repeat([1 2;
        3 4;
        5 6;], inner = [1, 1], outer = [1, 2])
 => [1 2 1 2;
     3 4 3 4;
     5 6 5 6;]

# Index in original to indices in result
#1 -> 1, 7
#2 -> 2, 6
# ...
#4 -> 4, 8

# sum_over_i(outer[i] * size(A, i - 1))???


repeat([:, :, 1][1 2;
                 3 4;],
       [:, :, 2][5 6;
                 7 8;], inner = [1, 1, 2], outer = [1, 1, 1])
 => [:, :, 1][1 2;
              3 4;]
    [:, :, 2][1 2;
              3 4;]
    [:, :, 3][5 6;
              7 8;]
    [:, :, 4][5 6;
              7 8;]

#1 -> 1, 5
#2 -> 2, 6
#3 -> 3, 7
# ...
#5 -> 9, 13
# ...
#8 -> 12, 16

# i -> i + intervenors_along_previous_dimensions, ... distances

repeat([1 2], inner = [1, 1], outer = [2, 1])
 => [1 2;
     1 2;]

sizes = [1, 2]

# If you had full index set, [i, j]
# you would do:
# i, j -> ??
#1, 1 -> 1, 1 AND 2, 1
#1, 2 -> 2, 1 AND 2, 2

repeat([1 2;
        3 4;], inner = [2, 2], outer = [2, 2])
 => [1 1 2 2 1 1 2 2;
     1 1 2 2 1 1 2 2;
     3 3 4 4 3 3 4 4;
     3 3 4 4 3 3 4 4;
     1 1 2 2 1 1 2 2;
     1 1 2 2 1 1 2 2;
     3 3 4 4 3 3 4 4;
     3 3 4 4 3 3 4 4;]

# Input Index -> Output Indices
#1 -> 1, 2, 5, 6, 9, 10, 13, 14, 33, 34, 37, 38, 41, 42, 45, 46
#2 -> 3, 4, 7, 8, 11, 12, 15, 16, 35, 36, 39, 40, 43, 44, 47, 48
#3 -> 17, 18, 21, 22, 25, 26, 29, 30, 49, 50, 53, 54, 57, 58, 61, 62
#4 -> 19, 20, 23, 24, 27, 28, 31, 32, 51, 52, 55, 56, 59, 60, 63, 64
#
# Observations
# * All sequences have same gaps
#  * Estimate sequence 1, 2, 5, 6, 9, 10, 13, 14, 33, 34, 37, 38, 41, 42, 45, 46
# * Offsets are 0, 2, 16, 18
# * Need rule to produce gap sequences and offsets

# i, j -> i:SF1:max_i, j:SF2:max_j ???
# inner increases offsets, but decreases distance between successors
# outer "maintains" offsets???
# i -> inner[1] * sizes[1]
#

repeat([1, 2], inner = [1, 1], outer = [2, 1])
 => [1;
     2;
     1;
     2;]

sizes = [2, 1]