ILE: make --seed actually reproducible on GPU

[oshaughnessy-junior]

Problem

--seed does not make GPU ILE runs reproducible.

bin/integrate_likelihood_extrinsic_batchmode implemented it as a bare
numpy.random.seed(opts.seed). The samplers, however, draw through the array
backend they were configured with — self.xpy on an instance, xpy_default at
module scope — and that backend is cupy whenever the job runs on a GPU:

• self.xpy.random.uniform in mcsamplerGPU.draw_simplified (the inverse-CDF
  proposal — this is the draw that decides the answer)
• xpy_default.random.uniform in mcsamplerAdaptiveVolume.sample_from_bins
• self.xpy.random.choice in the fair-draw / extrinsic-resample paths of
  mcsamplerGPU, mcsamplerAV, mcsamplerPortfolio, mcsamplerEnsemble

cupy keeps its own global generator, per device, which numpy.random.seed never
touches. Reported measurement: two byte-identical GPU invocations with
--seed 101 gave lnL = 75.857 (n_eff 1.02) and lnL = 71.687 (n_eff 2.04), a
4.17 nat spread. Reproduced here on the ILE-GPU-Paper demo at 0.287 nats.

Beyond being unbisectable, this silently invalidates any paired / replicate-seed
comparison design run on GPU, because the "same seed" arms are not in fact
paired. It was found while trying to run a stencil comparison campaign, which
had to be restricted to CPU as a result.

None of the samplers hold their own Generator objects — they all use
module-global backend state — so global seeding is the correct fix, not
per-instance plumbing.

Seeding the RNGs was necessary but not sufficient

After seeding cupy, same-seed GPU runs still disagreed at ~1e-14.

cupy.bincount(..., weights=...) accumulates through float atomicAdd, whose
ordering follows GPU thread scheduling. Measured on an RTX 2080 Ti, 7/7 repeated
calls on byte-identical input disagreed at ~2e-15 relative (the unweighted
bincount is exact — integer atomics). RIFT hits this in
RIFT/likelihood/vectorized_general_tools.histogram, which builds the adapted
sampling CDF — so the ULP noise is injected into every subsequent draw.

What isolated it: pinning PYTHONHASHSEED did not help, and n_eff diverged
while sqrt(2 lnLmax) stayed bit-identical — same sample stream, so the defect
was in an accumulator, not the RNG.

Changes

• new RIFT/integrators/seeding.py — seed_everything(seed) seeds every
  backend a RIFT sampler can reach (python, numpy, cupy, torch) and returns/prints
  what it actually managed to seed. The failure mode this fixes was invisible,
  so the report is deliberate. Absent backends are reported, not raised: a
  CPU-only install has no cupy, and only mcsamplerNFlow needs torch.
• RIFT/likelihood/vectorized_general_tools.py — deterministic weighted-bincount
  branch (argsort + prefix-sum differencing), behind a module flag defaulting to
  False. seed_everything flips it, so unseeded production runs pay nothing.
  Measured cost 1.2–1.5x, on a call that fires once per parameter per adaptation —
  far off the likelihood hot path. The flag is pushed from seeding rather than
  pulled, so RIFT.likelihood keeps no dependency on RIFT.integrators.
• the four drivers that implement --seed: integrate_likelihood_extrinsic_batchmode,
  ..._batchmode_lisa, integrate_likelihood_extrinsic, ile_postproc_add_time.
• new test/integrators/test_seeding_reproducibility.py — GPU tests skip
  cleanly on a CPU-only machine.

Verification

ILE-GPU-Paper demo (.travis/ILE-GPU-Paper/demos), on ldas-pcdev13,
RTX 2080 Ti / sm_75, cupy 10.6 / CUDA 11.2, CUDA_VISIBLE_DEVICES=3.

	seed 101 run A	seed 101 run B	seed 202
before	lnL 73.806544 (n_eff 5.93)	73.519973 (n_eff 10.78)	74.085756
after	74.867448087045 (n_eff 2.0391475242)	bit-identical	74.991472

• GPU, --sampler-method adaptive_cartesian_gpu: seeds 101 and 202 each reproduce
  bit-identical output files across repeat runs; the two seeds differ from
  each other, so this seeds rather than freezes.
• GPU, --sampler-method AV: same — bit-identical at seed 101, differs at 202.
• GPU, fully adaptive (no --no-adapt-after-first, --n-max 4e5, 40
  iterations): every iteration diagnostic line bit-identical across the pair. The
  fix survives the adaptation feedback loop, not just a frozen proposal.
• CPU: still bit-identical at a fixed seed, still seed-sensitive. lnL moved by
  9e-13 nats — the rounding difference of the new summation order.

New suite: 9 passed on GPU; 5 passed / 4 skipped CPU-only.

Neighbouring integrator suites run base-vs-branch on both CPU and GPU are
identical (51 passed CPU; 49 passed / 2 failed GPU, both failures pre-existing on
rift_O4d). test_mcsampler_gpu.py, test_mcsampler_foridiots.py and
test_mcsamplerEnsemble*.py fail at collection on the base branch too —
pre-existing and unrelated.

Self-review findings (second commit)

An adversarial pass over the first commit turned up four things, all addressed:

1. The default path was no longer byte-for-byte the old code. The
   ascontiguousarray that materializes the unweighted broadcast_to view sat in
   histogram(), so an unseeded run allocated n_samples of weights it never
   previously needed. Moved inside the deterministic branch — with the flag off,
   _bincount_weighted is now literally the original xpy.bincount call.

2. Prefix-sum differencing has an accuracy cost the first commit did not state.
   A bin total is the difference of two partial sums both of order the grand total,
   so a bin's relative error is amplified by (total / bin). Measured against an
   exact rational reference at n_bins=100:

   weights	bin-total spread	deterministic	atomic
   exponential	~1	2e-14	7e-15
   exp(lnL)-peaked	~2e5	5e-11	3e-14

   Acceptable: this histogram is a proposal density, not an estimator — the
   importance weights correct for whatever it is, it is consumed as a 100-bin
   interpolated CDF, and --adapt-floor-level mixes a uniform component in on top.
   The atomic branch's extra accuracy is unusable anyway, being irreproducible.
   Now in the docstring and pinned by a test at 1e-9 so it cannot drift.

3. ile_postproc_add_time is dead code, not a driver this fixes. It reads
   opts.seed with no matching add_option, so it dies with AttributeError at
   that line. Adding --seed just moves the crash to the next undefined option
   (manual_logarithm_offset), which would be a misleading half-fix — reverted, and
   said so in a comment. The seed_everything swap stays so reviving it does not
   reintroduce the bug. Resurrecting that script is out of scope here.

4. Added the missing GPU coverage for the unweighted broadcast_to branch, raised
   the atomics probe to 32 repeats, dropped an unused exception binding.

Checked and found fine: no driver switches CUDA device after the seed call, so
seeding the current device is sufficient; out-of-range and empty inputs behave
identically in both branches.

Every output file after these changes is bit-identical to the first commit's — the
review moved no numbers. Final counts: 11 passed GPU, 6 passed / 5 skipped
CPU-only, 51 passed across neighbouring integrator suites (unchanged from base).

Known residual, deliberately not fixed here

statutils.bootstrap_lnZ_quantiles defaults rng_seed=None →
numpy.random.default_rng(None), which pulls OS entropy, so the lnZ_ci90
diagnostic is irreproducible even on CPU. It is reporting-only (it does not feed
lnL/lnZ) and is gated behind a large-sigma condition. Left alone by decision, to
avoid perturbing the existing numpy stream in statutils.

🤖 Generated with Claude Code