Get permutations from a vector


I have the same needs as the author of this topic, however in my case it would be better to write a custom permutations function within the Stan program.
My function is based on Python’s itertools.permutations which can be found here.

I wrote the following:

functions {
    array[] int permutations(int[] v, int r) {
        int n = size(v);

        if (r > n) {
            reject("r cannot be greater than the size of the vector");

        int j;
        int perm = 0;
        # number of possible permutations = n! / (n-r)!
        int perms = to_int(tgamma(n + 1)  / tgamma(n - r + 1)); 
        array[perms, r] int out;
        int indices[n] = linspaced_int_array(n, 1, n);
        int cycles[r] = reverse(linspaced_int_array(r, n - r, n));

        out[0,] = v;

        while (n) {
            perm += 1;

            for (i in reverse(linspaced_int_array(r - 1, 1, r))) {
                cycles[i] -= 1;

                if (cycles[i] == 0) {
                    indices[i:] = append_array(indices[i+1:], indices[i:i+1]);    
                    cycles[i] = n - i;
                } else {
                    j = cycles[i];
                    indices[i] = indices[-j];
                    indices[-j] = indices[i];

                    for (p in indices[:r]) {
                        out[perm, p] = v[p];

            return out;

I am not sure how to test and compare it with itertools.combinations. Do you see any errors or room for improvement?

1 Like

This is a notoriously tricky function to work. And the return will be too big if the input size is too big.

Step one is to get the function to work. I wasn’t sure what the r was doing or all the code involving cycles, but that’s not so important.

The simplest way to test a function is in transformed data with a print statement. For example, for this one, you could do something like

transformed data {
  print("case A:\n", permutations([...], 2);

You could turn that into a unit test by checking the results are what you expect.

The return type should be array[ , ] int rather than array[] int given that’s the type you’re returning with out.

If the intention is to swap values at indices, this won’t work:

indices[i] = indices[-j];
indices[-j] = indices[i];

The second statement here is redundant. Both indices[i] and indices[-j] will be set to indices[-j] when this executes. You need to assign to a temp value if you want to swap. You also have to be careful in Stan that arrays stay the size that they were declared. I’m not sure this statement preserves that:

indices[i:] = append_array(indices[i+1:], indices[i:i+1]);
1 Like