Enumerating tuples from a larger space such that all pairs of values are present at least once

Question

I am trying to write a function that will accept a list of N lists and returns a list of N-tuples that represent all pairs of values from the original N lists at least once. To illustrate with an example, suppose we have:

foo = { { A, B }, { C, D, E, F }, { X, Y } };

Pairwise[ls_List] := (* some implementation here *)

Pairwise[foo] ==> { {A,C,X}, {A,D,X}, {A,E,X}, {A,F,X}, {A,C,Y}
                    {B,C,Y}, {B,D,Y}, {B,E,Y}, {B,F,Y}, {B,D,X} }

This particular enumeration ensures that all possible pairs from each of the lists appear at least once:

A-C, A-D, A-E, A-F, A-X, A-Y, B-C, B-D, B-E, B-F, B-X, B-Y, C-X, D-X, E-X, F-Y, C-Y, D-Y, E-Y, F-Y

Unlike an exhaustive enumeration, the pairwise set has substantially fewer members. The trouble is, it's not clear to me how to efficiently produce this kind of enumeration without keeping track of which pairs I've already produced. Since my input lists are likely to be quite large, it would be nice to use an algorithm that avoids tracking any of the intermediary outputs or have to test against an increasingly large set of intermediaries.

It's not necessary that the results contain each pair-wise combination only once - but it should definitely not simply yield the full Cartesian product.

Answer 1

More a comment than an answer, (but I have not enough reputation): This is a well known problem, although by far not solved. (One might not expect, but is of very practical relevance in software testing.) You find a lot of interesting stuff (theory and algorithms) by googling "orthogonal array" or - even better - "mixed orthogonal array". Also the book "Orthogonal arrays, theory and applications" by Hedayat, Sloane, and Stufken, contains lots of references and is referenced itself a lot in papers on this topic. ... and of course Wikipedia: http://en.wikipedia.org/wiki/Orthogonal_array. And yes he is the "Encyclopedia of Integer Sequences"-Sloane. For Orthogonal Arrays see his website http://neilsloane.com/doc/cent4.html

Edit: A generalization of the orthogonal arrays is called "Covering Arrays" (for the difference see e.g. http://csrc.nist.gov/groups/SNS/acts/covering-orthogonal-arrays.html). "Combinatorial Testing" is also very helpful term for Google searching.

Another nice website is at NIST: http://csrc.nist.gov/groups/SNS/acts/index.html, where you can download free software on this topic. For a scientific paper on covering arrays, where the problems and algorithms are discussed see e.g. https://www.research.ibm.com/haifa/projects/verification/mdt/papers/AlgorithmsForCoveringArraysPublication191203.pdf. There also the additional problem of restricted covering arrays is discussed, where you forbid certain pairs (or more general tuples) in the result, because they are not supported by the software you want to test.

Edit: In the terminology of Hartman et al. your problem is a $CA(2, 3, 2, 3, 2)$, a covering array of strength 2 with three domain sets. Two domain sets have size 2, one has size 3.

Edit: I recommend "Automatic Generation of Combinatorial Test Data (Springer Briefs in Computer Science)" by Zhang, Zhang & Ma (Springer, 2014), which very nicely summarizes all important algorithms for test data generation (including constrained problems). It contains a very good reference list, too.

Answer 2

I haven't yet tried to imagine an algorithm for your problem but here is a visualization that might help you approach it.

pairs[ls_List] := Join @@ Tuples /@ Subsets[ls, {2}];

plot[ls_List, set_List] :=
  ArrayPlot[
    Outer[SubsetQ, set, pairs @ ls, 1] // Boole,
    FrameTicks -> All,
    FrameLabel -> {"Tuples", "Pairs"}
  ]

If you are not using Mathematica 10 then also define:

SubsetQ[a_,b_] := a ⋂ b === Sort[b]

For example with your foo (using Strings rather than defined Symbols such as E):

foo = {{"A", "B"}, {"C", "D", "E", "F"}, {"X", "Y"}};

plot[foo, Tuples @ foo]

This shows for each of the 16 tuples which of the 20 pairs are represented. Reading left-to-right there are three distinct patterns; these correlate to the Subsets[ls, {2}] groups used by pairs.

As a second example here is your target set:

target =
 {{"A", "C", "X"}, {"A", "D", "X"}, {"A", "E", "X"}, {"A", "F", "X"},
  {"A", "C", "Y"}, {"B", "C", "Y"}, {"B", "D", "Y"}, {"B", "E", "Y"},
  {"B", "F", "Y"}, {"B", "D", "X"}};

plot[foo, target]

If I have time and interest I'll revisit this problem later but I hope for now this provides at least some small insight.

Answer 3

I came up with a little algorithm that produces minimal lists of paths. I might be overlooking something, so please let me know if you find any mistakes.

I shall focus on the case where the input list has 3 sublists, which I shall call layers henceforth. Cases with more than three layers can be computed recursively.

The algorithm is based on three ingredients:

1. The output is symmetric under permutations of the input layers.
2. Suppose we have two layers L and l with Length[L] > Length[l]. We can cluster elements in L together in such a way that the number of clusters is equal to Length[l].
3. We can iterate over the layers by rotating them. If we do this in a particular way, we can ensure that all nodes are connected using a minimal number of iterations (resulting in a minimal number of connecting paths).

It's probably best to show some examples rather than showing the algorithm in full generality (I included a link to my code below).

Example 1

Let's start with the example provided: foo = {{a, b}, {c, d, e, f}, {x, y}}.

• First, we use the symmetric property to get foo = {{a, b}, {x, y}, {c, d, e, f}}.
• Second, we partition the third layer in a way that is compatible with the second layer {x,y}, i.e. {{c, d}, {e, f}}. We thus obtain foo = {{a, b}, {x, y}, {{c, d}, {e, f}}}.
• This partitioning allows us to define a tree that connects the first element in the first layer to all nodes in the third layer:

• Finally, we 'rotate' the first and second layer as follows:

This second configuration can be obtained by rotating the first and second layer, thus computing the tree associated with {{b, a}, {y, x}, {{c, d}, {e, f}}}.

• The output is given by all 8 paths drawn in the two pictures above:

{{a, c, x}, {a, d, x}, {a, e, y}, {a, f, y}, {b, c, y}, {b, d, y}, {b, e, x}, {b, f, x}}

Example 2

In this example, we assume that the pre-processing steps have been performed, such that: foo = {{a, b}, {c, d, e, f}, {X, Y, Z, W}}, where e.g. X is a cluster (i.e. a list itself).

The tree with its four permutations are then given by:

Notice that root nodes are guaranteed to connect to all others by construction. The trick is to get the second-layer nodes to connect to all third-layer nodes without doing unnecessary work. We avoid doing this unnecessary work because we cluster the rotations of the second layer in a way compatible with the first layer. This is similar to the clustering of the third layer w.r.t. the second layer in the preprocessing step.

For completeness, the output thus generated is:

{{a, c, X}, {a, d, Y}, {a, e, Z}, {a, f, W}, {a, d, X}, {a, e, Y}, {a, f, Z}, {a, c, W}, {b, e, X}, {b, f, Y}, {b, c, Z}, {b, d, W}, {b, f, X}, {b, c, Y}, {b, d, Z}, {b, e, W}}

Implemented this in full generality, but the code is bit long. You can check it out here.

Having said all this, I do really hope that there's a simpler (or at least more mathematica-native) way to do this.

P.S. I also analyzed performance. Below you can see a plot of the computing time vs. input size. The straight line is quadratic fit, so the time approximately scales as $T\sim n^2$. Moreover, the output size of this algorithm is given by the product of the sizes of the two largest layers.

Answer 4

This feels quite a bit like an odd mix of systems of distinct representatives and block designs, although this exact problem isn't coming to mind as a particular construction in any of these combinatorial contexts. It would probably help to pull out a book about matroids too -- that's not my forte either.

It's important to note that all of your sets will be transversals of your set system, at least, that's how I'm reading the question. Otherwise, the union of these sets will trivially contain all the pairs... I'm quite sure that's not what you want!

I think it will be easier if we assume your set system is as follows:

L[1]={x[1,1], x[1,2], x[1,3], ... , x[1,n[1]]}
L[2]={x[2,1], x[2,2], x[2,3], ... , x[2,n[2]]}
...
L[k]={x[k,1], x[k,2], x[k,3], ... , x[k,n[k]]}

It is important to note, then, that such a configuration depends only on ${n[1],...,n[k]}$. If you have some other lists that you want to apply this to, that's fine, you can just pick of the x[i,j] in each transversal as representing the $j^{\mathrm{th}}$ element of the $i^{\mathrm{th}}$ list.

We also may as well assume each n[i] is at least 2. If it's just 1, then add x[1,n[i]] to every transversal in the set system minus those. We can also assume the n[i] are increasing.

My approach to the problem will be induction

First of all, what if $k=1$? Then you just return the set of single-element sets:

{{x[1,1]},{x[1,2]},...,{x[1,n[1]]}}

Now, if $k=2$ you literally need every combination, don't you?

I don't have time right now to elaborate or post much code, but if you let T be a set of transverals that covers your system, and say T-i is that set of transversals, each missing the element from set L[i].

Now, it's easy enough to get a set T of size n[k]. Let S[i] be each set of T-i with every possible combination of elements from L[i]. That's at most n[k]^2 sets, and the union of all of these S[i] gives you what you want with at most k*n[k]^2 sets. In fact, you can skip the last one, giving you (k-1)*n[k]^2, although that's not really an improvement -- linear is linear.

As I mentioned, I don't really have much time to give this more thought, it's interesting, but (unless I've done something wrong in my rush to say something useful in the time I have for MMASE today) it's quadratic in n[k] and linear in k, worst case.

This would produce, for the given example:

S = { { A, B }, { C, D, E, F }, { X, Y } }
T = { { A, C, X }, { B, D, Y }, { A, E, X }, { B, F, Y } }

T-1 = { { C, X }, { D, Y }, { E, X }, { F, Y } }
T-2 = { { A, X }, { B, Y }, { A, X }, { B, Y } }
T-3 = { { A, C }, { B, D }, { A, E }, { B, F } }

S[1] = { { A, C, X }, { B, C, X },
         { A, D, Y }, { B, D, Y },
         { A, E, X }, { B, E, X },
         { A, F, Y }, { B, F, Y } }
S[2] = { { A, C, X }, { A, D, X }, { A, E, X }, { A, F, X }, 
         { B, C, Y }, { B, D, Y }, { B, E, Y }, { B, F, Y }, 
         { A, C, X }, { A, D, X }, { A, E, X }, { A, F, X }, 
         { B, C, Y }, { B, D, Y }, { B, E, Y }, { B, F, Y } }

Now in this case, you get way too many -- and in general, this is not an especially efficient way of doing it, but I'm fairly sure this bound is better than n1*n2*n3*...*nk for whatever purposes. It wouldn't be particularly hard to code, but I don't have the time.

This could probably be faster/more efficient by making the judicious choice of a different covering set of transversals for each i. That and throwing away duplicates (throwing away duplicates, by the way, seems to give you 12 in this case).

This is assuming no use of anything even remotely intelligent like deleting duplicates in T-2. This could be optimized and I'm not of any particular belief that the bound k*n[k]^2 is tight in any way. Consider this answer an extended comment more than an answer, i just want to get this bound out there for everyone to chew on.

Belisarius is asking if there are m (my k) groups of 2, what is it, and in this case you might get lazy and call them {x1,x1} through {xk,yk}. The scheme I propose, with duplicates knocked out, I think is going to give you the following:

{x1,x2,x3,...,x(k-1),xk}
{y1,y2,y3,...,y(k-1),yk}
{y1,x2,x3,...,x(k-1),xk}
{x1,y2,y3,...,y(k-1),yk}
{x1,y2,x3,...,x(k-1),xk}
{y1,x2,y3,...,y(k-1),yk}
{x1,x2,y3,...,x(k-1),xk}
{y1,y2,x3,...,x(k-1),yk}
...
...
{x1,x2,x3,...,y(k-1),xk}
{y1,y2,y3,...,x(k-1),yk}

Obviously this list is incredibly redundant! Probably because all the n[i] are 2, it is exaggeratedly inefficient.

Answer 5

Update

I've suceeded to find some solutions to particular n-lists (n>3) cases and also I have found a ridiculous short algorithm (2 lines of code) which seems to give the optimal n-tuples lists (==> all the possible pairs are present only once) for a large class of configurations where the input n-lists have all the same number of elements.(And I think there is some room for improvement.) For now, the short automated tests I made gave me these results:

(*{length of n-lists, number of elements in each list,
{{max number of tuples, number of n-tuples in the solution},
number of pairs},yes}*)
{3,3,{{27,9},27},yes}
{3,5,{{125,25},75},yes}
{3,7,{{343,49},147},yes}
{3,9,{{729,81},243},yes}
{3,11,{{1331,121},363},yes}
{3,13,{{2197,169},507},yes}
{3,15,{{3375,225},675},yes}
{3,17,{{4913,289},867},yes}
{3,19,{{6859,361},1083},yes}
{4,5,{{625,25},150},yes}
{4,7,{{2401,49},294},yes}
{4,11,{{14641,121},726},yes}
{4,13,{{28561,169},1014},yes}
{4,17,{{83521,289},1734},yes}
{4,19,{{130321,361},2166},yes}
{5,5,{{3125,25},250},yes}
{5,7,{{16807,49},490},yes}
{5,11,{{161051,121},1210},yes}
{5,13,{{371293,169},1690},yes}
{5,17,{{1419857,289},2890},yes}
{5,19,{{2476099,361},3610},yes}
...

For example the last line means that I found an optimal solution for an input list of 5 lists, each containing 19 elements. The corresponding maximimum number of 5-tuples lists for this config is 2476099, the number of 5-tuples given by the solution is 361 where all the possible 3610 pairs are given only once.

If the OP or someone is interested in the algo. and test code I can publish it here later. (Actually it is rather clear that it is not a Mathematica related post anymore.)

---

Previous

I've been playing with this problem during some time now, and what I observe is that indeed one can optimise the solution at different nested levels ... (it is hard to explain and I will give an example later to illustrate that). What does that mean is that for a given problem you can give different solutions more or less optimised (where optimisation consists in giving the less number of N-tuples list). For the given problem : { { A, B }, { C, D, E, F }, { X, Y } }, there are a maximum of 16 3-lists, the OP gave a solution with 10 lists, the optimal is 8 lists.

Anyway, for now here is some code to produce optimal solution for 3-lists. (Sorry, the code is very ugly for now, but I focused on the solution.)

q2[x_, y_] := Table[Thread[{x, RotateLeft[y, inc]}], {inc, 0, Length@x - 1}];

mk3[{lis1_, lis2_, lis3_}] := Module[
{lenlis1 = Length@lis1, lenlis2 = Length@lis2, lenlis3 = Length@lis3, list12,
templist1, templist2, lenlist12},
list12 = If[lenlis1 == lenlis2,q2[lis1, lis2],
(templist1 = q2[lis1, Take[lis2, lenlis1]];
 templist2 = Thread /@ Thread[{lis1,Take[lis2, -(lenlis2 - lenlis1)]},List, -1];
 MapAt[Sequence @@ # &, Thread[{templist1, Sequence @@ templist2}], {All, 1}]
 )
];
lenlist12 = Length@list12;
If[lenlist12 == lenlis3,
MapAt[Sequence @@ # &, 
Thread /@ Thread[{list12, lis3}], {All, All, 1}],
Join[MapAt[Sequence @@ # &, 
 Thread /@ Thread[{list12, Take[lis3, lenlist12]}], {All, All, 1}],
MapAt[Sequence @@ # &, 
 Thread /@ 
  Thread[{list12[[1]], Take[lis3, -(lenlis3 - lenlist12)]}, 
   List, -1], {All, All, 1}]]]]

When the input is 3 lists with the same number of elements, then the code is a one-liner :

mk3[{s1_, s2_, s3_} /; Equal @@ (Length /@ {s1, s2, s3})] := Table[
Join @@@ Thread[{Thread[{s1, RotateLeft[s2, inc]}],
{First@RotateLeft[s3, inc]}}, List, 1],
{inc, 0, Length@s1 - 1}]

Applications :

sol = mk3[{{a, b}, {x, y}, {c, d, e, f}}]
TableForm /@ sol // TableForm

gives

sol = mk3[{{a, b}, {w, x, y, z}, {m, n, o, p}}];
TableForm /@ sol // TableForm

gives

sol = mk3[{{a, b, c, d}, {w, x, y, z}, {m, n, o, p}}];
TableForm /@ sol // TableForm

gives

You can try of course other combinations, the only constraint for the input in mk3is that the lists are sorted by increased length size.

I will edit this post later to give some examples of other optimised solutions but for longer input N-lists (N>3) and with more detailed explanations.

Answer 6

That's a possible solution. I was working exactly on the same problem (nice combination!!!) just last week, but in my case N-tuples have length 6 and not 2 as in your case. However, here I'll post the excerpt of my code for the case of your interest. Just for curiosity: the solution I found has been achieved after many failures. I started thinking about what was required as output, but no solution came out, then I moved to think to the "dual" problem, that's to say: what I don't want in the output. Many times, in finding a computational solution it helps to think in the opposite way: do not look for what you want, try to discard what you don't want. When you find it just have to make the Complement with all possible solutions.

The problem can be simplified as follows: given an alphabet of symbols (or numbers) and a set of rules (subsets of symbols, with different length), we have to enumerate all pairs with elements from the alphabet, discarding those indicated by the rules. I call "rules" the sub-lists on the original list, {{A, B}, {C, D, E, F}, {X, Y}} in your example.

These sub-lists are three rules which I interpret as follows:

Each rule says the element that cannot be in the same pair. Indeed, in the solution you don't have {A, B}, {C, D}, {C, E}, ... Once we have re-written the problem, we can go ahead with a possible solution.

I calculate all possible pairs from the first rule, then from the second rule and so I do the Intersection of pairs from the first two rules. So I obtain all pairs that satisfy the first and the second rule. Then, I go on with the third rule and apply the Intersection with the pairs from the previous step. Going ahead till the last rule, I get all pairs that satisfy, at the same time, all rules. Honestly, I'm not happy with the code I wrote for this solution, I think it's not elegant as usually the Mathematica code is, but I had no time to improve the elegance of the code, even because the calculation time and the memory required are not so bad. Here is the code. Note: I use numbers because some capital letters are system' symbols (C, D and E) and I don't like to work with strings.

alphabet = {1, 2, 3, 4, 5, 6, 7, 8};
rules = {{1, 2}, {3, 4, 5, 6}, {7, 8}};
With[{A = First[rules]},
  res = Join[Subsets[Complement[alphabet, A], {2}],
    Flatten[
     Map[Flatten[
        Outer[Sort[Join[#1, #2]] &, {#}, 
         Partition[Complement[alphabet, A], 1], 1], 1] &, 
      Partition[alphabet, 1]], 1]]];
Map[
  Function[A,
   res = Intersection[res, Join[Subsets[Complement[alphabet, A], {2}],
      Flatten[
       Map[Flatten[
          Outer[Sort[Join[#1, #2]] &, {#}, 
           Partition[Complement[alphabet, A], 1], 1], 1] &, 
        Partition[alphabet, 1]], 1]]]], Rest[rules]];

alphabet

{1, 2, 3, 4, 5, 6, 7, 8}

rules

{{1, 2}, {3, 4, 5, 6}, {7, 8}}

res

{{1, 3}, {1, 4}, {1, 5}, {1, 6}, {1, 7}, {1, 8}, {2, 3}, {2, 4}, {2, 5}, {2, 6}, {2, 7}, {2, 8}, {3, 7}, {3, 8}, {4, 7}, {4, 8}, {5, 7}, {5, 8}, {6, 7}, {6, 8}}

ReplaceAll[res, {1 -> "A", 2 -> "B", 3 -> "C", 4 -> "D", 5 -> "E", 6 -> "F", 7 -> "X", 8 -> "Y"}]

{{"A", "C"}, {"A", "D"}, {"A", "E"}, {"A", "F"}, {"A", "X"}, {"A", "Y"}, {"B", "C"}, {"B", "D"}, {"B", "E"}, {"B", "F"}, {"B", "X"}, {"B", "Y"}, {"C", "X"}, {"C", "Y"}, {"D", "X"}, {"D", "Y"}, {"E", "X"}, {"E", "Y"}, {"F", "X"}, {"F", "Y"}}

Here is a test with a longer alphabet and more rules

alphabet = Range[200];
rules = Table[RandomChoice[alphabet, RandomInteger[{2, 10}]], {100}];
AbsoluteTiming[
 With[{A = First[rules]},
  res = Join[Subsets[Complement[alphabet, A], {2}],
    Flatten[
     Map[Flatten[
        Outer[Sort[Join[#1, #2]] &, {#}, 
         Partition[Complement[alphabet, A], 1], 1], 1] &, 
      Partition[alphabet, 1]], 1]]];
 Map[
  Function[A,
   res = Intersection[res, Join[Subsets[Complement[alphabet, A], {2}],
      Flatten[
       Map[Flatten[
          Outer[Sort[Join[#1, #2]] &, {#}, 
           Partition[Complement[alphabet, A], 1], 1], 1] &, 
        Partition[alphabet, 1]], 1]]]], Rest[rules]];]

{1.872003, Null}

If we want to check the result

AbsoluteTiming[Flatten[MapThread[Cases[rules, {___, #1, ___, #2, ___}] &, Transpose[res]]]]

{0.872001, {}}

no one pair has two numbers mentioned in a rule. Note that for bigger alphabets the check requires much more time than the calculation. However, it has been useful the first time to be sure the result is exactly what we are looking for. Hopefully this is a suitable solution for your problem.

Answer 7

Here is a randomized algorithm that computes a solution:

Pairwise[ls_List]:=Module[{tuples, cumulativeNumberOfUniquePairs},
    tuples = RandomSample@Tuples[ls];
    cumulativeNumberOfUniquePairs = Length /@ FoldList[Union[#1,Subsets[#2,{2}]]&, {}, tuples];
    Pick[tuples, Differences[cumulativeNumberOfUniquePairs], _?Positive]
]

The algorithm starts from a shuffled list of all possible tuples of the input list ls. It then computes the number of pairs that each tuple contributes to the list of unique tuples and then only picks those tuples that actually have a positive contribution.

Because it is a randomized algorithm, it only produces a good result if you are lucky:

In[1]:= Pairwise[{{"A", "B"}, {"C", "D", "E", "F"}, {"X", "Y"}}]
Out[2]= {{B,C,Y},{A,E,Y},{B,E,Y},{B,E,X},{A,D,X},{A,C,X},{B,F,Y},{B,D,Y},{A,F,X}}

Answer 8

Here is a LinearProgramming approach:

pairCover[sets_] := Module[
    {t = Tuples[sets], pairs=Flatten[Tuples/@Subsets[sets,{2}], 1], lp},

    lp = Quiet[
        LinearProgramming[
            ConstantArray[1, Length[t]],
            Table[Boole[ContainsAll[pairs[[k]]] /@ t], {k, Length[pairs]}],
            Table[{1, 1}, {Length[pairs]}],
            0,
            Integers
        ],
        LinearProgramming::lpip
    ];
    Pick[t, lp, 1]
]

The vector that is minimized is the bit vector of the possible tuples. The constraints are that each pair occurs at least once.

A couple examples:

pairCover[{{a, b}, {c, d, e, f}, {x, y}}]
pairCover[{{a, b}, {c, d, e, f}, {x, y, z, t}}]

{{a, c, y}, {a, d, y}, {a, e, y}, {a, f, x}, {b, c, x}, {b, d, x}, {b, e, x}, {b, f, y}}

{{a, c, x}, {a, d, y}, {a, d, z}, {a, d, t}, {a, e, x}, {a, e, z}, {a, e, t}, {a, f, y}, {a, f, z}, {b, c, y}, {b, c, z}, {b, c, t}, {b, d, x}, {b, e, y}, {b, f, x}, {b, f, t}}