I'm not sure what distribution is desired. Here is one that chooses uniformly among all distinction permutations of partitions of all integers m <= n
into nonnegative parts no greater than k
. This is done by transposing the Young tableaux for partitions of 2n
into at most k + 1
parts, and subtracting 1
to get the parts to be between 0
and k
. We then have to select those whose sums are at most n
.
ClearAll[getparts, weights];
mem : getparts[n_, k_] := mem = getpartsC[
PadRight[IntegerPartitions[2 n, k + 1], {Automatic, k + 1}],
n];
getpartsC = Compile[{{partitions, _Integer, 2}, {n, _Integer}},
Select[
Function[p, Total@Transpose[UnitStep[p - #] & /@ Range[n]] - 1] /@ partitions,
Last[#] >= 0 &]
];
mem : weights[n_, k_] := mem = Multinomial @@@ (Tally /@ getparts[n, k])[[All, All, 2]];
randsamp = Function[{n, howmany}, Ordering /@ RandomReal[1, {howmany, n}]];
ClearAll[nparts];
nparts[n_, k_, nsamp_] :=
With[{parts = getparts[n, k]},
With[{p = RandomChoice[weights[n, k] -> parts, nsamp]},
Compile[{{pp, _Integer, 2}, {samp, _Integer, 2}},
MapThread[#1[[#2]] &, {pp, samp}]][p, randsamp[n, nsamp]]
]];
A one-time cost: computing the basic partitions:
getparts[32, 20] // Dimensions // AbsoluteTiming
(* {32.0588, {43202, 32}} *)
Generating samples is fairly quick (it seems to be an order of magnitude faster than the Combinatorica`
code of Kellen Meyers, although there is a claim in a comment, without supporting code, that suggests it might be about the same). A million samples should take a total of about six minutes.
Table[nparts[32, 20, 32], {10^3}] // Dimensions // AbsoluteTiming
(* {0.35359, {1000, 32, 32}} *)
Getting them in chunks is faster, if you have the memory. A thousand 32x32 arrays is not much, but a million takes more memory than I have RAM, and therefore would be much slower to do all at once.
ArrayReshape[
nparts[32, 20, 32*10^3],
{10^3, 32, 32}] // Dimensions // AbsoluteTiming
(* {0.070147, {1000, 32, 32}} *)
Comparison of this approach and the Combinatorica`
approach (MakeArray
):
10^3 AbsoluteTiming[Do[MakeArray[], {10^3}]]
10^3 AbsoluteTiming[Do[nparts[32, 20, 32], {10^3}]]
10^3 AbsoluteTiming[ArrayReshape[nparts[32, 20, 32*10^3], {10^3, 32, 32}];]
(*
{3146.06, 1000 Null}
{357.585, 1000 Null}
{78.877, 1000 Null}
*)
Note that you need to add 30 sec. to the nparts
results for the one-time computation of the basic partitions.