Utility Corner: Combinations
Vector 18:1, page 31
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 31
Here is another Zark newsletter extract. We were running short of issues to work from, but many thanks must go to Gary Bergquist who has once again supplied us with some more copies so that we can continue this series of reprints for some time to come. Below we reprint the solutions to the COMBIN problem, which throws up some immense variety, particularly in the resulting timings. Note: Jonathan Barman points out to me that whilst writing this article up he came across a dynamic function that is shipped with Dyalog APL version 9.0 (cmat is shipped in the dfns workspace), which performs the same task. Jonathan found that it was about 4 times slower than the fastest solution presented below. (It is however faster than the other solutions.) However, this proves once again that there is always worth in revisiting old problems with new technology at hand. There may be new, more elegant ways of solving an old problem, which may also turn out to be faster. Please don’t be afraid to try and shed new light on the problems being discussed in this column.
Utility Corner: Combinations¶
(The purpose of this column is to make you more productive by introducing you to utility functions. Think of utility functions as APL functions that have names instead of symbols. By expanding your function vocabulary, you’ll be able to write APL code that’s more concise, more efficient, and more readable.)
In last issue’s LIMBERING UP column, we asked you to write a function named COMBIN that returns a matrix of all combinations defined by its arguments. For example,
In this example, COMBIN returns all combinations of 5 items taken 3 at a time. In general, for M COMBIN N, the result is an integer matrix with M!N rows and M columns, in which each row is a distinct combination of the numbers ⍳N.
What makes this problem interesting is its hidden complexity. While simple to explain, it’s not so simple to solve.
We received a wide variety of solutions. For this problem, there seems to be an inverse relationship between simplicity and efficiency. The simplest functions were slowest. The fastest functions were obscure. Take your pick: clarity or speed.
(Note: In the following presentation, we look liberties. The functions received have been changed to make them work in either origin, to be consistent with one another, and in some cases to be a bit faster or a bit clearer [we think].)
The simplest approach was to recurse. Here’s a recursive nonlooping solution from Howard Peelle:
∇ R←M COMBIN N
[1] ⍝ Doubly recursive algorithm.
[2] →(M=1,N)/L1,L2
[3] R←⎕IO,1+(M-1)COMBIN N-1
[4] R←R,[⎕IO]1+M COMBIN N-1
[5] →0
[6] L1:R←(N,1)⍴⍳N ⍝ 1!N
[7] →0
[8] L2:R←(1,N)⍴⍳N ⍝ N!N
∇
Maria Wells sent us another recursive solution that includes a bit of nesting logic to keep you on your toes:
∇ R←M COMBIN N;A;P
[1] ⍝ Recursive nested algorithm:
[2] →(M=1,N)/L1,L2
[3] A←+\(P←1+N-M)⍴1
[4] R←⊃,[⎕IO]/(⍳P),¨A+(M-1)COMBIN¨N-A
[5] →0
[6] L1:R←(N,1)⍴⍳N ⍝ 1!N
[7] →0
[8] L2:R←(1,N)⍴⍳N ⍝ N!N
∇
As mentioned above, these simple functions paid the price in efficiency. Howard’s function was 14 times slower and Maria’s 10 time slower than the fastest function submitted. (All functions were timed on an APL*PLUS II system for the expression 4 COMBIN 15, which returns 1365 rows. Timings vary for different arguments and different APL systems.)
The next function, also submitted by Howard Peelle, uses a "truth table" algorithm. Here’s the approach: each row in the result of COMBIN can be viewed as the result of the expression T/⍳⍴T, where T is an N-element bit vector with M 1s. Find all the different combinations of N bits (there are 2*N of them); select just those with M 1s; and use the resulting bit vectors to compress the vector ⍳N.
∇ R←M COMBIN N
[1] ⍝ Truth table algorithm.
[2] ⍝ All combinations of N 0s and 1s:
[3] R←(N⍴2)⊤(⍳2*N)-⎕IO
[4] ⍝ Pick those with exactly M 1s:
[5] R←,⍉⌽(M=+⌿R)/R
[6] R←((M!N),M)⍴R/(⍴R)⍴⍳N
∇
Because of the 2*N component in this function, it slows down dramatically as N gets larger. It’s also prone to WS FULLs. For 4 COMBIN 15, the function was 92 times slower than the fastest function received. (When a function gets that slow, you find yourself rooting for the WS FULL.)
Another "select-from-all-combinations" approach, which was used in a function from Craig Tien, is suggested by viewing the result of COMBIN after subtracting ⍳M from each row. For example,
The result is a table of base-3 digits (or, in general, base-P digits where P←1+N-M). Translate the vector ⍳P*M (in origin 0) into its base-P digits. Select just those sets of digits that are in ascending order (eliminating, for example, 0 1 0, 0 2 0, and 0 2 1). Add each selected set to ⍳M.
∇ R←M COMBIN N;P
[1] ⍝ The nonlooping counting method.
[2] P←1+N-M
[3] ⍝ Create a table with elts between 0
[4] ⍝ and P-1 to add to vector ⍳M.
[5] R←(M⍴P)⊤(⍳P*M)-⎕IO
[6] ⍝ Keep only those in ascending order:
[7] R←⍉(∧⌿(1 0↓R)≥¯1 0↓R)/R
[8] ⍝ Add vector ⍳M to each one:
[9] R←R+(⍴R)⍴⍳M
∇
This function, like the one before it, has a tendency to get slow and is WS FULL prone on large problems. However, it grows according to (1+N-M)*M rather than 2*N, so it’s not too bad for values of N and M that don’t differ much. For the 4 COMBIN 15 test, it was 17 times slower than the fastest function received.
Henri Brudzewsky used a similar approach but looped to avoid the WS FULL.
∇ R←M COMBIN N;A;I;II;J;MM;P
[1] ⍝ The looping counting method.
[2] MM←M!N ⍝ Number of combinations
[3] A←⍳M ⍝ First combination
[4] R←(MM,M)⍴A ⍝ Prepare the output mat
[5] I←⎕IO ⍝ Row counter
[6] P←M/1+N-M ⍝ Base values to increment
[7] J←P⊤II←0 ⍝ Unpacked counter
[8] MM←MM+⎕IO ⍝ Termination inded
[9] L:→(I≥MM)/0 ⍝ Exit when done
[10] R[I;]←A+J ⍝ E.g. 1 2 3 4+0 0 1 2
[11] II←II+1 ⍝ Next increment
[12] J←P⊤II ⍝ Unpack increment
[13] J←⌈\J ⍝ Eliminate previous values
[14] II←P⊥J ⍝ Re-adjust increment
[15] I←I+1 ⍝ Inded of next row
[16] →L
∇
On line 13, ⌈\ is used to quickly ignore all sets of bas-P numbers that are not in ascending order. While WS FULL-resistant, this function is somewhat slower that the previous one (for 4 COMBIN 15), running 27 times slower than the fastest function submitted.
The last two functions are the fastest and the most obscure. What they have in common is their lop-by-column approach. Instead of working on the rows of the result of COMBIN (as do the previous three solutions) or on the submatrices of the result (as do the two recursive solutions), they build one column at a time.
The function from Jack Zinn loops from the first column to the last. Jack noticed that each column in the result of COMBIN contains numbers from the set C+⍳P, where P←1+N-M and C is the column indeed (in origin 0). For example, from 3 COMBIN 5, the first column contains numbers from ⍳3, the second from 1+⍳3, and the third from 2+⍳3.
The first column from 3 COMBIN 5 is 6 3 1/⍳3. Notice that 6 3 1 is 2!4 3 2. The second column is 3 2 1 2 1 1/1+⍳3. The numbers 3 2 1 come from 1!3 2 1, and their pattern from the indices 1 2 3,2 3,3 or (0+⍳3),(1+⍳2), (2+⍳1). From these observations, and the power of induction (or luck), Jack developed the following function:
∇ R←M COMBIN N;C;I;IP;P
[1] ⍝ Working from first column, formula:
[2] ⍝ R←A[I]/B[I]
[3] ⍝ where A is say 3!5 4 3; B is say
[4] ⍝ 1+⍳3; I starts at ⍳P and then
[5] ⍝ becomes a function of the previous
[6] ⍝ element of I, where the segment
[7] ⍝ starts at that value of I and
[8] ⍝ continues to the end of ⍳P
[9] IP←(I←⍳P←1+N-M)-⎕IO ⍝ Origin 0 inds
[10] R←((M!N),M)⍴0 ⍝ Build all-0 result
[11] M←M-~C←⎕IO ⍝ Index of last iteration
[12] L1:R[;C]←((M-C)!⌽IP+M-C)[I]/(IP+C)[I]
[13] →(M<C←C+1)/0 ⍝ Exit if done
[14] I←⊃,/(I-⎕IO)+⍳¨(P+⎕IO)-I ⍝ Next inds
[15] ⍝ I←∊(I-⎕IO)+⍳¨(P+⎕IO)-I ⍝ If APL2
[16] →L1
∇
This function is more than twice as fast as any of the previous functions (on 4 COMBIN 5), but is still 5 times slower than the fastest function submitted.
The function from Ken Gradney loops from the last column to the first. Given APL’s right-to-left orientation, it’s perhaps fitting that this function should be the fastest. The moral: screw your head on backwards before working with APL.
For 3 COMBIN 5, here’s how Ken’s function proceeds:
In each iteration, you replicate the previous result by row-indexing, and then stick a new column to the left of it. For example, the second array above comes from the first via:
Again, after a healthy dose of induction, the general pattern becomes apparent (to Ken anyhow). Here’s the fastest function:
∇ R←M COMBIN N;D;E;F;G;P
[1] ⍝ Formula, working from last column:
[2] ⍝ R←(F/D),R[((1↑F)-F) PLUSIOTA F;]
[3] ⍝ where D starts as (M-1)+⍳1+N-M and
[4] ⍝ is decremented by 1 after each col;
[5] ⍝ F starts as (1+N-M)⍴1 and is cumu-
[6] ⍝ latively summed from the end after
[7] ⍝ each column. The algorithm for
[8] ⍝ E←A PLUSIOTA F is:
[9] ⍝ E←F/A-+\¯1↓0,F
[10] ⍝ E←E+⍳⍴E
[11] ⍝ giving (A[1]+⍳F[1]),(A[2]+⍳F[2]),...
[12] ⍝ The +\¯1↓0,F term for this col hap-
[13] ⍝ ens to be the A term for the next.
[14] E←(⍳P←N-R←M-1)-⎕IO
[15] D←R+⍳P
[16] R←(P,1)⍴D
[17] P←P⍴1
[18] L1:→(⎕IO>1↑D←D-1)⍴0
[19] P←+\P
[20] G←+\¯1↓0,F←⌽P
[21] E←F/E-G
[22] R←(F/D),R[E+⍳⍴E;]
[23] E←G
[24] →L1
∇
When asked what insight led to this algorithm, Ken said, "PLUSIOTA. The MONIOTA and PLUSIOTA function are the two most versatile and efficient utility functions in my bag of tricks."
For those of you who haven’t used these functions, here’s a listing of each:
∇ R←MONIOTA LEN
[1] ⍝ R←(⍳LEN[1]),(⍳LEN[2]),(⍳LEN[3]),...
[2] ⍝ In APL2: R←∊⍳¨LEN
[3] R←LEN/¯1↓0,+\LEN
[4] R←(⍳⍴R)-R
∇
∇ R←A PLUSIOTA LEN
[1] ⍝ R←(A[1]+⍳LEN[1]),(A[2]+⍳LEN[2]),...
[2] ⍝ In APL2: R←∊A+⍳¨LEN
[3] R←LEN/A-¯1↓0,+\LEN
[4] R←R+⍳⍴R
∇
MONIOTA 4 3 5
1 2 3 4 1 2 3 1 2 3 4 5
10 20 30 PLUSIOTA 4 3 5
11 12 13 14 21 22 23 31 32 33 34 35
Both functions work properly in either origin.