Successive sampling without replacement: draw one item with probability proportional to its weight, remove it, renormalise over the survivors, repeat.
What this distribution is, and one thing it is not. Drawn to exhaustion it is the
weighted shuffle, also called the Plackett-Luce order. The FIRST draw is exactly
proportional to weight; later draws are proportional among the survivors. It is not
inclusion-probability-proportional-to-size: an item with twice the weight is not twice as
likely to appear somewhere in the first k, and no construction with this interface can
make it so. Most callers who reach for "weighted sampling without replacement" want this
one, but the distinction is worth knowing before relying on it.
type(pf_weighted_draw) :: d
integer :: item
logical :: ok
call d%init(weights, seed)
do
call d%next(item, ok)
if (.not. ok) exit ! the population is exhausted
! ... judge item; exit on acceptance ...
end do
The tree is seed-independent, and that is the whole reason this type exists. It is built
from the weights alone; the seed only steers the descent. So a second sequence over the same
weights costs O(k log n) to restore rather than an O(n) rebuild, which is what makes an
outer loop of many short sequences -- the shape this type was designed for -- nearly free.
Use %reseed for that, holding the seed fixed and passing the outer iteration as stream.
Cost: %init is O(n), %next is O(log n), %reset and %reseed are O(k log n)
in the draws already taken.
Serial by nature, and there is no threads=. Draw k+1 cannot be produced until draw
k has been removed, so there is no parallelism inside one sequence to expose; an argument
that was accepted and ignored would be worse than its absence. Independent sequences ARE
parallel, and that is where the threads go -- see the note below.
To run many sequences in parallel, give each thread its own sampler, built inside the
parallel region. %next mutates the tree, so one sampler cannot serve two threads.
type(pf_weighted_draw), allocatable :: dd(:)
integer :: nt, tid, j
nt = 1
!$ nt = omp_get_max_threads()
allocate(dd(nt))
!$omp parallel do default(shared) private(tid)
do tid = 1, nt
call dd(tid)%init(weights, seed) ! nt builds, wall-clock of ONE
end do
!$omp parallel do default(shared) private(tid, item, ok) schedule(dynamic)
do j = 1, n_outer
tid = 1
!$ tid = omp_get_thread_num() + 1
call dd(tid)%reseed(seed, stream=int(j, int64))
do
call dd(tid)%next(item, ok)
if (.not. ok) exit
end do
end do
Never put this type in an OpenMP private() clause. A private copy is a FRESH object,
not a copy of yours: measured on gfortran 15.2.1 and ifx 2026.1.1, the tree comes back
allocated to the right shape with UNINITIALISED contents and every scalar reset to its
default. Nothing aborts and every drawn item still looks valid, so the failure is a wrong
answer with no symptom. The shared per-thread array above avoids the privatisation machinery
entirely; firstprivate() also copies correctly on both compilers, if you prefer it.
Because each sequence is named by its coordinates, the result does not depend on the
schedule or the thread count. schedule(dynamic) above costs nothing in reproducibility.
"Zero-weight" means a weight below wd_min_weight (2.05e-307), not exactly zero. Below
that bound -log(u)/w would overflow, so such a weight cannot be ordered at all; and against
weights of order one its draw chance is of order 1e-307, which no finite sample can tell
from zero. Thresholding rather than testing > 0 also makes the classification independent
of the floating-point model, since a build that flushes denormals reads them as 0.0.
Zero-weight items are kept out of the tree and handed out last, in uniform random order, so draining a population is always a genuine permutation of every item. Keeping them out is not an optimisation: a zero leaf would be indistinguishable from a spent one, and the descent's liveness test reads exactly that.
Prepares the sampler over weights. O(n). stream defaults to 0.
%init with no stream; the stream is 0.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(pf_weighted_draw), | intent(inout) | :: | self |
the sampler |
||
| real(kind=real64), | intent(in) | :: | weights(:) |
one weight per item |
||
| integer(kind=int64), | intent(in) | :: | seed |
the seed |
%init with an integer(int32) stream.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(pf_weighted_draw), | intent(inout) | :: | self |
the sampler |
||
| real(kind=real64), | intent(in) | :: | weights(:) |
one weight per item |
||
| integer(kind=int64), | intent(in) | :: | seed |
the seed |
||
| integer(kind=int32), | intent(in) | :: | stream |
which sequence; sign-extends |
%init with an integer(int64) stream. The form the other two delegate to.
Deliberately intent(inout) with an explicit called-twice guard, not intent(out).
intent(out) would reset the object on entry and make a second %init silently succeed,
quietly discarding a sequence in progress; %reseed is the supported way to reuse a
sampler and is O(k log n) where a rebuild is O(n), so a caller reaching for %init
twice is nearly always reaching for the wrong one.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(pf_weighted_draw), | intent(inout) | :: | self |
the sampler |
||
| real(kind=real64), | intent(in) | :: | weights(:) |
one weight per item |
||
| integer(kind=int64), | intent(in) | :: | seed |
the seed |
||
| integer(kind=int64), | intent(in) | :: | stream |
which sequence |
Draws the next item. O(log n). ok is .false. once every item has been drawn.
%next into an integer(int32) item index.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(pf_weighted_draw), | intent(inout) | :: | self |
the sampler |
||
| integer(kind=int32), | intent(out) | :: | item |
the item drawn, or 0 when exhausted |
||
| logical, | intent(out), | optional | :: | ok |
|
%next into an integer(int64) item index.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(pf_weighted_draw), | intent(inout) | :: | self |
the sampler |
||
| integer(kind=int64), | intent(out) | :: | item |
the item drawn, or 0 when exhausted |
||
| logical, | intent(out), | optional | :: | ok |
|
Restarts the SAME sequence. O(k log n).
Restarts the SAME sequence: the next %next returns what the first one did.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(pf_weighted_draw), | intent(inout) | :: | self |
the sampler |
Starts a DIFFERENT sequence over the same weights. O(k log n); stream defaults to 0.
%reseed with no stream; the stream is 0.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(pf_weighted_draw), | intent(inout) | :: | self |
the sampler |
||
| integer(kind=int64), | intent(in) | :: | seed |
the new seed |
%reseed with an integer(int32) stream.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(pf_weighted_draw), | intent(inout) | :: | self |
the sampler |
||
| integer(kind=int64), | intent(in) | :: | seed |
the new seed |
||
| integer(kind=int32), | intent(in) | :: | stream |
the new stream; sign-extends |
%reseed with an integer(int64) stream. The form the other two delegate to.
Works on a pristine sampler, one that has never drawn: the journal is empty, the replay is a no-op and only the coordinates change. That is not an edge case to tolerate but the normal path for the per-thread-array idiom, where every sampler is freshly built and then immediately reseeded on its first outer iteration.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(pf_weighted_draw), | intent(inout) | :: | self |
the sampler |
||
| integer(kind=int64), | intent(in) | :: | seed |
the new seed |
||
| integer(kind=int64), | intent(in) | :: | stream |
the new stream |
Items not yet drawn; exact, O(1).
Items not yet drawn, zero-weight ones included. Exact and O(1) -- never a weight sum.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(pf_weighted_draw), | intent(in) | :: | self |
the sampler |
items still available
type, public :: pf_weighted_draw private !> The segment tree: `2*p2` nodes, leaf `j` at `st(p2 + j - 1)`, every internal node the !! sum of its two children. Node `1` is unused padding's parent -- the root is `st(1)`. real(real64), allocatable :: st(:) !> Leaf `j` holds the weight of item `leaf_item(j)`. **Unallocated when no weight is !! zero**, which is the common case; leaf `j` is then item `j` and the map is skipped. integer(int64), allocatable :: leaf_item(:) !> The zero-weight items, in input order. Unallocated when there are none. integer(int64), allocatable :: zero_item(:) !> Undo journal: the leaf node zeroed by each draw so far, in order. !! !! **One entry per DRAW, not one per tree write.** Undoing a draw restores its leaf and !! then recomputes that leaf's ancestors from their children, exactly as the forward !! direction does -- so the `O(log n)` ancestor writes need not be journalled at all. That !! is 1 entry where the obvious design has `1 + log2(n)`, which at `n = 10**6` is a 21x !! difference in what a fully drained sampler holds, and it is why no capacity policy or !! rebuild fallback is needed. It is exact rather than approximate: every node is a pure !! function of its two children, so replaying the leaves reproduces the built tree bit for !! bit -- `test_reset_restores_tree_exactly` asserts precisely that. integer(int64), allocatable :: jr_pos(:) !> Undo journal: the weight each of those leaves held. real(real64), allocatable :: jr_val(:) integer(int64) :: p2 = 0 !! leaves in the tree; the least power of two `>= npos` integer(int64) :: npos = 0 !! items with a strictly positive weight integer(int64) :: nzero = 0 !! items with weight exactly zero integer(int64) :: live = 0 !! positive-weight items not yet drawn integer(int64) :: ndrawn = 0 !! draws taken in this sequence, zero-weight tail included integer(int64) :: jr_n = 0 !! journal entries in use integer(int64) :: wseed = 0 !! the seed steering the descent integer(int64) :: wstream = 0 !! which sequence of that seed integer(int64) :: zkey = 0 !! derived seed ordering the zero-weight tail logical :: ready = .false. !! `%init` has run; guards every other entry point contains procedure, private :: init_base => wd_init_base !! `%init` with no stream procedure, private :: init_s32 => wd_init_s32 !! `%init` with an `int32` stream procedure, private :: init_s64 => wd_init_s64 !! `%init` with an `int64` stream !> Prepares the sampler over `weights`. `O(n)`. `stream` defaults to 0. generic :: init => init_base, init_s32, init_s64 procedure, private :: next_i32 => wd_next_i32 !! `%next` into an `int32` item procedure, private :: next_i64 => wd_next_i64 !! `%next` into an `int64` item !> Draws the next item. `O(log n)`. `ok` is `.false.` once every item has been drawn. generic :: next => next_i32, next_i64 procedure :: reset => wd_reset !! Restarts the SAME sequence. `O(k log n)`. procedure, private :: reseed_base => wd_reseed_base !! `%reseed` with no stream procedure, private :: reseed_s32 => wd_reseed_s32 !! `%reseed` with an `int32` stream procedure, private :: reseed_s64 => wd_reseed_s64 !! `%reseed` with an `int64` stream !> Starts a DIFFERENT sequence over the same weights. `O(k log n)`; `stream` defaults to 0. generic :: reseed => reseed_base, reseed_s32, reseed_s64 procedure :: remaining => wd_remaining !! Items not yet drawn; exact, `O(1)`. end type pf_weighted_draw