I implemented Bertrand's solution. One other option would be to add a keyword argument to permuted_ols, in order to return the scores for all permutations. Then, we could compare h0 to the theoretical F distribution.
With the current version of the score, h0 is the distribution of "max(score)" and we (at least I) do not know its theoretical distribution, and so it is difficult to check.