Adjective Noun

My recent obsession has been around combinatorics. For those of you who may be unfamiliar, combinatorics is a branch of mathematics closely related to graph theory. If I had to explain it in a short sentence, I'd probably say it's about the different ways in which a set of elements can be enumerated or constructed. That's a gross generalisation, but it will do for now.

There are a whole host of combinatoric algorithms, and Raku has 2 of them in the core language: permutations and combinations. There's good reason why just those 2... They are among the most common, and most useful, but that's not to say the other's aren't useful, and when I found myself needing one of those other algorithms, it led me on my aforementioned obsession.

The first one I want to talk about is "combinations with repetitions". This algorithm could be described as... At a given ice cream shop, how many different ways can I order 2 scoops. Order of choices doesn't matter, so 'Vanilla and Chocolate' is the same as 'Chocolate and Vanilla'

As a general rule, when order doesn't matter, you're talking combinations. When order matters, you're talking permutations

Now, there exists a way to do this in Raku on RosettaCode, but I want to state that I did come up with a solution by myself first based on a something I read in the Python documentation, and it also helped me later realise that - upon seeing it - the RosettaCode snippet was incorrect.

So back to Python for a minute... It has a combinations_with_replacement function in the itertools core module. Lets see what it looks like.

>>> from itertools import *
>>> list(combinations_with_replacement('ABCD', 2))
[('A', 'A'), ('A', 'B'), ('A', 'C'), ('A', 'D'), ('B', 'B'),
 ('B', 'C'), ('B', 'D'), ('C', 'C'), ('C', 'D'), ('D', 'D')]

In the itertools documentation for this function, it mentions that the result can be "expressed as a subsequence of product() after filtering entries where the elements are not in sorted order". In Raku, using the cross (Cartesian product) meta-operator ([X]), I came up with this nifty one-liner.

> sub cwr(@l, $k) { ([X] ^@l xx $k).unique(:as(~*.sort)).map({ @l[|$_] }) }
> cwr(<A B C D>, 2)
((A A) (A B) (A C) (A D) (B B) (B C) (B D) (C C) (C D) (D D))

I started by creating $k copies of my list indices, then create a Cartesian product of those lists, keeping unique ones (based on the stringified sorted order). I then use those indices to get the elements from the original list.

For the couple of benchmarks I ran (on admittedly small datasets), doing .unique(:as(~*.sort)) was slightly faster than doing something like .grep({ [≤] $_ }). In a pinch, this little snippet will do the trick, but it's also quite clear that I'm generating a bunch of data that I just throw away, so it can never be truly efficient.

Now take a look at the Raku snippet on RosettaCode for comparison. At the time of writing, it looked like this.

[X](@S xx $k).unique(as => *.sort.cache, with => &[eqv])

It certainly looks similar enough, and initially when I tried it out it seemed to work... However I quickly realised it had a flaw.

> sub ros(@S, $k) { [X](@S xx $k).unique(as => *.sort.cache, with => &[eqv]) }
> ros([0,1,2,3], 2)
((0 0) (0 1) (0 2) (0 3) (1 1) (1 2) (1 3) (2 2) (2 3) (3 3))
> ros([1,1,1,1], 2)
((1 1))
> cwr([1,1,1,1], 2)
((1 1) (1 1) (1 1) (1 1) (1 1) (1 1) (1 1) (1 1) (1 1) (1 1))

And here's Python just for good measure

>>> list(combinations_with_replacement([1,1,1,1], 2))
[(1, 1), (1, 1), (1, 1), (1, 1), (1, 1), (1, 1), (1, 1), (1, 1), (1, 1), (1, 1)]

Now I suppose you could argue that it's a combination, so order doesn't matter, but to push my ice cream analogy... Say your ice cream shop only has one flavour, but it has four buckets of that flavour. This algorithm is concerned with the different ways you can take two scoops in terms of which buckets you scoop from, so this RosettaCode snippet is slightly broken.

There's also a recursive version on RosettaCode, which I've included below.

proto combs_with_rep(UInt, @) {*}
multi combs_with_rep(0,  @) { () }
multi combs_with_rep(1,  @a) { map { $_, }, @a }
multi combs_with_rep($,  []) { () }
multi combs_with_rep($n, [$head, *@tail]) {
    |combs_with_rep($n - 1, ($head, |@tail)).map({ $head, |@_ }),
    |combs_with_rep($n, @tail);
}

say combs_with_rep(2, [1, 1, 1, 1]);

# OUTPUT: ((1 1) (1 1) (1 1) (1 1) (1 1) (1 1) (1 1) (1 1) (1 1) (1 1))

Apart from the minor difference of taking the list as the second argument, this function performs correctly, but it's slower than my one-liner (at least in the few benchmarks I ran).

I committed to finding a faster and more efficient algorithm. Most of the other snippets on RosettaCode were recursive functions. I knew that iterative code was generally more performant than recursive, so I kept looking for a iterative solution. I noticed the C++ version, and converted it to Raku. It was faster, but eventually I came upon another algorithm which - when converted to Raku - benched even faster.

I'm sure those of you of the more Computer Science persuasion could have told me where to look, but several sites referenced Donald Knuth's The Art of Computer Programming books. Specifically, "Fascicle 2: Generating All Tuples and Permutations" and "Fascicle 3: Generating All Combinations and Partitions". I had a look and it seems the books don't straight-up give you some code, but rather more-or-less describe an algorithm. I suspect most the algorithms in use for this sequence are interpretations of the algorithm described.

So far, the fastest algorithm I found (as far as pure Raku benchmarks are concerned) is the following

sub cwr(@list, int $k) {
    gather {
        my @idx = 0 xx $k;
        take @list[@idx];
        my int $e = @list.end;
        loop {
            if @idx[$k - 1] < $e {
                @idx[$k - 1]++;
            }
            else {
                loop (my int $j = $k - 2; $j ≥ 0; $j--) {
                    last if @idx[$j] != $e;
                }
                last if $j < 0;
                @idx[$j]++;
                loop ($j += 1; $j < $k; $j++) {
                    @idx[$j] = @idx[$j - 1];
                }
            }
            take @list[@idx];
        }
    }
}

This algorithm does not take into account what should happen when $k ≤ 0 or @list is empty, but those can be added fairly trivially. Upon gazing at this code, your first thought might be "Egads man! Why are you using c-style loops", and the reason should be obvious. I benched it and it was faster than using a Range.

UPDATE - Jun-2020: This may no longer be the case, as the extraordinary lizmat has made several optimisations to Ranges in Raku since this article was published

So far, this is the fastest algorithm I benched in pure Raku, but can it go faster? It can if we move beyond pure Raku, and into the world of NQP. NQP is the sub-language that forms the building blocks of the Raku language. It's more difficult to write, but you'll find that most expensive operations in the Raku core are written in NQP (including the existing permutations and combinations built-ins).

Writing these algorithms in NQP was a challenge for me. I hadn't written NQP before, so I mainly copied what I'd seen in the Rakudo code base, and referred to the list of NQP Opcodes page when necessary. The reward for my efforts was functions that ran much faster. I converted the few different algorithms I found to NQP, but the the above one was also (marginally) the fastest in NQP.

This post is already quite long enough, so I don't want to dump a whole page of NQP code here, but while my mind still has a hankering for combinatorics, I figure I might tackle a few more algorithms and make a module out of it. I'm gonna keep it off the ecosystem until it's a bit more fleshed out, but if you're interested in combinatorics, and/or a deft hand with NQP, pull requests are welcome.

Lastly, I would remiss to mention that Perl has a Algorithm::Combinatorics module, which has just about any combinatoric algorithm you could need written in fast XS, and it can be used just fine in Raku via Inline::Perl5.

> use Algorithm::Combinatorics:from<Perl5> 'combinations_with_repetition'
> combinations_with_repetition(<A B C D>, 2)
[[A A] [A B] [A C] [A D] [B B] [B C] [B D] [C C] [C D] [D D]]

Once imported, it's combinations_with_repetition function is at least twice as fast as my NQP algorithm. Which is to say, if you have a C compiler installed, and have Perl built with the right flags to support Inline::Perl5, you can install that module and use it today.

For the rest of you who need/want a fast native combinatorics library, I hope to implement as many of those algorithms as I can in NQP to make a Raku equivalent of Algorithm::Combinatorics. NQP still won't top C for performance, but Raku will allow very nice functionality, such as lazy evaluation, and Iterator shortcuts like count-only (which I've already implemented).

use Math::Combinatorics 'multicombinations';
use Algorithm::Combinatorics:from<Perl5> 'combinations_with_repetition';

sub time-it($desc, &func) {
    say "$desc: {func()} (%s seconds)".sprintf: now - ENTER now;
}

time-it 'Raku', { multicombinations(^16, 10).elems }
time-it 'Perl', { combinations_with_repetition(^16, 10).elems }

#`[ OUTPUT:
Raku: 3268760 (0.0043160 seconds)
Perl: 3268760 (5.1210621 seconds)
]

For algorithms that can find the "nth" iteration, then the skip methods can also be implemented for fast indexing into the sequence.

I'm not sure about some of the names, though. For example, combinations-with-replacement is quite a mouthful. I've seen it referred to as multicombinations in some circles - so that's what I'm using - but I'm not entirely sure if it means the same thing. If you're familiar with combinatorics, let me know if that name makes sense.

I've purposely labeled this article "Part 1" to force gently remind myself to keep working on this stuff. I'll probably be tackling some permutation of the permutations algorithm next.

To be continued...