I’m just wondering if anyone has implemented a Stan function that returns all permutations of a vector?

Background is that I am trying to improve the maxi-challenge component of my drag race predict-a-looza model. I would like to make the model able to handle an arbitrary number of winners and losers and I think something like the following code would work, if only there was a `get_permutations`

function.

```
functions {
real rank_logit_lp(vector abilities){
real out = 0;
for (i in 1:(rows(abilities) - 1)){
out += abilities[i] - log_sum_exp(abilities[i:]);
}
return out;
}
real maxi_lp(vector abilities, int n_winner, int n_loser){
int n_safe = rows(abilities) - n_winner - n_loser;
vector[n_winner] winners = abilities[:n_winner];
vector[n_safe] safe = abilities[n_winner+1:n_loser-1];
vector[n_loser] losers = abilities[n_loser:];
vector[n_winner] winner_orders[tgamma(n_winner+1)] = get_permutations(winners);
vector[n_safe] safe_orders[tgamma(n_safe+1)] = get_permutations(safe);
vector[n_winner] loser_orders[tgamma(n_loser+1)] = get_permutations(losers);
real order_lps[size(winner_orders), size(safe_orders), size(loser_orders)];
for (iw in size(winner_orders)){
for (is in size(safe_orders)){
for (il in size(loser_orders)){
order_lps[iw, is, il] =
rank_logit_lp(append_row(append_row(winner_orders[iw], safe_orders[is]), loser_orders[il]));
}
}
}
return log_sum_exp(to_array_1d(order_lps));
}
}
```

The idea is to find the log probability of the winners ranking higher than the safe contestants, who in turn beat the losers, given all the contestants’ abilities. Based on this paper I think this can be done by assuming an underlying rank logit model and considering every full ordering compatible with the known rankings. The problem is that I can’t work out how to enumerate all the orderings.