Optimal way to compute permutations in julia

Question

Consider a list [1,1,1,...,1,0,0,...,0] (an arbitrary list of zeros and ones). We want the whole possible permutations in this array, there'll be binomial(l,k) permutations (l stands for the length of the list and k for the number of ones in the list).

Right now, I have tested three different algorithms to generate the whole possible permutations, one that uses a recurrent function, one that calculates the permutations via calculating the interval number [1,...,1,0,0,...,0] to [0,0,...0,1,1,...,1] (since this can be seen as a binary number interval), and one that calculates the permutations using lexicographic order.

So far, the first two approaches fail in performance when the permutations are approx. 32. The lexicographic technique works still pretty nice (only a few miliseconds to finish).

My question is, specifically for julia, which is the best way to calculate permutations as I described earlier? I don't know too much in combinatorics, but I think a descent benchmark would be to generate all permutations from the total binomial(l,l/2)

Answer 1

As you have mentioned yourself in the comments, the case where l >> k is definitely desired. When this is the case, we can substantially improve performance by not handling vectors of length l until we really need them, and instead handle a list of indexes of the ones.

In the RAM-model, the following algorithm will let you iterate over all the combinations in space O(k^2), and time O(k^2 * binom(l,k))

Note however, that every time you generate a bit-vector from an index combination, you incur an overhead of O(l), in which you will also have the lower-bound (for all combinations) of Omega(l*binom(l,k)), and the memory usage grows to Omega(l+k^2).

The algorithm

"""
Produces all `k`-combinations of integers in `1:l` with prefix `current`, in a
lexicographical order.

# Arguments
- `current`: The current combination
- `l`: The parent set size
- `k`: The target combination size
"""
function combination_producer(l, k, current)
    if k == length(current)
        produce(current)
    else
        j = (length(current) > 0) ? (last(current)+1) : 1
        for i=j:l
            combination_producer(l, k, [current, i])
        end
    end
end

"""
Produces all combinations of size `k` from `1:l` in a lexicographical order
"""
function combination_producer(l,k)
    combination_producer(l,k, [])
end

Example

You can then iterate over all the combinations as follows:

for c in @task(combination_producer(l, k))
    # do something with c
end

Notice how this algorithm is resumable: You can stop the iteration whenever you want, and continue again:

iter = @task(combination_producer(5, 3))
for c in iter
    println(c)
    if c[1] == 2
        break
    end
end

println("took a short break")

for c in iter
    println(c)
end

This produces the following output:

[1,2,3]
[1,2,4]
[1,2,5]
[1,3,4]
[1,3,5]
[1,4,5]
[2,3,4]
took a short break
[2,3,5]
[2,4,5]
[3,4,5]

If you want to get a bit-vector out of c then you can do e.g.

function combination_to_bitvector(l, c)
    result = zeros(l)
    result[c] = 1
    result
end

where l is the desired length of the bit-vector.