Fully Unique Combinations from Multiple Arrays

General Usage
, ,
1

Hello!

I have a model with n terms and the ability to remove k of them. However each of the k removal options can only act on an arbitrary list of model terms. I want a memory-efficient way to generate every unique combination without any repetitions. For example:

Given vectors a = [1;2;3] and b=[3;4], I would like to get the output of:
[1, 3] [1,4] [2, 3] [2, 4] [3, 4]

I have found that filter!(allunique,unique,(Iterators.product(a,b))) does produce the desired result but I’m worried about the memory consumption since my actual problem involves ~40 terms and 12 removals. Combinatorics.combinations(unique([a;b]),2) gets close but includes [1, 2] in the result. Is there a way to get fitler!() to work the Combinations or skip generating such results in the first place?

Thank you in advance

2

Hi,

You could probably use a simple for loop (equivalent), like

julia> [(x, y) for y in b, x in a if x != y]
5-element Vector{Tuple{Int64, Int64}}:
 (1, 3)
 (1, 4)
 (2, 3)
 (2, 4)
 (3, 4)

or more generally

julia> as = (a, b);

julia> [x for x in Iterators.product(as...) if allunique(x)]  # (Presumably you could add a couple of (Iterators.)reverse to get the same order as above)
5-element Vector{Tuple{Int64, Int64}}:
 (1, 3)
 (2, 3)
 (1, 4)
 (2, 4)
 (3, 4)

which doesn’t first collect any of the product.

julia> @btime filter!(allunique, unique(Iterators.product($a, $b)));
  204.437 ns (7 allocations: 752 bytes)

julia> @btime [x for x in Iterators.product($as...) if allunique(x)];
  58.460 ns (3 allocations: 224 bytes)

Edit: Or use collect(Iterators.filter(...)) instead of filter!(..., unique(...)):

julia> @btime collect(Iterators.filter(allunique, Iterators.product($as...)))
  46.216 ns (2 allocations: 192 bytes)
5-element Vector{Tuple{Int64, Int64}}:
 (1, 3)
 (2, 3)
 (1, 4)
 (2, 4)
 (3, 4)

(or don’t collect if you don’t have to).

Edit2: I was assuming here that the entries of as (a and b) are already internally unique (all(allunique, as)). If so, then also allunique(Iterators.product(as...)) is satisfied, so unique(Iterators.product(as...)) is equivalent to collect(Iterators.product(as...)), but slower. If not, then uniquifying the as (uas = unique.(as)) will be more efficient than the same for the product (unique(Iterators.product(as...))).

3

You can also take a look at the package SmallCombinatorics.jl.

4

I’m not sure that I fully understand the problem. How many lists do you have and how many elements does each list have roughly? And does the order of the elements in the output lists matter? For instance, what output would you get in your example above if b is the same list as a?

EDIT: I now think that given n collections c_1 to c_n you are looking for all sets of size n in the collection of all sets that are of the form \{x_1,\dots,x_n\} with x_1\in c_1, …, x_n\in c_n.

5

So I have 12 arrays of lengths ranging from 5 to 17 filled with integers ranging from 1 to 41. Each array has only unique elements but some of the elements are shared between arrays. I would like each of the outputs to get 1 element from each array and if there are two elements that only occur in a single array, I don’t want them together. I do not believe order will matter.

The output arrays will be used as indices for solving different sets of linear systems of equations and I’m trying to wasting computational time on combinations that won’t be solvable. I appreciate everyone’s thoughts thus far and hope this helps clarify what I’m hoping for

6

You should write down simple code that does what you want, albeit slowly, together with small examples, input and output copy-pasted from the REPL.

Then, I’d like to know the output size you desire. (if your output happens to be ~1e15 many combinations, then a memory efficient iteration won’t help you.)

Your explanation is extremely unclear and your purported “definition by implementation” does not work:

julia> a = [1;2;3]; b=[3;4]; filter!(allunique,unique,(Iterators.product(a,b)))
ERROR: MethodError: no method matching filter!(::typeof(allunique), ::typeof(unique), ::Base.Iterators.ProductIterator{Tuple{Vector{Int64}, Vector{Int64}}})
The function `filter!` exists, but no method is defined for this combination of argument types.

An example that I want to see included is a=[1,2];b=[1,2];: Do you want [[1,2], [2,1]] or do you want [[1,2]]?

7

I apologize, I don’t currently have code that works for my goal, even slowly and for the typo in my initial post. I had meant to type:

julia> a = [1;2;3]; b=[3;4]; filter!(allunique,unique(Iterators.product(a,b)))

However, I’ve realized even this breaks if I introduce c=[1;2;4].
I would want something like the following:

julia> a = [1;2;3]; b=[3;4]; c=[1;2;4]; f(a,b,c)
4-element Vector{Tuple{Int64, Int64, Int64}}:
 (1, 2, 3)
 (1, 3, 4)
 (1, 2, 4)
 (2, 3, 4)
julia> a = [1;2]; b=[1;2]; f(a,b)
1-element Vector{Tuple{Int64,Int64}}:
 (1, 2)

Just using combinations, there are more than 30e6 but my hope is that more robust combinatorics and cleverness can reduce this number. These are the actual arrays for my problem:

a = [3,4,6];
b = [19,20,39,40];
c = [12,17,23];
d = [13,15,16,17,21,24,25];
e = [12,13,15,16,17,25,26];
f = [29,30,38,41];
g = [13,14,15,16,17];
h = [23,24,26,29,41];
i = [23,24,25,29,30];
j = [32,34,35,36,41];
k = [31,32,33,34,35,40,41];
l = [31,32,33,34,36,39,40,41];

Thank you for the continued help and thoughts.

8

Do you mean something functionally equivalent to

julia> unique(([sort(SVector(x)) for x in Iterators.product([1;2;3], [3;4], [1;2;4]) if allunique(x)]))
4-element Vector{SVector{3, Int64}}:
 [1, 2, 3]
 [1, 2, 4]
 [1, 3, 4]
 [2, 3, 4]

julia> unique(([sort(SVector(x)) for x in Iterators.product([1;2], [1;2]) if allunique(x)]))
1-element Vector{SVector{2, Int64}}:
 [1, 2]

(but more (memory) efficient)?

9

Yes! Thank you! This produces the desired results for subsets of my problem but is too expensive for my full problem

10

This is not particularly efficient in terms of runtime since we still first generate all possible combinations, and then only keep the valid ones, but for memory consumption this should be fine:

using StaticArrays

function generate(as)
    N = length(as)
    T = eltype(eltype(as))
    s = Set{SVector{N, T}}()
    _generate!(s, Iterators.product(as...))  
    # Function boundary for type stability 
    # (We will have to recompile for every value of N though, but this is relatively fast)
   return s
end

function _generate!(s::Set{SVector{N, T}}, c) where {N, T}
    for x in c
        allunique(x) || continue
        sx = sort(SVector{N, T}(x))
        push!(s, sx)
    end
    return s
end

julia> as = [a, b, c, d, e, f, g, h, i, j, k, l]  # (or tuple (of (S)Vectors), or ...)
12-element Vector{Vector{Int64}}:
 [3, 4, 6]
 [19, 20, 39, 40]
 [12, 17, 23]
 [13, 15, 16, 17, 21, 24, 25]
 [12, 13, 15, 16, 17, 25, 26]
 [29, 30, 38, 41]
 [13, 14, 15, 16, 17]
 [23, 24, 26, 29, 41]
 [23, 24, 25, 29, 30]
 [32, 34, 35, 36, 41]
 [31, 32, 33, 34, 35, 40, 41]
 [31, 32, 33, 34, 36, 39, 40, 41]

julia> @time generate(as)
  7.380655 seconds (2.21 M allocations: 871.420 MiB, 0.85% gc time, 7.67% compilation time)
Set{SVector{12, Int64}} with 2707644 elements:
  [6, 13, 17, 19, 21, 25, 26, 32, 36, 38, 40, 41]
  [4, 14, 15, 23, 24, 26, 29, 30, 33, 36, 40, 41]
  [6, 15, 20, 23, 25, 26, 29, 30, 31, 33, 35, 38]
  [4, 14, 16, 21, 23, 24, 25, 29, 34, 35, 36, 40]
  [3, 13, 14, 17, 20, 24, 25, 26, 29, 34, 39, 40]
  [4, 12, 16, 17, 20, 21, 23, 30, 32, 36, 38, 41]
  [6, 12, 17, 19, 23, 24, 29, 30, 35, 38, 39, 41]
  [4, 12, 13, 16, 24, 26, 29, 34, 36, 38, 39, 40]
  [6, 12, 14, 16, 25, 26, 29, 30, 31, 35, 39, 41]
  [3, 15, 23, 24, 25, 30, 33, 35, 36, 38, 39, 41]
  [6, 12, 13, 15, 21, 23, 30, 34, 35, 38, 39, 40]
  [4, 13, 14, 15, 17, 19, 23, 25, 32, 36, 38, 40]
  [4, 12, 14, 16, 17, 20, 26, 30, 31, 33, 38, 41]
  [3, 12, 13, 23, 25, 26, 31, 32, 38, 39, 40, 41]
  [3, 13, 15, 20, 21, 23, 25, 26, 31, 33, 35, 41]
  [6, 12, 13, 15, 17, 26, 30, 33, 34, 35, 39, 41]
  [4, 12, 16, 17, 20, 23, 25, 26, 29, 32, 33, 34]
  [3, 14, 15, 19, 21, 23, 25, 32, 33, 34, 38, 41]
  [6, 13, 14, 17, 20, 21, 23, 24, 29, 31, 35, 36]
  [4, 16, 17, 21, 24, 25, 32, 35, 36, 38, 40, 41]
  [3, 12, 16, 17, 20, 23, 24, 29, 33, 34, 38, 41]
  [3, 14, 17, 21, 23, 24, 25, 31, 34, 35, 38, 40]
  [4, 13, 14, 17, 21, 25, 29, 32, 36, 38, 39, 40]
  [6, 12, 14, 16, 17, 20, 25, 29, 34, 35, 38, 41]
  [6, 14, 17, 19, 21, 23, 26, 29, 31, 35, 38, 41]
  ⋮

julia> @time s = generate(as);  # Compilation only adds many small allocations
  6.904371 seconds (68 allocations: 759.835 MiB, 0.95% gc time)

julia> (length(s) * sizeof(eltype(s))) / 2^20  # Minimal size of output, in MiB
247.8922119140625

julia> Base.summarysize(s) / 2^20  # Including Set overhead
388.00011444091797

So we are (in total) allocating about three times what we actually want to output.

(Edit: if you assume that the entries of a, b, … are small (cf. foobar_lv2 and matthias314’s solutions below, where the entries are assumed to lie between 1 and 64), you can of course also just use smaller containers than Vector{Int}s to reduce memory:

julia> as8 = Vector{Int8}.(as);  # (or UInt8, but Int8 prints nicer)

julia> @btime generate($as8);
  5.913 s (67 allocations: 101.84 MiB)

The bit mask approach below is more efficient though.)

1 Like
11

I think you should should use an appropriate representation (one such combination is a 64 bit number), and then store the thing and reuse it.

For example:

julia> items = [a,b,c,d,e,f,g,h,i,j,k,l]; 

julia> function gen!(res, abc, idx, current)
       if idx > length(abc) 
           push!(res, current)
           return
       end
       choices = abc[idx]
       for other in choices
           (current & 1<<other) != 0 && continue
           gen!(res, abc, idx+1, current | (1<<other))
       end
       end
julia> function gen(abc)
       res = Set{UInt64}()
       current = UInt(0)
       idx = 1
       gen!(res, abc, idx, current)
       res2 = collect(res)
       sort!(res2)
       res2
       end
julia> sizeof(gen(items))
21661152
julia> @time gen(items);
  0.935008 seconds (64 allocations: 111.821 MiB, 9.18% gc time)

Thats mere 21 megabytes, and runs in < 1 second for me, with ~ 2.7e6 elements. Just keep the resulting thing in a global const.

3 Likes
12

That’s what the type SmallBitSet from SmallCollections.jl does. Everything looks like sets, and the bit manipulations are under the hood. This may make the code more readable.

using SmallCollections

const S = SmallBitSet{UInt64}

Base.hash(s::S, h::UInt) = fasthash(s, h)  # faster hash (not compatible with `AbstractSet`)

function gen!(res, abc, idx, current)
    if idx > length(abc)
        push!(res, current)
        return
    end
    choices = abc[idx]
    for other in choices
        new = push(current, other)  # add element to set
        new != current && gen!(res, abc, idx+1, new)
    end
end

function gen(abc)
     res = Set{S}()
     gen!(res, abc, 1, S())
     res
end

items = S.([a,b,c,d,e,f,g,h,i,j,k,l])

(I skip the conversion from Set to Vector.) The run time is as before:

julia> @b gen($items)  # foobar_lv2 (with result as a set)
1.453 s (55 allocs: 70.509 MiB, 0.34% gc time, without a warmup)

julia> @b gen($items)  # this version
1.411 s (55 allocs: 70.505 MiB, 0.29% gc time, without a warmup)
5 Likes
13

That’s up to the context of whatever OP wants to do with the result. I wrote the proposal that way for 2 reasons:

  1. The thing to optimize for is probably memory consumption of the final object and iteration speed of the final object. (I would hope that OP generates the thing once, and uses it often)
  2. The sorting of intermediate results is imo good practice when affordable in order to get a canonical and reproducible output. This simplifies debugging and testing of alternative implementations. The iteration order of Set is kinda unstable otherwise (but it is reproducible at least, since there are no ref-types involved – for ref-types, hashes and hences iteration order depends on addresses which depend on malloc internal state, which is not reproducible unless you go to great lengths).
14

In case overlapping indices are somewhat common, this seems to perform better:
Edit: actually I think this is faster because it vectorizes.

function disjoint_combinations(A, B)
	res = Set{S}()
	# Base.sizehint!(res, length(A)*length(B)) # reduces allocation count but slightly increases runtime and allocated memory
	(length(A) == 0 || length(B) == 0) && return res
	nAB = length(first(A)) + length(first(B))
	for a in A
		for b in B
			ab = union(a,b)
			if length(ab) == nAB
				push!(res, ab)
			end
		end
	end
	return res
end

gen_abraemer(abc) = reduce(disjoint_combinations, abc)


item_sets = [S.(a),S.(b),S.(c),S.(d),S.(e),S.(f),S.(g),S.(h),S.(i),S.(j),S.(k),S.(l)]

I get

julia> @b gen_matthias($items)
649.157 ms (55 allocs: 70.501 MiB, without a warmup)
julia> @b gen_abraemer($item_sets)
236.072 ms (310 allocs: 116.693 MiB, without a warmup)
# with Base.sizehint!: 249.912 ms (62 allocs: 188.273 MiB, without a warmup)

There is more opportunity for optimization here: It matters in what order the sets are combined. If one merges sets that have more overlap early, then that reduces the number of combinations considered. Maybe one could make reasonably simple “greedy scheduler” that tries to take good route by merging the sets with the largest overlap first or something.
E.g. with the current example:

julia> @b gen_abraemer($(reverse(item_sets)))
148.111 ms (329 allocs: 124.012 MiB)

almost twice as fast!

2 Likes
15

I would first sort the sets to choose from lexicographically, then iterate over the cartesian product while maintaining a list of prefixes you have “exhausted.” then if you would start building a combination with some exhausted prefix you can cut off that entire subtree

16

I don’t think that sorting guarantees you any kind of ‘good’ ordering in general. The problem is invariant under renaming the constituents after all, so any order is in some sense arbitrary.

Here is an attempt at choosing a good order for the reduction:

function gen_greedy(sets)
	sets = Set.(sets)
	while length(sets) > 2
		ind = argmax(c -> count(in(sets[c[1]]), sets[c[2]]), 
					 ((a,b) for a in eachindex(sets) for b in eachindex(sets) if a < b))
		sets[ind[1]] = disjoint_combinations(sets[ind[1]], sets[ind[2]])
		deleteat!(sets, ind[2])
	end
	return disjoint_combinations(sets[1], sets[2])
end

With this I get:

julia> @b gen_greedy($item_sets)
182.278 ms (218 allocs: 71.289 MiB)

So faster and less memory allocated but overall slower then the reversed. So this probably not quite the right strategy for merging.

Edit: Kinda obvious in hindsight. This implementation does not compute the ‘overlaps’ between set of collection of different size correctly. So after merging the first 2 sets, the result has 0 overlap with all other sets because it consists of collections of 2 elements while the other set only have collections of size 1. I think we should be able to compute the overlaps between the sets in the beginning and then construct the reduction path solely based on that info but probably don’t have more time now to think about it.