I did not see it mentioned in the article, but I wonder why or if the team rejected a permutation function on the index. The shuffled dataset would then be Xs[i] = X[f(i)]
One reason cited in TFA for the half-shuffled files approach is that it's easy to rotate old data out of and new data into the half-shuffled files.
The downside is that it still needs one random read access per element, so cache friendly hierarchical algorithms, like the one described in the post, are probably still faster for on disk data.
For example, using the permuted the index mod M rather than random draw in the first pass could avoid the whole issue of oversized piles and the extra wasted space per pile.