Skip to content

Catch rift_O4d up with oshaughnessy-junior/rift_O4d - #161

Merged
oshaughn merged 179 commits into
oshaughn:rift_O4dfrom
oshaughnessy-junior:rift_O4d
Aug 7, 2026
Merged

Catch rift_O4d up with oshaughnessy-junior/rift_O4d#161
oshaughn merged 179 commits into
oshaughn:rift_O4dfrom
oshaughnessy-junior:rift_O4d

Conversation

@oshaughnessy-junior

Copy link
Copy Markdown

Summary

Bring oshaughn/rift_O4d up to date with the accumulated work on oshaughnessy-junior/rift_O4d.

This includes the recent integrator and JAX work, calibration and sample-rate fixes, test coverage, container survey/warmup tooling, and documentation updates.

O4c integration check

The payloads of the recent oshaughn/rift_O4c merges are already integrated into the junior rift_O4d branch (patch-equivalent commits):

  • Virgo calibration correction convention
  • honoring --srate-resample-time-marginalization
  • exact requested export-rate recovery

Comparison

  • Base: oshaughn/rift_O4d (85b451dc)
  • Head: oshaughnessy-junior/rift_O4d (420f0791)
  • 91 files changed, 12,260 insertions, 222 deletions

No additional code changes were made solely for this catch-up PR.

oshaughnessy-junior and others added 30 commits July 10, 2026 11:34
The cache was sieved by ifos=channel[:1] -- the site LETTER -- so the three ET detectors, which
share the 'E' site, all matched the same frames and reading E2:/E3:FAKE-STRAIN was served from the
E1 frame -> "channel not found / Wrong name".  After the site-letter sieve, narrow to the exact IFO
when the cache carries full-IFO observatories (E1/E2/E3/C1 from <IFO>-*.gwf frames).  Standard
single-letter observatories ('H','L','V') never equal the 2-char IFO, so ordinary H1/L1/V1 caches
are untouched.  Surfaced by the finite-size 3G sky-loc pool run: all CE+ET (C1+E1+E2+E3) jobs held
on output transfer because the ILE died reading E2.  Verified: CEpET finite reads all 4 detectors
and produces .dat + .xml.gz.
…-size (Path D)

Extend the AD-compatible jax_ile likelihood to the slow-rotation and
frequency-dependent finite-size detector-response models, so the full high-SNR
extrinsic analysis runs as one differentiable JAX graph (gradient-based NUTS).

Design: both models generalize the baseline the same way -- one small summed
"band" index (sidereal harmonic a=(p,n) / basis weight p) contracted with a
per-sample coefficient vector C[a]; the baseline is the A=1, C=[F] case. So a
single generalized kernel (_accumulate_unit_banded) handles both, dispatched by
a one-line feature check at the top of _accumulate_unit -- baseline path
byte-for-byte unchanged, every marginalization variant inherits the features.

- response_slowrot.py / response_freqresponse.py: pure-JAX ports of
  rotation_coefficients_vector / response_coefficients (differentiable in
  RA/DEC/psi; GMST(tref) a host constant). Validated vs numpy to ~1e-16.
- banded.py + wrapper.py: build banded JAXLikelihoodData from the cupy packed
  banks (pack_rotation_arrays / pack_freqresponse_arrays, reused verbatim) and
  one-call build_*_data_from_precompute builders.
- Validation ladder (test/jax/test_jax_slowrot*.py): JAX-nearest reproduces the
  cupy NoLoopWithRotation / FreqResponseNoLoop references to max|rel| ~1e-14;
  distmarg gradient vs finite-diff ~1e-7 (rot) / ~1e-5 (freqresp); jit/grad/
  hessian finite (Fisher forecast falls out). Baseline unchanged (2.27e-13).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
test/jax/demo_slowrot_highsnr_nuts.py: injects a finite-size BNS into CE+ET
(reusing slowrot_fs_lib), builds the banded JAX finite-size likelihood, and
drives multistart_nuts at SNR 100/300/1000 -- reporting n_eff, 90% sky area and
ring-aware sky diagnostics. Demonstrates NUTS achieving high n_eff (1835 at
SNR=100) where the production AdaptiveVolume MC gives n_eff~1. Runs on GPU in
the JAX container with this branch on PYTHONPATH.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…_nuts

Production fix for high-SNR NUTS on the extrinsic posterior. The distance-
marginalized angular posterior is a thin, curved sky ring entangled with the
psi/incl/phiref degeneracies; numpyro's DEFAULT DIAGONAL mass matrix leaves the
Hamiltonian geometry wildly anisotropic, so NUTS hits max_tree_depth (~1024
leapfrog steps) on essentially every sample and STALLS (a single SNR=300 chain
spun >2.5 h with no progress).

multistart_nuts now defaults to dense_mass=True (adapts the full 5x5 inverse-
mass matrix ~ posterior covariance during warmup, whitening the geometry) and
exposes max_tree_depth (int or (warmup,sampling) tuple) to bound the pre-
adaptation warmup window. Validated: the previously-stalled SNR=300 finite-size
CE+ET case now completes all 6 chains in ~25 min (was: never), converged, MAP
recovered 1.6 deg from the injected truth -- a ~30x speedup.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…NUTS demo

Separate the two 'n_eff' notions the high-SNR demo conflated:
 * multistart_nuts now also returns theta_per_chain, so callers can compute a
   POSTERIOR effective sample size (within-chain ESS summed over chains) -- the
   honest 'did NUTS resolve the posterior' diagnostic. AV gives ~1 effective
   posterior draw; NUTS gives hundreds (SNR=100: post_ESS=536 of 3000).
 * the demo now reports post_ESS as the headline and labels the Gaussian-mixture
   importance EVIDENCE neff as a separate, high-SNR-limited number (a Gaussian
   mixture cannot wrap the thin curved sky ring -- a documented jax_ile
   limitation, not a sampling failure), plus SNR-scaled finer sky binning.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…high-SNR NUTS

multistart_nuts gains rotate_phase=True: sample phase_p=phiref+psi, phase_m=
phiref-psi (each on [0,4pi)) instead of (psi,phiref), mapping back
psi=(phase_p-phase_m)/2, phiref=(phase_p+phase_m)/2 -- the JAX mirror of
production RIFT bin/integrate_likelihood_extrinsic_batchmode --internal-rotate-
phase. The (2,+/-2) likelihood depends on 2psi+/-2phiref, so the curved
polarization/phase degeneracy becomes AXIS-ALIGNED in (phase_p,phase_m) and the
dense mass matrix is near-diagonal.  Constant-Jacobian rotation -> flat prior
preserved.  Combined with sky_coords='network' (baseline-frame sky, straightens
the time-delay ring) this is the full high-SNR extrinsic reparameterization.

Validated on the finite-size CE+ET SNR=300 case: SKY effective sample size rises
from ~34 (naive dense-mass NUTS) to 948 with the rotation -- a ~28x sampling win
-- MAP recovered 1.6 deg from truth.  Demo (demo_slowrot_highsnr_nuts.py) now
enables sky_coords='network'+rotate_phase and reports sky_ESS separately.  Also
adds demo_slowrot_reparam.py (phimarg+Fisher-whiten path, slower) and
demo_slowrot_flowmc.py (flowMC comparison).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Report per-dimension posterior ESS (sindec/skyx/psi/incl/phiref) so the sky
sampling quality is separable from the orientation sector. Shows the phase-
rotation + network-coords win is concentrated in the SKY (sky_ESS ~700-1200
across SNR 100/300/1000) while the residual high-SNR ESS loss is confined to the
intrinsically-sharp psi/incl/phiref orientation directions.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ecovery)

multistart_nuts gains polish_seeds=True: gradient-ascend (+ Fisher-inverse
Newton, reusing _map_polish_4) each raw pilot seed onto its local MAP before
running NUTS.  At extreme SNR the sky posterior is a ~1/SNR-thin ring that a
finite prior pilot cannot land ON -- like AV, the best pilot draw sits many nats
below the peak (SNR=1000: seed lnL ~24500 below 0.5<d|d>), so a chain started
there samples the WRONG arc and the recovered sky sits degrees off truth.  A few
hundred exact-gradient + Newton steps climb from the broad basin onto the true
needle -- the whole point of having analytic gradients (which AV lacks).

Validated (finite-size CE+ET, phase-rot + network + dense mass): SNR=1000 sky
recovery MAP 9.01deg -> 1.61deg from truth (matching SNR 100/300), logZ nan ->
4.91e5 (finite, at the peak), sky_ESS 751.  So the full SNR 100/300/1000 ladder
now recovers the sky where the production AdaptiveVolume MC gives n_eff~1 and
never finds the peak -- the high-SNR 3G demonstration.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…calization

The histogram sky_area_90 is floored by the ~5 deg^2 bin size, so it cannot show
a compact high-SNR peak shrinking (it reports ~25-30 deg^2 at every SNR even when
the peak is sub-deg^2).  Add a covariance-ellipse 90% area (pi*chi2_0.9*sqrt(det
Cov) in a local tangent plane) that resolves sub-deg^2 localization from a few
thousand samples, so the demo shows the 3G sky-localization power tightening with
SNR.  Reported alongside the histogram value (which bounds a curved ring arc).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…amples)

test/jax/make_3g_figdata.py: runs the finite-size CE+ET injection at an SNR
ladder with the full high-SNR reparameterization stack (network sky coords +
phase rotation + dense mass + gradient seed-polish), and writes the .npz data
for the paper's 3G subsection figures -- 90% sky area vs network SNR (KDE +
covariance estimators) and per-SNR sky samples, plus full extrinsic samples at a
representative SNR (distance drawn per angular sample from its analytic
conditional).  Consumed by analyses/slowrot_finite-size/3g/plot_3g.py.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…art_nuts

Multi-start NUTS places one chain per mode, but at high SNR the sub-dominant
time-delay-ring images and amplitude-degeneracy branches sit many nats below the
true peak (exp(-O(100)) suppressed). Pooling the chains with EQUAL weight
over-represents those negligible modes, so a naive credible region / corner plot
is dominated by junk (e.g. the injected inclination ends up in a low-density gap
while the sampler's true peak is exactly on it).

multistart_nuts now returns post_weight: each chain weighted by its mode's peak
likelihood exp(peak_lnL_k - max_k peak_lnL)/n_k (a robust mode-evidence proxy for
comparable-width modes). theta stays the raw pooled draws; callers use post_weight
for credible regions / areas / corner plots. make_3g_figdata resamples by it, so
the effective posterior collapses the negligible modes -> a clean, compact
recovery on the injection (dominant-mode draws are 100% within 2deg of truth,
inclination recovered) instead of the mode-spread mess.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ruth seed) + 3G figure generator

Building the 3G paper demo surfaced three distinct high-SNR fidelity issues in the
differentiable freqresponse/rotation ILE, all fixed here:

* core.py: add 4-point cubic Lagrange gather (`_gather_cubic`, registered "cubic")
  mirroring production `_cubic_Q_window_numpy`. Linear interpolation UNDERSHOOTS the
  razor-sharp high-f rholm peak (worse than nearest); cubic matches cupy NoLoop to
  numerical precision (44533.7 exactly vs 42801 linear / 43411 nearest).

* wrapper.py: build the tvals grid as arange(-Nw,Nw)*deltaT (spacing EXACTLY deltaT)
  in all three build_*_from_precompute, not linspace(-iwh,iwh,npts) (spacing
  deltaT*npts/(npts-1)). The linspace half-sample drift biases the sky ~1 deg at
  high SNR; the 1e-13 baseline validation missed it (same linspace both sides).

* samplers.py: multistart_nuts gains `extra_seeds` (caller seeds prepended before
  the gradient polish; chain loop now iterates len(seeds)). At very high SNR + fine
  deltaT the true peak is thinner than a pilot draw can resolve, so a blind pilot can
  settle on a secondary mode nats below the global peak; a known seed near the true
  basin (production: intrinsic grid + coarse extrinsic pass) guarantees one chain
  characterizes it.

* make_3g_figdata.py: 3G figure-data generator honours SLOWROT_FMAX (rholm deltaT),
  SLOWROT_SELFCONSISTENT (self-consistent b_p*W_p injection), SLOWROT_SEED_TRUTH, and
  the MAP-mode-robust sky-area estimator; interp="cubic" throughout.

Diagnostics: debug_fs_consistency, debug_jax_vs_cupy_inj, debug_orientation_degeneracy,
debug_timeshift (Qmax/incl/time-shift/oversample sweeps that isolated the systematics:
the residual high-SNR ~0.16 deg sky offset is the recovery's rholm cubic time
interpolation of per-detector fractional-sample delays, removed by oversampling deltaT).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Jax slowrot : slow rotation finite-size/slow rotation response
New test/integrators/benchmark_integrators.py: analytic targets with known truth
(correlated Gaussian in any D, Rosenbrock, Gaussian mixture), a uniform adapter
over the heterogeneous sampler call/return conventions (default/AC/GMM/AV/NF), and
the paper's quality metrics (integral bias, fractional error, n_eff, Kish n_ESS,
efficiency n_eff/N_eval, N_eval/wallclock to target neff, and 1-D marginal JS
divergence vs analytic truth). Backend (CPU numpy vs GPU cupy) is chosen by the
caller's CUDA_VISIBLE_DEVICES, as in production ILE. Includes a warm_start hook for
cold-vs-warm bootstrap comparisons and a CLI.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Adds warm-start support to mcsamplerAdaptiveVolume so an integrate_log() call can
begin from a concentrated live-volume instead of carving it cold from the full
prior every time.  New methods:
  - bootstrap_from_samples(samples[,loglkl])      reliable: seed covers support
  - bootstrap_from_gaussian(mean,cov)/from_fisher  fast UNIMODAL seed (Fisher oracle)
  - bootstrap_from_gaussian_mixture(means,covs,wt) multimodal / flow-oracle seed
  - save_state(path)/load_state(path)              lightweight (~8KB) state reuse

Design notes:
  * VARAHA's live volume only ever CONTRACTS, so a warm start is only unbiased if
    its seeded volume is a superset of the true support.  bootstrap_from_samples
    and _from_gaussian_mixture cover the support; a single Fisher is unimodal only
    (documented) -- on a 3-component mixture it otherwise biases the integral low.
  * The seeded fractional volume V is set GEOMETRICALLY (n_bins/prod(nbins)) so the
    final normalization stays unbiased regardless of how the seed was produced.
  * Occupied bins are dilated by axis-neighbor layers as a coverage safety margin.
  * Cold-start behavior is unchanged (self._warm is None): CI 4-sigma gate intact.

Benchmark (test/integrators/test_AV_bootstrap.py, GPU A100): on unimodal targets
warm-start reaches the target neff in ~1.3x fewer likelihood evals and ~1.3x higher
efficiency; on the multimodal mixture the mixture/samples oracle keeps the integral
unbiased (cold bias_ln -0.10 -> warm +/-0.01).  Gains are bounded because VARAHA
already contracts quickly; the larger value is amortizing carving across a
production ILE fleet and enabling correctness-preserving proposal reuse.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…wiring/GPU

Makes the portfolio 'oracle' mechanism actually work and adds a robust oracle.

New:
  - unreliable_oracle/fisher_gaussian.py: FisherGaussianOracle -- proposes from a
    box-truncated N(mean, cov=Gamma^-1); the cheap Fisher substitute for a trained
    flow the paper describes.  Optional adapt=True refits (mean,cov) from weighted
    history. Returns a proper proposal density.
  - unreliable_oracle/__init__.py: package exports.
  - test/integrators/test_portfolio_oracle.py: quantitative needle-in-haystack test.

Fixes to the existing stubs (they raised / mis-sampled before):
  - puffball: 'if ln_weights:' on an ndarray raised; also selected ROWS (params)
    instead of COLUMNS (samples) in the high-weight cut. Now guards None and
    indexes samples.
  - hill_climber: np.random.choice(range(len(sample_array))) used ndim, not the
    sample count; ln_weights mask was misaligned with the truncated history.

Fixes to mcsamplerPortfolio oracle wiring + backend consistency:
  - ORACLE BLOCK: numpy.append results were DISCARDED, so oracle proposals never
    actually entered the training data.  Now proposals + their true lnL are
    appended to rvs_train/log_weights_train (host), and members adapt on them.
    Oracles only propose -> they can never bias the integral.
  - portfolio_default_weights: guard sum(n_ess-1)==0 (all members degenerate on a
    hard target's first chunk) -> was producing nan weights -> negative per-member
    sample counts -> crash.
  - GPU path: the portfolio aggregates on the host but members return cupy; the
    old isinstance(type(x),type(y)) guard never fired, colliding numpy/cupy.  Now
    member draws + lnL are brought to host and the running estimate is forced onto
    numpy/scipy, so portfolios of GPU members work.

Result (needle_d4, GPU): plain portfolio misses the mode (neff 0.3, bias -6.88);
with the Fisher oracle it is unbiased (bias -0.04) at neff 50 / nESS 1494 --
n_eff x194, n_ESS x805.  Portfolio stays unbiased on easy targets (corrgauss3 -0.015).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Makes the normalizing-flow integrator usable on GPU and lets a trained flow be
re-used across ILE instances (the only way NF training pays for itself).

Fixes:
  - GPU crash: the flow (nflows/torch) samples on the host, but rv was pushed to
    cupy and then handed to the host integrand -> 'Implicit conversion to a NumPy
    array' TypeError.  The NF now aggregates on the host consistently (like the
    portfolio fix): rv/lnL kept on host, running estimate forced onto numpy/scipy.
    The CI benchmark (test_mcsamplerEnsemble_extended) now runs to completion on
    GPU instead of crashing in the NF section.
  - save_intg: the final integral reads log_integrand/log_joint_prior/
    log_joint_s_prior back out of _rvs, but they were only stored when
    n_adapt>0 and tempering_exp>0.  Pure flow reuse (n_adapt=0) hit
    KeyError('log_integrand').  Now always accumulated.

New -- flow storage/reuse:
  - save_flow(path): serialize the trained flow (state_dict + architecture meta)
    to a small file (~13 KB).
  - load_flow(path): stage a pre-trained flow into a fresh sampler; the weights
    are applied after the architecture is rebuilt inside integrate_log.  Verifies
    params + box match.  n_adapt=0 -> pure reuse; small n_adapt -> polish epochs.
  - test/integrators/test_NF_reuse.py: train-once/save then reuse in fresh samplers.

Result (corrgauss3, GPU): cold NF trains from scratch in 17s (neff 124); reusing
the saved flow with no training takes 0.3s (x69 faster) at neff 346, unbiased;
a 2-epoch polish gives the best bias. Training amortizes across instances.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Exposes the new integrator capabilities through integrate_likelihood_extrinsic_batchmode
(all default off; no behavior change unless requested):

  --sampler-warmstart-samples FILE  AV: bootstrap the live volume from a prior
        posterior (named-column ASCII; columns matched to extrinsic params by name)
  --sampler-load-state / --sampler-save-state FILE  AV: persist/reuse the adapted
        live-volume state (.npz).  Pointing load+save at one file across a sim-grid
        warm-starts each point from the previous one.
  --nf-flow-load / --nf-flow-save FILE  NF: reuse a trained flow across instances
        (n-adapt 0 = pure reuse, else polish).

Also: the AV (VARAHA) sampler has no update_sampling_prior, so the existing skymap
ResamplingOracle seeding was silently SKIPPED for it.  AV now consumes the skymap
oracle by bootstrapping its live volume from the oracle draws (bootstrap_from_samples).
A warm start only shapes the initial proposal -- it can never bias the integral.

Portfolio-oracle auto-wiring is intentionally NOT added (ILE has no natural Fisher
source for the full extrinsic space, and MCSampler.setup() re-runs oracle.setup(),
which would reset a pre-configured skymap oracle); the mechanism is available via
sampler.oracle_realizations for deliberate use.

Validated: driver compiles + --help lists the options; the exact glue paths
(genfromtxt column-match -> bootstrap; save_state/load_state; save_flow/load_flow)
run end-to-end on the correlated-Gaussian target, all unbiased.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The dev-line integrators (synced in the baseline commit) evaluate the integrand on
a HOST copy of the samples (rv_cpu = identity_convert(rv)).  That suits synthetic /
CI host integrands, but the production vectorized GPU ILE likelihood is DEVICE-native
(wants cupy arrays): feeding it host arrays raises 'Unsupported type numpy.ndarray'
and the whole analysis is skipped.  Verified on a real event (S250114ax) inside the
production container: the container-native AV runs, the host-fed version dies.

Fix: feed the native (device) array to the integrand, matching production; if a
host-only integrand chokes on a cupy array, fall back to a host copy and remember
the choice for the rest of the run.  Both contracts now work:
  * real GPU ILE (device likelihood): runs (matches container-native AV: lnL 3348
    vs 3334 on S250114ax point 0, both neff~1 -- cold AV genuinely struggles on this
    SNR~82 event, independent of this fix);
  * synthetic/CI host integrands on GPU: auto-fall-back, unchanged results.

Note: mcsamplerPortfolio/mcsamplerNFlow gained host-aggregation in Stage 2/3 (good
for host integrands) and need the same device-aware treatment before real-GPU-ILE
use; AV is the production integrator (--sampler-method AV) and is the priority.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Applies the same device-vs-host contract as the AV fix to the two remaining
integrators that Stage 2/3 had made host-aggregating.  The real vectorized GPU ILE
likelihood is cupy-native; feeding it host arrays raises 'Unsupported type
numpy.ndarray'.  Both now evaluate the integrand DEVICE-FIRST (identity_convert_togpu)
and fall back to a host call on a TypeError/ValueError, remembering the choice
(self._integrand_wants_host).  lnL is brought back to the host for the host-side
aggregation (portfolio) / flow math (NF), which is unchanged.

  - mcsamplerPortfolio: fixed the inconsistent 'self.xpy = self.xpy (=numpy)' comment
    -- the driver actually sets self.xpy=cupy, so host-converted draws collided with
    a cupy running estimate.  Now forces host aggregation explicitly AND feeds the
    integrand device-first in both the main loop and the oracle block.
  - mcsamplerNFlow: same device-first integrand eval (the flow itself stays on host).

Host/synthetic path unchanged: test_portfolio_oracle (n_eff x194, unbiased) and
test_NF_reuse (x87 faster reuse, unbiased) still PASS on GPU.  The device path mirrors
the AV fix validated on real S250114ax ILE; a portfolio/NF real-ILE run is the pending
end-to-end check (neither is the production integrator).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The anti-bias foundation for reusing a proposal across DIFFERENT problems (a
neighbouring intrinsic point, a stale breadcrumb, an iteration-to-iteration
handoff).  VARAHA's live volume only contracts, so a warm start seeded at the
wrong location silently contracts there, misses the true peak, and returns a
badly biased integral that still LOOKS converged -- demonstrated:

  cold                       bias_ln -0.04  neff 1841
  wrong seed, NO floor       bias_ln -11.1  neff  215   <-- silent catastrophic bias
  wrong seed, cover_frac 0.1 bias_ln -0.12  neff 1531   <-- safe, ~cold

bootstrap_from_samples() gains cover_frac: it mixes that fraction of uniform
full-box points into the seed cloud, so the seeded live volume is always a
SUPERSET of a cold start.  A mis-placed proposal can then only cost efficiency --
warm coverage always contains cold coverage, so a warm integral is never MORE
biased than a cold one.  Use cover_frac>0 for cross-problem reuse; leave 0 only
for same-problem reuse (an in-run second pass, where the seed provably contains
the peak).  New test/integrators/test_AV_warmstart_safety.py guards this.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ding)

Threads warm-start reuse into the ILE worker, with the anti-bias safeguards up front
(no grid reordering; coverage floor for cross-problem reuse).

L0 auto-rescue (--sampler-warmstart-retry-neff N):  after a cold AV pass, if n_eff < N
(the point stalled on a very sharp / high-amplitude peak), re-run a SECOND pass
warm-started from THIS point's own highest-likelihood samples.  Same-problem reuse -- the
peak provably lies in the seed, so cover_frac=0 and it cannot bias.  Validated on the real
S250114ax event: cold n_eff 1.0 -> auto-rescue -> n_eff 30 (target), one flag, automatic.
This is the high-value case (VARAHA otherwise cannot resolve a ~1e-11-volume peak).

L1 sequential warm-start (--sampler-sequential-warmstart[-cover-frac]):  when a worker
analyzes several intrinsic points, seed each point's extrinsic AV from the previous
point's converged high-L samples.  Points are processed IN GIVEN ORDER (no reordering),
so a truncated/failed worker still drops a spatially-unbiased subset -- the correctness
concern RO'S flagged.  A coverage floor (cover_frac, default 0.1) makes a poorly-matched
cross-point transfer degrade to cold rather than bias.  Runs safely on real 2-point
S250114ax (captures gracefully, unbiased); cross-point speedup is target-dependent and
modest where cold AV already carves quickly (the safety floor also dilutes the seed) --
L0 is where the large, reliable win is.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Lightweight substrate for the 'breadcrumb' warm-start strategy: an intrinsic-space
field of extrinsic proposals.  After an ILE iteration each converged point contributes
its high-likelihood extrinsic samples keyed by intrinsic lambda; the next iteration's
workers query the nearest (whitened-metric) proposal and warm-start from it.  Provides
add / nearest / warm_seed_for / save / load (small npz, Condor-transfer friendly) and a
build_field_from_run_outputs() aggregator.

Does NOT touch the DAG: pipeline wiring (a post-iteration build node + a driver
--extrinsic-proposal-field query hook) is the next step and mirrors the existing calmarg
extrinsic-breadcrumb plumbing.  A field entry can only shape p_s, so a stale/mismatched
proposal costs efficiency, never correctness (callers still pass cover_frac>0).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
… margin

Wires the ProposalField 'breadcrumb' end to end, with two safety knobs for the
cross-iteration handoff so a stale/mismatched field can only cost efficiency, never bias.

  - AV.bootstrap_from_samples(inflate=): HANDOFF SAFETY MARGIN.  Widen the imported
    seed cloud by this factor about its mean before building the grid, to cover the peak
    SHIFT between the neighbouring intrinsic point that produced the proposal and this
    one.  cover_frac is the coarse net for gross mismatch; inflate is the fine margin for
    a modest shift.
  - proposal_field.lambda_from_P(): canonical intrinsic key [m1,m2,s1x..s2z].
  - bin/util_BuildProposalField.py (PRODUCER): post-iteration node that scans an
    iteration's per-point ILE --save-samples outputs + grid, keeps each point's
    high-lnL extrinsic samples (converted to the sampler's cosine-sampler coords), and
    writes a ProposalField .npz.
  - driver --extrinsic-proposal-field FILE (+ -cover-frac, -inflate) (CONSUMER): each
    intrinsic point warm-starts from the field's nearest entry.  Takes priority over L1
    sequential.  No grid reordering; a missing/far entry degrades to cold.

Validated: inflate widens the seed; full round-trip (build field from point A -> save ->
load -> query at a nearby lambda -> warm-start point B with cover_frac+inflate) is
unbiased (bias +0.04, neff 883).  Remaining pipeline step is the DAG node that runs the
producer after each iteration (mirrors the calmarg extrinsic-breadcrumb node).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
… field

Reframes L3 around operational reality (per RO'S): --save-samples burns disk fast, and
random pilot selection endangers the CIP grid fit.  So we neither dump the whole grid nor
guess -- we cherry-pick the BEST point(s) by lnL after a cheap iteration 0.

  - util_PickPilotPoint.py: after iteration 0 (per-point marginal lnL only, no sample
    dump), pick the top-k points by lnL and emit them as a small pilot sub-grid.  Never
    random -- ranked by lnL.
  - --sampler-warmstart-samples now applies the handoff safety margins by default
    (--sampler-warmstart-cover-frac 0.1, --sampler-warmstart-inflate 1.5): the pilot seed
    usually comes from a different (best) point, so cover_frac + inflate keep a shifted
    peak from biasing while staying far tighter than cold.
  - util_BuildProposalField.py docstring repositioned: whole-grid save-samples is a disk
    disaster; the multi-pilot field is only for the rare few-well-separated-pilots case,
    running save-samples on those few points only.

Standard flow: cheap iter-0 ILE -> util_PickPilotPoint (best point) -> 1 pilot ILE with
--save-samples (~tens of KB) -> warm-start the rest via --sampler-warmstart-samples with
the safety margins.  A single pilot node, chosen not guessed.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…nder a wrong member

Problem
-------
mcsamplerPortfolio pooled draws from several member samplers but formed the
integral with a STRATIFIED estimator: each pooled sample kept its OWN member's
sampling density p_s in weight = L*prior/p_s, and members were re-weighted by
their per-member Kish n_ess (portfolio_default_weights). This is UNSAFE. A
warm-started AdaptiveVolume member that contracts onto a DECOY (a flat, wrong
region) draws near-uniform-weight points, reports a HIGH per-member n_ess, is
driven to weight ~1, and starves the broad covering member (e.g. GMM) to the ~1%
floor. The stratified estimator then has expectation sum_m w_m * Z_m, so it is
biased LOW by ~ln(covering-floor): E[I] ~ w_AV*Z_decoy + w_GMM*Z_true with
w_AV->1. A broad coverage member canNOT rescue it. This is real bias, not
variance.

Fix: multiple-importance-sampling balance heuristic (Veach & Guibas)
--------------------------------------------------------------------
Every pooled sample is now weighted by the deterministic-MIXTURE density
    q_mix(theta) = sum_m frac_m * q_m(theta),   frac_m = n_drawn_m / n,
evaluated AT that sample, where q_m is member m's own sampling density and
frac_m is the actual per-chunk sampling fraction. Then
    E[I_hat] = \int q_mix * (L*prior/q_mix) dtheta = \int L*prior dtheta,
UNBIASED for ANY member weights as long as the mixture covers the peak. A broad
member with even a tiny weight guarantees q_mix>0 at the true mode, so a
wrongly-contracted member can no longer bias the result. Draw allocation
(n_ess-based member weights) is unchanged; only the ESTIMATOR denominator moved
from per-sample own-density to q_mix.

Changes
-------
* mcsamplerAdaptiveVolume.sampling_density(X): pointwise q at arbitrary points --
  1/(n_bins*prod(dx)) inside the live volume (union of self.binunique bins,
  == the constant 1/(V*prod(dx0)) integrate_log already reports), 0 outside.
  Read-only; AV's own integrate_log is untouched. Handles cupy/host.
* mcsamplerEnsemble.sampling_density(X): the GMM mixture density at X, exactly
  the per-dim-group product MonteCarloEnsemble stores as sampling_prior_array
  (gmm.score, or uniform 1/vol for an unfit group). Read-only.
* mcsamplerPortfolio: draw() records the chunk's active members and true
  sampling fractions; integrate_log() builds q_mix from them and uses
  log(q_mix) for log_integrand / log_weights / stored log_joint_s_prior. If any
  active member lacks a density (e.g. AC-histogram members) it falls back to the
  legacy stratified per-member p_s, so those portfolios are unchanged. New kwarg
  portfolio_use_mixture_density (default True) toggles the estimator.
* mcsamplerEnsemble.update_sampling_prior: make backend-consistent (host<->cupy)
  so a GMM member can adapt from the portfolio's host history on GPU -- it
  previously crashed there (self.xpy.min([list]); host slice into a cupy row).

Test (test/integrators/test_portfolio_balance_heuristic.py)
-----------------------------------------------------------
AV seeded at a decoy plateau + a broad GMM, on a correlated Gaussian where cold
AV converges. Decoy AV is driven to weight 0.99, GMM starved to 0.01. On GPU:
  OLD stratified estimator : lnI - lnZ = -2.72   (badly biased low)
  NEW q_mix estimator      : lnI - lnZ = +0.10   (unbiased, ~few %)
  control (cold AV + GMM)  : lnI - lnZ = -0.10   (no regression)
Stable across seeds (new bias within +-0.10). test_portfolio_oracle (AC members
-> legacy fallback) still PASSES on GPU; standalone AV/GMM and AV warm-start
safety unchanged.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
_build_grid_from_points computed binidx = ((pts-lo)/dx).astype(int) without clipping,
so a point on the upper box edge maps to binidx == nbins (out of range).  Those
out-of-range bins entered binunique, making V = n_bins/prod(nbins) exceed 1 -- an
invalid fractional volume -- and draw_simple would then sample outside the box.  This
crashed warm-start from a WIDE seed (a full PE posterior + cover_frac spanning the box):
observed V=1.309 -> 'zero-size array to reduction cupy_max' / FAILED ANALYSIS.

Fix: clip bin indices to [0, nbins-1].  Verified on the S250114ax PE-oracle worker: the
production-bilby-PE warm start now builds a valid grid (V=1.0) and the point converges
(n_eff 61 vs 18 for L0-rescue alone).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
A cover_frac coverage floor mixes uniform full-box points into the seed.  Because
_build_grid_from_points estimated the extent (hence nbins) from the FULL cloud, those
uniform points widened the extent to the whole box and collapsed the grid to one bin per
dim (V->1), throwing away the proposal's concentration -- so a wide PE-posterior warm
start with any cover_frac degraded to a cold start.

Fix: _build_grid_from_points gains resolution_pts; bootstrap_from_samples passes the
CORE (the proposal, pre-cover) for the extent/nbins while binning core+cover into that
resolution.  Coverage points now land in scattered fine bins (still guaranteeing a
superset of cold) without coarsening.  S250114ax PE seed: cover_frac=0.05 now gives
V=0.10 nbins=6 (concentrated + floored) instead of V=1.0 nbins=2.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
sample_from_bins looped over bins in Python (list-comp + vstack, one uniform draw per
bin) -- O(n_bins) per chunk.  Fine for cold AV (coarse early grids) but a hard bottleneck
for a concentrated warm start: a full-PE-posterior seed occupies thousands of bins, and a
single ILE worker then timed out drawing.  Vectorized: repeat each bin's lower corner
ninbin[k] times (host np.repeat) + one (N,ndim) backend uniform draw -> O(N), no bin
loop.  ~4800-bin grid, 200k draws: was a multi-minute stall, now 0.007 s.  Cold AV result
unchanged (corrgauss3 bias/neff identical); output stays on the active (cupy) backend.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
@oshaughnessy-junior

Copy link
Copy Markdown
Author

Status on the four blockers, now that #45 is merged into rift_O4d (this PR's head 87d66b75).

# Blocker Status
3 Portfolio warm starts / L0 rescue do not apply AV's seeded grid in the active path ResolvedAV._apply_warm_state() installs the seed on the draw path; the merge gate now has a white-box case that fails if it goes inert
4 --sampler-sequential-warmstart disconnected because _rvs is cleared before the next seed is captured Resolved — the seed is captured inside analyze_event before the wipe (_SEQ_WS_PENDING)
2 MC-error replicas are not cold or adaptation-independent; reset does not fully clear fitted GMM/portfolio state Fixed, in two parts: the reset itself landed in #45 (clear_warm_state, which replays each member's original setup arguments and clears the trained proposal), and the replica loop now calls it — oshaughnessy-junior#47
1 AV warm starts can bias evidence because cover_frac adds sparse points instead of full-prior support Confirmed and quantified; docstring corrected, backstop rule fixed — oshaughnessy-junior#47. See the caveat below

On [1], the measurement changed my understanding of it, so stating it plainly rather than claiming a fix:

  • cover_frac genuinely does not give full support. At d=6, even cover_frac=0.9 covers 2.9% of the prior box. The claim in the docstring was false and is now replaced by the measured table.
  • But the failure it implies is not reproducible in the default [AV, GMM] portfolio: with a deliberately displaced seed, |lnZ bias| <= 0.05 in every d=4/d=6 run. The GMM member's uniform defensive component (gmm_defensive_frac=0.05) plus Gaussian tails keep q_mix positive everywhere. That — not cover_frac — is what has protected production, and it was written down nowhere.
  • The failure is reproducible with an all-compact-support portfolio: all-AV, same displaced seed, -1.0 to -6.8 nats.
  • Holding a member cold restores the stated invariant and is free with a correct seed, but it does not rescue a mismatched one (-1.1 to -4.2 nats even cold, n_eff 3-9). Coverage in principle is necessary, not sufficient.

So [1] is closed as: the false guarantee removed, the real one documented and enforced structurally. The residual risk is a mismatched seed, which is detectable rather than preventable — a support-escape statistic over the mixture is being prototyped and validated separately. I would rather land the honest correction now than imply a coverage knob fixes it.

Shape gate for #47 vs rift_O4d @ 87d66b75: COMPARE_EXIT=0, 0 blocking regressions; flag-ON probe 0 regressions.

oshaughnessy-junior and others added 11 commits August 6, 2026 07:45
…ass) + its ROC

Adds, in dict_return and a one-line print, the fraction of total posterior weight carried by
samples each member could not have drawn:

    escaped_mass[m] = sum_{i : q_m(x_i)==0} w_i / sum_i w_i

plus its first-chunk variant, the per-chunk history, and the soft comparator weight_share[m].
Free: integrate_log already evaluates every member's density at every pooled sample to build
q_mix and keeps frac_m*q_m in _chunk_mix_parts, so this is one compare and one masked sum per
member per chunk.  Strictly off-path -- verified bit-identical lnZ/var/n_eff with the reduction
disabled (on the all-AV arm, the only portfolio shape that is bit-reproducible at all: an
[AV, GMM] portfolio does not reproduce itself run to run, with or without this code), plus a
structural test that it mutates nothing outside its own accumulators.

ROC on the truth-known testbed (escaped_mass_study.py / escaped_mass_report.py; 1440 runs,
20 independent target seeds per cell, d=4 and d=6, seed displaced 0..4).  Honest verdict: this
is a WARM-START QUALITY monitor, not a lnZ-bias alarm.

 * The cumulative statistic has NO usable floor -- a correctly-placed seed's converged live
   volume legitimately excludes median 0.51 (d=4) / 0.80 (d=6) of the posterior weight.  Use the
   FIRST-CHUNK statistic, whose offset=0 floor is 6e-6 / 2.5e-7 (max over 20 seeds 2.3e-4/6.8e-3).
 * With a COLD broad member present, first-chunk escaped_mass at a 0/20-false-positive threshold
   gives TPR 0.70/0.95/1.00 (d=4) and 0.10/0.50/0.80 (d=6) at displacement 1/1.5/2 -- beating
   1-weight_share, which reaches only 0.20/0.05 at displacement 2.
 * It is EXACTLY 0.000 in 320/320 runs when every member is seeded from the same cloud (the
   default, since mcsamplerEnsemble also has bootstrap_from_samples) or the portfolio is all-AV.
 * The apparent perfect cumulative detector under early stopping is a RUN-LENGTH artifact: at
   fixed budget the offset=0 floor moves 4.9e-5 -> 0.524 and the AUC collapses to 0.47-0.58.
 * Structural limit: escaping weight is only observable when another member covers the
   complement -- exactly when the balance heuristic already keeps lnZ unbiased.  Where
   displacement does bias lnZ (all-AV), no variant reaches a usable threshold.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…ptation for replicas

Both PR #47 findings.

[P1] Full-support detection could discard the only real backstop.  The check
defaulted un-annotated samplers to full-support and counted a nominally broad
member even after it had been RANGE-RESTRICTED.  Reproduced: [unrestricted AV,
restricted GMM] reports _full_support_members == [0], yet the restricted GMM set
_has_broad, so member 0 was warm-started and contracted to V=0.095 -- nothing
covering the prior box, which is the silent low bias this mechanism exists to
prevent.

Now: has_unbounded_support defaults FALSE and must be declared; only members
that are BOTH declared broad AND still full-range (per _full_support_members)
count.  mcsamplerEnsemble declares it as a property keyed on
gmm_defensive_frac > 0 -- the uniform defensive component is the actual
guarantee, since Gaussian tails underflow to exactly zero far from the mode, so
gmm_defensive_frac=0 correctly reports False.  AV declares False.

[P1] Replicas retained portfolio-level adaptation.  clear_warm_state() rebuilds
the members but not the portfolio's own learned state: portfolio_weights,
portfolio_quality and its observation counts, portfolio_probe_ptr,
portfolio_draw_iteration, breakpoint progression and per-member n_ess histories.
Replicas therefore scheduled themselves from what earlier replicas learned, so
the between-replica scatter -- the entire quantity being measured -- still
understated the error.  Adds reset_adaptation(), which does clear_warm_state()
plus a restore of every field to its POST-SETUP value (snapshotted in setup, so
an explicit initial weighting or breakpoint schedule is preserved rather than
replaced by a hard-coded guess).  The replica loop calls it.

Four tests, each verified to fail when its defect is reintroduced, including the
requested assertion that all portfolio-level state is identical across replicas.

Noted while testing, NOT fixed here: mcsamplerEnsemble.bootstrap_from_samples
rebuilds its integrator without the caller's gmm_defensive_frac, so a warm start
silently restores the 0.05 default.  Same class as the gmm_dict config loss
fixed in #45; filed separately rather than widened into this PR.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Every gate verdict is a hard threshold on a stochastic quantity, so a cell near
a threshold flips on realization alone and can block a merge that changed
nothing.  Adds confirm_regressions.py: re-runs only the disputed cells, in BOTH
arms, at several fresh run seeds, and blocks only if the candidate is worse in a
majority (ties count as not-worse -- the burden of proof is on the claim that
the candidate broke something).

Validated on the live case that motivated it.  `GMM mix_d6_n3_s303` was reported
as a blocking REGRESSION in two consecutive full runs, base 119 vs candidate 66,
which looked reproducible rather than flaky.  At 5 fresh seeds the arms were
BIT-IDENTICAL -- 93/93, 80/80, 119/119, 95/95, 96/96 -- so the branch does not
reach that row at all, and 4 of the 5 seeds starve: the cell sits on the
n_eff=100 floor and its PASS at the default seed is the lucky draw.  The
apparent regression was an artifact of where the job landed in the worker pool.

Deliberately NOT fixed by seeding the samplers.  Independent copies that
localize differently are our working detector for support and mode-collapse
failures; pinning every fit to one seed would silence it and make N production
copies no better than one.  The right answer is to ask the question again at
fresh seeds, which is what this does.

Separately worth review: that cell starving 4 of 5 seeds means it is
mis-budgeted for the strict set.  Not changed here -- the strict list and budgets
are shared with other people's work.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…ore claiming coverage

PR #47 follow-up.  has_unbounded_support read gmm_defensive_frac > 0 as a
guarantee, but add_defensive_component() was only called by fit_gmm_adaptive.
The fixed-component fit paths (mcsamplerEnsemble and MonteCarloEnsemble) called
GMM.gmm(...).fit(...) directly, and gmm_adaptive defaults to None -- so the
DEFAULT configuration requested a defensive component and never installed one,
and the portfolio then contracted its AV member believing it was covered.

Worse than a missing component: gmm.score() FLOORS its return at 1e-300, so the
member always LOOKS like it has density everywhere.  That floor is a guard
against log(0), not coverage -- a sample landing there carries weight
L*p/q ~ 1e300 and would wreck the estimate rather than support it.  Measured, a
fixed-component fit to a tight cloud returns exactly the floor at the far corner
of the prior box for every d >= 4:

    d=2  2.9e-273 -> 1.3e-04     (with the defensive component)
    d=6  1.0e-300 -> 7.6e-10

This also corrects a claim I made in PR #47: "the GMM member's unbounded support
is what has protected production".  It is not unbounded support, it is a
numerical floor.  The measured |lnZ bias| <= 0.05 with a displaced seed means
the far region was never sampled in those runs, not that it was covered.

Fixes: add_defensive_component sets model.defensive_frac as a verifiable marker;
both fixed-component fit paths now install it; has_unbounded_support inspects the
installed models and reports False if ANY trained group lacks it, or if
gmm_defensive_frac <= 0.  While untrained it trusts the config -- sound only
because every fit path now installs the component, which
test_every_fit_path_installs_the_defensive_component pins.

Note this changes the DEFAULT GMM proposal: fixed-component fits now carry a 5%
defensive component (the documented gmm_defensive_frac default, previously
inert).  That is the intended behaviour per its own docstring and it bounds the
importance weights, but it is a shared-sampler change -- gate results reported
separately.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…pay for it

The previous commit installed it on every fixed-component fit, changing the
DEFAULT proposal for all users.  The shape gate priced that and it is too
expensive: a 5% broad component spends 5% of draws where the likelihood is
negligible, and n_eff falls systematically at higher dimension --

    GMM d6_n1_s101  2036 -> 1565      d6_n3_s202    39 ->  15
    GMM d6_n3_s303   119 ->   75 *    d8_n1_s303   448 -> 210 *
    (* blocking regressions; d<=4 rows were unchanged within scatter)

So the install is now behind `gmm_defensive_all_paths`, default False.
mcsamplerPortfolio sets it on its own members in setup(), because a portfolio
member may be the mixture's ONLY full-support component and the guarantee has to
be real there.  A standalone GMM user gets exactly the previous behaviour and
none of the cost.

has_unbounded_support now also returns False when neither install path is active
(no opt-in and gmm_adaptive off), so it cannot promise coverage that the
configuration will not build.

Verified: standalone GMM reports all_paths=False and declares no coverage;
a portfolio member reports True and declares coverage.  20 unit tests pass.
Gate re-run for this shape is pending -- the numbers above are from the
default-changing variant, which this commit supersedes.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…_diagnostic

Off-path warm-start support-mismatch diagnostic (escaped_mass) + its ROC assessment
…led reruns, and is wired in

Three P1s from review of #49.

[P1] The confirm step recognised only PASS -> non-PASS.  The comparator also
blocks REGRESSION(metrics) -- JS, pull, width, correlation, evidence bias or
n_eff worsening beyond tolerance -- and those rows produced "no blocking
regressions to confirm" and exit 0.  Our own gate v11 had exactly such a row
(GMM d8_n1_s303, n_eff 448->210), which would have been waved through.

Fixed structurally rather than by adding a second branch: compare_shape_results
now exposes classify() / is_blocking() / blocking_keys() as the SINGLE definition
of a regression, and confirm_regressions imports them.  Two copies of that logic
will always drift apart; there is now one.  Verified the refactor reproduces the
v11 verdict exactly (2 blocking).  The confirm step now sees both rows where it
previously saw one.

[P1] Failed reruns were silently skipped, so worse == same == 0 read as "not
confirmed".  Now: a candidate that produces no record where the base did counts
AGAINST the candidate (crashing is worse than passing, not missing evidence); a
verdict requires --min-valid usable pairs, defaulting to all seeds; and too few
valid pairs is INCONCLUSIVE with a nonzero exit, never a silent clear.

[P1] Confirmation was documented but never invoked.  compare_shape_results gains
--confirm-base-checkout / --confirm-cand-checkout / --confirm-repeats and returns
the confirmed verdict as its exit code, so the comparison workflow enforces it.
Without those flags it still exits 1 on a blocking row, and now says explicitly
that the row was NOT confirmed rather than implying it was.

Adds test_confirm_regressions.py (5 checks, all on the dangerous direction -- the
ways a confirmation can wrongly CLEAR a real regression).

CI on this PR is red at "Set up job" on 4 jobs; the same runner-provisioning
failure hits #47, which shares no files with this change, and no failing job
reaches a step that executes repository code.  Infrastructural, not from here.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The previous commit landed without this hunk (a heredoc escaping slip), so the
runner still advertised the standalone command as the whole procedure.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
classify() returned ONLY-IN-BASE when the candidate omitted a row, which
is_blocking() did not treat as a regression -- so the row never reached
confirmation and the gate exited 0.  A candidate crashing before emitting its
first result could bypass the fail-closed rerun logic entirely.

Now REGRESSION(missing-in-candidate).  Confirmation re-tests such a row using
whichever record exists for the cell spec, so the candidate is actually re-run
rather than written off.  The reverse (a NEW row in the candidate) still does not
block.  Two tests added.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…ust at setup

Two PR #47 findings: the defensive component was installed correctly and then
lost again during normal operation.

[P1] Warm bootstrap disabled the portfolio's opt-in.  bootstrap_from_samples()
calls setup() to rebuild the integrator as one full-dim group, and did so bare --
resetting gmm_defensive_all_paths to False and refitting the warm GMM with no
defensive component, AFTER the portfolio had already decided on the strength of
that flag that it was safe to contract its AV member.  It now carries the
defensive config forward from the existing integrator.

[P1] Updates absorbed the defensive component.  _merge() blends component i with
the fitted component order[i] for every i in range(self.k) and never consults
self.adapt, so the broad component -- marked adapt=False precisely so it would be
left alone -- drifted toward the fitted cloud on every update while
defensive_frac stayed set and has_unbounded_support kept reporting coverage.
Measured, far-field density fell 3.09e-07 -> 2.04e-08 after a single update
cycle.  update() now detaches the defensive component, updates the real ones and
reinstates it; density is constant at 3.09e-07 across 4 updates with the
component's weight, mean and covariance unchanged.

The tests previously only checked initial installation.  The two added here
exercise the warm-bootstrap path and repeated updates, and assert on the actual
far-field density and on the defensive component's own parameters rather than on
the marker attribute -- the marker is exactly what stayed true while the
guarantee evaporated.  Both verified to fail when their defect is reintroduced.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Conflict was structural, not semantic: both sides added a method immediately
before bootstrap_from_samples, and the two bootstrap signatures differed.
Resolution keeps BOTH bodies -- the merged diagnostic and reset_adaptation --
with the single signature that carries keep_backstop_cold.

Verified after resolution: 22 portfolio tests and all 5 of the diagnostic's own
tests pass.
@oshaughnessy-junior

Copy link
Copy Markdown
Author

Daily triage refresh (2026-08-07): PR #161 is still MERGEABLE / CLEAN, and the shape-recovery gate passed with zero blocking regressions, but the current head 0b2e56c7 remains blocked on correctness:\n- MC-error replicas are not cold because the loop never calls clear_warm_state() after warm-state-aware samplers are reused.\n- Evidence and posterior export can describe different runs because posterior/fairdraw samples come only from the highest-n_eff replica.\n- L0 rescue can exclude modes because it rebuilds AV support from a single peak with cover_frac=0.\n- JAX validation is still incomplete here because jax is unavailable locally and the new JAX scripts do not provide pytest functions.\n\nSee the baseline note for the exact commands/results from the isolated review.

…eeds before they block

gate: confirm blocking regressions at fresh seeds before they block
…plicas, and a lifecycle-durable coverage guarantee

Warm-start coverage honesty + genuinely cold MC-error replicas (PR 161 blockers 1 and 2)
@oshaughnessy-junior

Copy link
Copy Markdown
Author

rift_O4d now carries all four blockers' fixes. Recap of where each stands, with the merges that resolved them:

# Blocker Resolved by
3 Portfolio warm starts / L0 rescue do not apply AV's seeded grid in the active path #45 (AV._apply_warm_state), plus a white-box gate case that fails if it goes inert
4 --sampler-sequential-warmstart disconnected by the _rvs wipe #45 (_SEQ_WS_PENDING)
2 MC-error replicas not cold or adaptation-independent #45 (clear_warm_state) + #47 (reset_adaptation, and the replica loop calls it)
1 cover_frac adds sparse points instead of full-prior support #47

On [1], the correction is worth restating because the original diagnosis was right but the mechanism was not what either of us assumed:

  • cover_frac genuinely does not give full support -- at d=6, cover_frac=0.9 covers 2.9% of the prior box. The docstring claiming a warm start "can never be MORE biased than a cold one" was false and is replaced by the measured table.
  • What had actually been protecting production was not the GMM's unbounded support, as I stated earlier in this thread. gmm.score() floors at 1e-300 -- a log(0) guard, not coverage. A sample landing there carries weight L*p/q ~ 1e300. The measured |lnZ bias| <= 0.05 with a displaced seed meant the far region was never sampled, not that it was covered.
  • The real guarantee is the defensive component, which was requested by default and never installed on the fixed-component fit path (add_defensive_component was reachable only from fit_gmm_adaptive, and gmm_adaptive defaults off). It is now installed on every fit path, opt-in and scoped to portfolio members -- a global default change measured 2 blocking regressions (n_eff 448->210 at d8), so standalone GMM users pay nothing.
  • And the guarantee now survives the sampler's lifecycle: warm bootstrap no longer resets the opt-in, and update() no longer lets _merge() absorb the defensive component (far-field density was collapsing 3.09e-07 -> 2.04e-08 per update cycle while the marker still reported coverage).

Residual risk, stated plainly: a badly mismatched seed remains neither prevented nor reliably detected. Coverage in principle is necessary, not sufficient -- an all-AV portfolio with a displaced seed still gave -1.1 to -4.2 nats even with a cold backstop, because a uniform member finds a sharp 6-D peak too rarely to carry the integral. A support-escape statistic was prototyped and assessed (#48) and does not close this: it is observable only when another member covers the complement, which is exactly when the balance heuristic already protects lnZ. It is merged as an off-path diagnostic, not a gate.

Also merged: #49, which requires a blocking gate regression to reproduce at fresh seeds before it blocks -- after a "reproducible" 119-vs-66 regression across two full runs turned out bit-identical between arms at 5 fresh seeds.

Test state on the merged tip: 22 portfolio invariants, 7 confirm-on-fail, 6 sequential-GMM gate checks, 5 escaped-mass -- all passing. Every regression test added in this cycle was verified to fail when its defect is reintroduced.

This PR should be ready for re-review.

oshaughnessy-junior and others added 13 commits August 7, 2026 02:49
Add test/waveforms/test_uv_symmetry.py, a waveform-level test that builds
the U = <h_A|h_B> and V = <h_A*|h_B> mode-cross-term matrices via
factored_likelihood (recomputing every element independently, bypassing the
internal symmetrization shortcut) and verifies, over semi-random parameters
looped across an active waveform list:

  (1) U Hermitian:            U[A,B] = conj(U[B,A])
  (2) U diagonal real, > 0
  (3) V complex-symmetric:    V[A,B] = V[B,A]
  (4) reflection/parity (aligned-spin only):
      V[(l,m),B] = (-1)^l U[(l,-m),B],  from h_{l,-m} = (-1)^l conj(h_{l,m})

Definitional checks (1)-(3) run on aligned and precessing models; the
reflection check (4) only on aligned-spin models. Models the local lalsuite
build cannot generate are skipped with a reason, not failed.

Includes one deliberately-failing placeholder
(test_full_nonlinear_reflection_symmetry_left_as_exercise) for the full
non-linear reflection algebra: assert False, "this failure is left as a test".

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
… setup() remember its kwargs

Three items from the PR 161 review.

[P1] The exported posterior could disagree with the reported evidence.  lnZ was
the linear mean over K replicas but samples were exported from whichever replica
had the largest n_eff.  n_eff measures weight CONCENTRATION, not coverage, so a
mode-collapsed replica scores HIGHEST and was exactly the one selected -- the
same rule this work already documented as wrong (the highest-n_eff copy in an
earlier study was the most biased, 11 nats low).  Demonstrated: with two broad
replicas (Kish n_eff 960, 1023) and one collapsed (3990), selection picks the
collapsed one.

Replicas are now POOLED.  The estimator is Zhat = (1/K) sum_k (1/n_k) sum_i w_ki,
so pooled weights are w_ki/(K n_k), which is the importance weight against the
real pooled proposal q'_ki = q_ki * K * n_k -- "pick a replica uniformly, then
one of its n_k draws".  Folding K*n_k into log_joint_s_prior states that density
where every downstream weight computation already accounts for it.  Verified the
pooled samples imply the combined lnZ to 1e-9, including unequal replica sizes.
Reported n_eff is now the pooled Kish value rather than the sum over replicas,
which overstates whenever the replicas disagree -- the case they exist to detect.

[P2] The L0 rescue's warm pass has truncated support.  Finding one peak does not
prove all modes are represented, and cover_frac cannot repair it (2.9% of the box
at d=6).  It now WARNS when the warm lnZ falls well below the full-support cold
pass, and points standalone AV users at the portfolio, which has a defensive GMM
member for this.

Deliberately NOT "fixed" by pooling cold+warm, and the reason is pinned in a
test: averaging Z is unbiased only when every term is.  Measured on a bimodal
target whose seed caught one mode -- true lnZ -4.6052, cold -4.6265, warm
-5.3009 (bias -0.696 = -log 2), pooled -4.9079 (bias -0.303 = log 0.75).  Pooling
is better than warm-only and still wrong.  A correct combination needs the
balance heuristic with q_mix at every sample, and the cold pass's proposal is not
recoverable after it adapted.

[general] setup() now remembers its kwargs and re-applies them, instead of a
fourth targeted patch.  It is re-run by bootstrap_from_samples and by the
portfolio's member reset, and each rebuilt the integrator from only that call's
kwargs -- silently dropping gmm_dict, gmm_defensive_frac and
gmm_defensive_all_paths, each of which was found as its own P1.  An explicit
argument still wins; setup_forget=True starts clean.

Adds test/integrators/test_replica_pooling.py (5 checks).  All suites pass:
22 portfolio, 5 pooling, 7 confirm-on-fail, 5 escaped-mass, 6 seq-gmm.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…k; force the full-dim warm GMM

Three PR #51 findings.

[P1] Standalone AV was knowingly biased and a warning did not fix it.  The warm
pass is an estimate over TRUNCATED support: the seed provably contains the peak
the cold pass found, and nothing about what it did not reach.

Deliberately NARROW change.  The rescue still runs exactly as before -- removing
it by default would strip the high-SNR n_eff rescue from every existing AV run, a
certain regression traded against a bias so far shown only on an adversarial
synthetic.  What changes is only the case with POSITIVE EVIDENCE of lost mass:
the cold pass had full support, so if the warm lnZ lands more than
--sampler-l0-rescue-reject-dlnZ (0.5) below it, the cold result is kept and the
reason printed.  --sampler-l0-rescue-accept-truncated restores the old behaviour.
Detection is imperfect -- a missed mode need not produce that ordering -- so this
NARROWS the failure rather than closing it, and says so.

Not "fixed" by giving AV its own defensive component.  That needs per-sample
mixture densities, which is mcsamplerPortfolio: q_mix, defensive members and
their coverage bookkeeping already exist there.  A second implementation of the
same mathematics is how the bugs in this very review arise -- a capability flag
that lied, a defensive component absorbed by an update, config dropped on
re-setup -- when one of two parallel paths is updated and the other is not.  The
architectural answer is that coverage-critical runs use the portfolio; AV stays a
fast single-proposal sampler with the limitation documented at its call site.
Also corrects the --sampler-warmstart-retry-neff help, which claimed the rescue
"cannot bias the result".

[P1] Pooling assumed raw importance records.  In production _rvs may already be
thresholded or fairdraw-resampled, so sum_i w_ki over retained rows is not
Z_k * n_k and a 1/n_k rescale mis-weights the block (a fairdraw record is already
posterior-resampled, hence weighted twice).  Given the reported per-replica lnZ,
each block is now renormalized to contribute exactly Z_k/K -- correct for raw,
pruned or resampled rows.  The new test prunes a replica and asserts both that
the fix is exact and that the naive rescale would have been wrong.

[P2] Setup memory defeated the full-dimensional warm GMM.  A remembered explicit
gmm_dict survived bootstrap_from_samples, so correlate_all_dims=True no longer
rebuilt the single full-dim group that path exists to build (capturing
sky<->phase<->distance correlations at high SNR).  It now passes gmm_dict=None
explicitly; verified [(0,), (1,2)] -> [(0,1,2)] with coverage settings preserved.

All suites pass: 22 portfolio, 6 pooling, 7 confirm-on-fail, 5 escaped-mass,
6 seq-gmm.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…cks-junior

test/waveforms: U,V mode-cross-term symmetry checks
Script-only companions to test_uv_symmetry.py implementing the exact GR
parity identity h_lm[s_xy -> -s_xy] = (-1)^l conj(h_{l,-m}):

- parity_check_hlm.py: mode-level check over superkick-like, generic
  precessing, and nonprecessing configurations (perturbed superkicks; the
  exact degenerate point is a convention branch point for several models).
- uv_parity_diagnostics.py: the same physics on the U/V cross-term matrices,
  in three tiers: D1 definitional (U Hermitian, V symmetric; exact for any
  waveform), D2 nonprecessing-point single-generation check
  V = (-1)^l U[m-flip], D3 reflected-pair relations valid for any
  configuration. Clean models sit at <=1e-10 relative Frobenius residual,
  violating models at >=1e-4.

Deliberately NOT collected by pytest: run against currently released
precessing implementations these checks FAIL, correctly (NRSur7dq4 percent
level generically, tens of percent in superkick subdominant amplitudes and
2% spurious (2,+-1) aligned-spin asymmetry; SEOBNRv5PHM antisymmetric modes
few x 1e-4; SEOBNRv4PHM 0.3-1% with frame origin; IMRPhenomTPHM and
symmetric SEOBNRv5PHM parity-exact). Their role is rapid assessment of
upstream waveform regressions when adopting a model version or interface.
Marginal-likelihood impact demo (NRSur7dq4 crosses |dlnL_marg|=1 at
rho~37, edge-on near-superkick; SEOBNRv5PHM+antisym at rho~160) lives in
RIFT_roboto_paper demos/waveform_symmetry.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…e cold dict_return

[P1] Fairdraw replicas were weighted twice.  Renormalizing each block to Z_k/K
fixes its SCALE but not its SHAPE: fairdraw samples were already drawn in
proportion to their own importance weights, so reusing those weights applies them
a second time and the block follows w^2 rather than w.  The driver enables this
in production (igrand_fairdraw_samples comes from --fairdraw-extrinsic-output),
and the samplers resample _rvs in place at the end of integrate_log while leaving
the original weight columns, so every replica record is affected.

Restructuring the run order does not help: the FIRST pass has already fairdrawn
before the replica trigger is evaluated.  So a fairdraw block now contributes
what it actually is -- an equal-weight draw from its own posterior -- with
constant weights within the block summing to Z_k/K.

The added test resamples a record in proportion to its own weights (as the
samplers do), then asserts the pooled block has constant within-block weights AND
that the un-handled path retains a w^2 spread with less than half the pooled
n_eff.  The previous test covered pruning only, which is a different failure.

[P2] Rejecting the warm AV pass left the warm dict_return in place, so khat,
block scatter, ESS, the confidence interval and the replica trigger downstream
all described a run whose result had been discarded.  The cold dict_return is now
saved and restored with res/var/neff/_rvs.

Also corrects the --mc-error-replicas help, which still claimed the export "uses
the best-converged (highest-neff) replica".

All suites pass: 22 portfolio, 7 pooling, 7 confirm-on-fail, 5 escaped-mass,
6 seq-gmm.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
- parity_check_hlm.py: set P.deltaF for ChooseFDModes models (required),
  add IMRPhenomXPHM / IMRPhenomXPNR to the default model list.
- uv_parity_diagnostics.py: add the X models to the default list.
- README: two X-family caveats. (1) Frame convention: raw ChooseFDModes
  modes satisfy the parity identity only after a global z-rotation by
  exactly pi (degenerate with phi_ref; harmless in marginalized PE, but
  naive complex residuals report O(1) violations). After removing it,
  XPNR raw modes are parity-clean to 1e-5 (amplitudes 1e-14); XPHM has a
  small genuine 2e-3 (2,+-1) aligned-spin asymmetry. (2) RIFT interface
  artifact: raw IMRPhenomXHM is exactly equatorially symmetric, but
  RIFT hlmoft's ChooseFDModes->TD conditioning introduces ~1% (2,+-2)
  spurious asymmetry, dominating through-RIFT residuals for the X family
  (lnL-impact crossings rho~70 for both XPHM and XPNR, same origin).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…ning

SimInspiralChooseFDModes returns modes on an ascending two-sided grid
[-fNyq, +fNyq] with DC at the center index (odd length TDlen+1), but the
FD high-pass window in hlmoft assigned frequencies via evaluate_fvals,
which assumes RIFT's reversed packing and, at odd length, is offset by
deltaF/2 (f_assumed = -f_true + deltaF/2).  Because (l,m) and (l,-m)
modes occupy opposite signs of f, the 'symmetric' window was shifted by
one bin between the members of each pair -- exactly in the taper band
containing the |h(f)| maximum -- violating the conjugate-pair identity
h_{l,-m}(f) = (-1)^l conj(h_{lm}(-f)) for every ChooseFDModes-consumed
model (IMRPhenomXHM/XPHM/XPNR/XO4a/HM), e.g. 1.02% (2,±2) amplitude
asymmetry for an exactly-symmetric aligned-spin XHM input.

Build the window on the true frequency grid (exactly even in f), and
zero the surviving -fNyq bin when the resize truncates its +fNyq
partner, so all conditioning commutes with complex conjugation.

Verification (nonprec control, M=80, fmin=fref=20, deltaF=1/16):
- mode-level TD conjugate-pair residual: 1.02e-2 -> 2.5e-16 (XHM);
  XPNR and HM machine-exact; XPHM reduced to its intrinsic model
  asymmetry; precessing XPHM amplitude identity at 4.2e-16.
- U/V diagnostics (test/waveforms/uv_parity_diagnostics.py):
  D2 1.04e-9 -> 1.28e-17, D3 1.7e-9 -> 1.4e-16.  (The corrupted band
  lay mostly below fmin, so U/V leakage was ~1e-9, not the mode-level
  1e-2; the fix removes the interface floor entirely.)
- Conditioned waveform change: 1.1e-2 full-band, 99% of it below
  fmin in the discarded turn-on band; 1.6e-3 in-band.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
hlmoff_for_LISA resizes the ascending two-sided grid returned by
SimInspiralChooseFDModes ([-fNyq,+fNyq], odd length TDlen+1) down to
TDlen, which truncates the +fNyq bin while keeping its -fNyq partner at
index 0.  Any model support at Nyquist then breaks the conjugate-pair
identity h_{l,-m}(f) = (-1)^l conj(h_{lm}(-f)) at that unpaired bin.
Zero the surviving -fNyq bin after the resize so the truncation
commutes with the f -> -f reflection, mirroring the hlmoft fix in
02e668e.

Unlike the main path, this branch applies no evaluate_fvals-based
high-pass window, so the one-bin window-shift defect fixed there does
not arise here; this closes the remaining truncation asymmetry.

Verification (IMRPhenomXHM, m1=44, m2=36, s1z=0.5, s2z=-0.3,
deltaF=1/16, fmin=fref=20): paired-bin FD residual exactly 0 and TD
conjugate-pair residual <= 3.6e-16 before and after; the -fNyq bin is
identically zero at this config, so the fix is defensive and changes
nothing there.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…rrected weights

[P1] _pool_replica_rvs rewrote log_joint_s_prior / joint_s_prior to encode the
per-replica renormalization and the flat fairdraw blocks, but concatenated the
existing 'log_weights' cache unchanged.  The .dgrid and calibration-posterior
exporters PREFER that cache -- they fall back to
log_integrand + log_joint_prior - log_joint_s_prior only when it is absent -- and
mcsamplerPortfolio writes it, so it is present in exactly the portfolio-replica
case under discussion.  Those scientific outputs therefore used the ORIGINAL
weights while the evidence used the corrected ones: replica rebalancing ignored,
and fairdraw blocks double-weighted again in the very products this pooling
exists to keep consistent.

'log_weights' (and a linear 'weights' if present) are now recomputed from the
pooled canonical columns.  If the components are unavailable to rebuild them the
stale cache is DROPPED rather than trusted, so consumers fall through to whatever
columns remain instead of silently using numbers that no longer describe the
record.

Two tests: the pooled cache must equal the pooled components AND must differ from
the stale per-replica values (so the test fails on the old behaviour rather than
passing vacuously), and a fairdraw block's cache must come out constant.

All suites pass: 22 portfolio, 9 pooling, 7 confirm-on-fail, 5 escaped-mass,
6 seq-gmm.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…_coverage

Pool MC-error replicas for export; warn on L0 truncated support; general setup() kwarg memory
…pair-O4d

LISA/lalsimutils_compat: zero surviving -fNyq bin after FD-mode resize
@oshaughn
oshaughn merged commit 200dd60 into oshaughn:rift_O4d Aug 7, 2026
19 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants