pf_weighted_permutation Interface

public interface pf_weighted_permutation

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).

Module Procedures