parquet_sampling Module

Element k of a uniform-looking permutation of 1 .. m, addressed by its coordinates.

The permutation counterpart of pf_random_at: a pure function of (seed, m, k), computed in constant time and constant memory, touching no array. Element k depends on no other element, so a loop over k may be run in any order, on any number of threads, and gives the same answer -- which is the property the whole module exists for, extended from draws to permutations.

!$omp parallel do
do k = 1, m
    perm(k) = pf_random_perm_at(seed, m, k)      ! same result at any thread count
end do

The first n values are a uniform random n-subset of 1 .. m, so a subset needs no separate machinery and no memory: ask for k = 1 .. n. Prefix consistency follows for free -- a size-3 subset is a prefix of a size-6 one from the same (seed, m) -- and a subset at n == m is the permutation, rather than merely agreeing with it.

m and k share their kind (integer(int32) or integer(int64)) and the result follows them; seed is always integer(int64). k outside [1, m] is clamped rather than reported, the same convention draw_or_1 uses and for the same reason: this is pure elemental and has no way to abort.

Uniformity comes in two grades, split at m = 20, and the split is visible in pf_random_perm_algorithm.

For m <= 20 the result is exactly uniform over all m! permutations -- one exactly uniform rank in [0, m!) composed with a bijection onto S_m, so this is a property of the construction rather than a measurement. 20 is where it stops because 20! is the last factorial an integer(int64) holds.

For m >= 21 a 64-bit seed cannot address m! permutations at all, so exact uniformity is not available to any construction with this signature. What is measured instead is that the result is indistinguishable from a uniform permutation under an all-cells chi-square, an order-4 tuple statistic, a parity test over the alternating group, and fixed-point, cycle-structure, position-uniformity, subset-membership and structural tests.

Consecutive m are independent, in both regimes. Asking for a permutation of 5 and one of 6 under the same seed gives two unrelated answers rather than two views of one draw.

Fills perm with the whole permutation of 1 .. size(perm) under seed.

The bulk form of pf_random_perm_at, and an identity rather than a second algorithm: perm(k) equals pf_random_perm_at(seed, size(perm), k) for every k, so the two may be mixed freely and a prefix of one is a prefix of the other. That is the same relationship pf_random_at and pf_random_fill_draws already have, and it is what keeps the two forms from drifting apart under a later optimisation.

perm is a rank-1 integer(int32) or integer(int64) array, intent(out); seed is integer(int64). A zero-sized perm is a defined no-op.

threads is optional and changes only how fast the array is filled, never what is in it -- and that is provable here rather than merely intended, which is unusual for a threading argument. perm(k) is a pure function of (seed, size(perm), k) and reads nothing another element writes, so the 1-thread and N-thread results are bit-identical and a program may switch thread counts while debugging without its permutation changing under it. Absent means automatic: as many threads as OpenMP offers, capped by parquet_set_random_threads, and 1 inside an OpenMP parallel region, since a nested region is the caller's business. An explicit threads= is honoured there too.

A work floor applies to both, and it is not a tuning detail. Below parquet_set_random_parallel_min_elements (default 1000) elements per thread the call runs serially, because threading a small permutation is measurably slower than not threading it -- so threads=8 on a 100-element array is honoured by being declined.

It is about 2.4x cheaper per element than calling the elemental form size(perm) times (machine B, gfortran 14.2.1, --profile release: 9.9 ns against 24.1 ns at m = 10**8), for two structural reasons rather than any tuning. Enumerating the domain in order hands the loop the (l, r) split that the scalar form has to recover with an integer division; and, the larger of the two, the width rule and the key schedule are derived once per call instead of once per element, which a pure elemental function has no way to avoid. The values are unchanged by either -- only the route to them is.

Fills idx with the first size(idx) elements of the permutation of 1 .. m under seed.

A uniform random subset of size size(idx) drawn without replacement, in uniform random order -- and it costs O(size(idx)) rather than O(m), which is the reason this permutation is coordinate-addressed at all: 1000 rows out of 10**12 reads 1000 elements and never materialises the population. pf_random_subset(idx, size(idx), seed) is pf_random_permutation(idx, seed), and a subset of size n is a prefix of one of size 2n.

idx is a rank-1 integer(int32) or integer(int64) array, intent(out); m is integer(int32) or integer(int64); seed is integer(int64). A zero-sized idx is a defined no-op and is not validated -- it asks for nothing, so no precondition applies to it. threads behaves exactly as it does on pf_random_permutation, including the work floor; see there. Note the floor is measured in elements PRODUCED, size(idx), not in m -- a 100-element subset of a population of 10**12 is 100 elements of work.

Three preconditions, all of which abort rather than truncate or wrap when idx is non-empty: m >= 1; size(idx) <= m, since a subset drawn without replacement cannot be larger than its population; and -- for an integer(int32) idx only -- m <= huge(int32), since an element may be any value in [1, m] and one above huge(int32) has nowhere to go. The last of those is why the integer(int32)-array/integer(int64)-m pairing is accepted at compile time and rejected at run time: refusing it in the interface would make a perfectly ordinary integer :: m fail to compile against an integer(int64) array.

Fills idx with size(idx) values drawn from 1 .. m with replacement.

The third member of the resampling trio, and the one whose construction is not a construction at all: drawing with replacement means size(idx) independent uniform integers in [1, m], with no dedup structure, no permutation and no sort. So this is the draw-axis integer bulk fill under another name, and the identity is exact and is asserted by the suite:

call pf_random_resample(idx, m, seed, stream)
call pf_random_fill_draws(seed, stream, idx, 1_int64, m)   ! the SAME values

What the name buys is the 1.4x-1.6x a caller loses by writing the obvious loop. Without it the natural code is do k = 1, n; idx(k) = pf_random_int_at(seed, i, 1, m, k); end do, which re-enciphers a Philox block for every value where the bulk form serves two draws from each one -- measured at 25.7 against 15.7 ns per value on machine B (gfortran 15.2.1, --profile release, 10M values) and 16.1 against 10.1 on machine A. The values are identical either way; only the route to them differs.

idx is a rank-1 integer(int32) or integer(int64) array, intent(out); m is integer(int32) or integer(int64); seed is integer(int64). stream is optional, integer(int32) or integer(int64), and defaults to 1 -- it selects which replicate this is, so replicate b is reproducible from (seed, b) alone, independent of how many replicates were asked for or in what order they ran. A zero-sized idx is a defined no-op and is not validated: it asks for nothing, so no precondition applies to it.

Two preconditions, both aborting rather than truncating or wrapping, and only when idx is non-empty: m >= 1; and -- for an integer(int32) idx only -- m <= huge(int32), since an element may be any value in [1, m] and one above huge(int32) has nowhere to go. The second is why the integer(int32)-array/integer(int64)-m pairing is accepted at compile time and rejected at run time, exactly as on pf_random_subset.

There is deliberately NO size(idx) <= m precondition, and its absence is the clearest statement of how this differs from its sibling. That bound belongs to a subset drawn without replacement; here size(idx) == m is the most ordinary bootstrap there is, and size(idx) far beyond m is perfectly meaningful. A guard copied across from pf_random_subset would refuse the procedure's main use.

stream is this procedure's replicate axis, and the siblings do not have one. pf_random_permutation and pf_random_subset are keyed by (seed, m) alone, so independent replicates of those come from pf_random_key(seed, b) instead. A resample is built on the draw axis, which carries a stream coordinate already, so it costs a caller one integer rather than a key derivation. Both routes are available here: stream = b and seed = pf_random_key(seed, b) are equally independent.

threads behaves exactly as it does on pf_random_permutation and pf_random_subset -- same parquet_set_random_threads default, same parquet_set_random_parallel_min_elements work floor, and the floor applies to an explicit request too, so a small resample stays serial however many workers are asked for. The result is bit-identical at every thread count, and by construction rather than by care: element k is a pure function of (seed, stream, k), so a thread filling elements lo .. hi is the serial fill started at draw lo. That is asserted rather than argued.

threads requires an explicit stream, and that is a language constraint rather than a choice. threads and stream are both integers in argument position 4, so a generic offering (idx, m, seed [, threads]) beside (idx, m, seed, stream) is rejected by the compiler -- "Ambiguous interfaces in generic interface 'pf_random_resample'" -- and giving the two dummies different keyword names does not rescue it, because keyword names do not make specifics distinguishable. Both spellings were compiled to confirm it. So write call pf_random_resample(idx, m, seed, 1, threads=8) for the default replicate; omitting stream produces a "no specific subroutine matches" error that does not explain itself. This is the same constraint the parquet_open_reader split answers one level along, where an optional argument differing only by kind could not disambiguate either.

Fills perm with a weighted random permutation of 1 .. size(weights).

The order family. perm(k) is the item drawn k-th by successive sampling: the first is proportional to weight, each later one proportional among the survivors. Drawn to exhaustion, that is the weighted shuffle. Same distribution as pf_weighted_draw -- see the note on realizations below.

perm is a rank-1 integer(int32) or integer(int64) array, intent(out), whose size must equal size(weights); weights is real(real64); seed is integer(int64); stream is optional and integer(int64) -- note the asymmetry with pf_weighted_draw, whose %init/%reseed take either kind. It is forced rather than chosen: an optional threads and an integer(int32) stream are not distinguishable as a positional fourth argument, so one generic cannot carry both, and threads= is worth more here than saving a cast. Write stream=int(j, int64) in a loop. threads behaves as it does elsewhere in this module, including the work floor.

It is a DIFFERENT REALIZATION from the sequential family, not a different distribution. pf_weighted_permutation(perm, w, seed) and a pf_weighted_draw drained under the same seed both draw from successive sampling, and for one seed they give different draws from it -- in the same way two different seeds would. They are INDEPENDENT at matched coordinates, not merely different, which is a stronger claim and one the two families had to be domain-separated to earn: each derives its own seed through pf_random_key, so neither shares a uniform with the other or with a caller drawing at coordinates it chose itself. Before that separation both read draw 1 of (seed, stream) and chose the lowest-weight item together on 66 % of the seeds where either did, against 0.078 % by chance -- with every marginal clean. test_families_independent is what holds this now; see wd_family_label. There is no prefix identity across the two families. Within the sequential family the prefix identity does hold; see pf_weighted_subset.

Construction: the exponential race. Give item i the key -log(u_i)/w_i for independent uniforms, and sort ascending. That is exactly successive sampling, and unlike the sequential form it is embarrassingly parallel and coordinate-addressed, so the answer is bit-identical at every thread count. Use this when you want many whole shuffles; use pf_weighted_draw when you want a few draws, or many short sequences over the same weights, where it is thousands of times cheaper.

To produce many shuffles, parallelise ACROSS them rather than within one. The key loop is memory-bound and scales sub-linearly, whereas independent shuffles are perfectly parallel. An OpenMP loop over sequences, each calling this with threads=1, beats a serial loop over threaded calls; the module's own auto-threading already resolves to serial inside an active parallel region, so this needs no special handling.

Zero-weight items come last, in uniform random order -- they can never be drawn while a positive weight remains, so a full shuffle is still a genuine permutation of every item.

Preconditions, all of which abort: size(perm) == size(weights); every weight finite and >= 0; at least one weight > 0; no weight so small that -log(u)/w overflows (below about 2e-307, which no real weighting reaches and which would otherwise tie several items at infinity); and, for an integer(int32) perm, size(weights) <= huge(1_int32).

Fills idx with the first size(idx) items of a weighted sequential draw over weights.

Defined as size(idx) calls to pf_weighted_draw%next, not as a second algorithm, so a subset of size k is a prefix of one of size 2k, and stopping a %next loop after k draws gives exactly this. That identity is asserted by the suite rather than intended.

idx is a rank-1 integer(int32) or integer(int64) array, intent(out); weights is real(real64), one per item; seed is integer(int64); stream is optional and is integer(int32) or integer(int64). The population is size(weights) and the number drawn is size(idx) -- note this differs from pf_random_subset, whose second argument is the POPULATION, because here the population is carried by the weights themselves.

Preconditions, all of which abort when idx is non-empty: size(idx) <= size(weights); every weight >= 0; at least one weight > 0; and, for an integer(int32) idx, size(weights) <= huge(1_int32). A zero-sized idx asks for nothing and is a defined no-op.

O(n + k log n), so it is the cheaper form whenever k is small against n; past roughly k = n/2 a caller wanting most of the population is better served by asking for the whole weighted permutation.

Successive sampling without replacement: draw one item with probability proportional to its weight, remove it, renormalise over the survivors, repeat.

Type-Bound Procedures

Exposes random_threads for testing. Test-only; no library code calls it.

Arguments

Return Value integer

Forces the Feistel round count. Test-only. n <= 0 restores the compiled-in value.

Arguments

Enables or disables the parity correction. Test-only; see parquet_debug_set_perm_rounds.

Arguments

Lowers the population ceiling the three int32-index refusals test against. Test-only.

Arguments

Sends m <= perm_exact_max through the Feistel. Test-only; see the two hooks above.

Arguments

Reports the permutation kernel's effective configuration. Test-only.

Arguments