#!/usr/bin/env python
"""
integrate_likelihood_extrinsic_jax
==================================

A JAX, automatic-differentiation-based driver for the RIFT ILE extrinsic
likelihood.  It mirrors the *structure* and the relevant CLI/output conventions
of ``integrate_likelihood_extrinsic_batchmode`` but replaces the Monte-Carlo
likelihood evaluation with the differentiable JAX likelihood in
``RIFT.likelihood.jax_ile`` (a reimplementation of the production
``DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop`` "fused" code branch,
validated to ~1e-13 against the numpy reference).

REUSED unchanged from production RIFT (deliberately not reinvented):
  * frame reading            -> lalsimutils.frame_data_to_non_herm_hoff
  * PSD handling             -> lalsimutils.get_psd_series_from_xmldoc / resample
  * waveform + inner products-> factored_likelihood.PrecomputeLikelihoodTerms
  * array packing / epoch    -> factored_likelihood.PackLikelihoodDataStructuresAsArrays

NEW (the point of this driver):
  * the extrinsic -> lnL map is pure-JAX and AD-compatible, so we can
    gradient-ascend the (narrow) likelihood peak, form a Fisher approximation,
    and run efficient (gradient-aware) importance sampling.
  * DISTANCE MARGINALIZATION is performed analytically inside the JAX
    likelihood when ``--distance-marginalization`` is passed (ILE-compatible
    flag, default OFF).  The bare factored likelihood diverges on amplitude
    slivers (template power -> 0, e.g. inclination -> pi); marginalizing
    distance against the volumetric prior regulates this and yields a smooth,
    bounded objective over the five angular/sky parameters.  Required by
    ``--mode nuts``.

Full ILE argument compatibility: every integrate_likelihood_extrinsic_batchmode
option is accepted so this can be dropped into a production command line.
Implemented options are used; the rest are silently ignored (and reported);
science-changing options we do not implement (calibration marginalization, ROM,
NR templates, supplementary-likelihood, zero-likelihood, maximize-only) cause a
hard failure rather than a misleading result.  --sim-xml / --sim-grid and
--n-events-to-analyze enable batch processing of multiple intrinsic templates.

Run modes (``--mode``):
  prior-mc    : brute-force importance sampling from the physical prior
                (robust, slow; the classic ILE estimator).
  laplace-is  : prior-seeded *adaptive* Gaussian importance sampling
                (efficient for narrow peaks).  [default]
  map         : gradient-ascend the angular peak + report Fisher (AD demo).

Output (matches ILE conventions):
  <output>_<event>_.dat :
     event_id m1 m2 s1x s1y s1z s2x s2y s2z  lnL sigma_lnL ntotal neff
     ``lnL`` = log marginal likelihood (evidence) over extrinsic parameters,
     ``sigma_lnL = sqrt(var)/Z`` (ILE's ``sqrt(var)/res``).  Masses in M_sun.
  <output>_<event>_samples.dat (--save-samples) : per-sample extrinsic params +
     loglikelihood with ILE-style column names.  Equal-weight FAIR DRAW rows
     (no weight column), so a collapsed integration whose weights admit no fair
     draw writes NO samples file and fails the event instead.
"""

from __future__ import print_function

import os
import sys
import json
from optparse import OptionParser, OptionGroup

import numpy as np

# The opt-in CuPy precompute and JAX likelihood share one GPU. Set this before
# importing JAX or any RIFT module that may initialize its allocator.
if os.environ.get('RIFT_GPU_PRECOMPUTE', '0') == '1':
    os.environ.setdefault('XLA_PYTHON_CLIENT_PREALLOCATE', 'false')

import jax
import jax.numpy as jnp
jax.config.update("jax_enable_x64", True)

# Configure the persistent cache before importing modules that construct ILE
# JITs.  RIFT selects a compatibility-keyed child directory, so a shared cache
# root is safe across heterogeneous GPU/JAX installations.
from RIFT.jax_cache import (argv_option, configure_persistent_cache,
                            import_bundle, runtime_compatibility)
_JAX_CACHE_DIR = configure_persistent_cache(jax, sys.argv[1:])
_JAX_CACHE_BUNDLE = (argv_option(sys.argv[1:], "--jax-cache-bundle")
                     or os.environ.get("RIFT_JAX_CACHE_BUNDLE"))
_JAX_CACHE_PROFILE = argv_option(sys.argv[1:], "--jax-cache-profile")
if _JAX_CACHE_BUNDLE:
    if _JAX_CACHE_DIR is None:
        raise RuntimeError("a cache bundle was requested but no writable JAX cache is available")
    import_bundle(_JAX_CACHE_BUNDLE, _JAX_CACHE_DIR.parent,
                  runtime_compatibility(jax), _JAX_CACHE_PROFILE,
                  destination=_JAX_CACHE_DIR)

import lal
import lalsimulation as lalsim

import RIFT.lalsimutils as lalsimutils
from RIFT.likelihood.jax_ile import (
    build_data_from_precompute, build_rotation_data_from_precompute,
    build_freqresponse_data_from_precompute,
    build_rotating_freqresponse_data_from_precompute)
from RIFT.likelihood.jax_ile.wrapper import bandlimited_storage_requirement
from RIFT.likelihood.jax_ile import anglemarg as _anglemarg
from RIFT.likelihood.jax_ile.samplers import angle_marg_eval_chunk as _angle_marg_eval_chunk
from RIFT.likelihood.jax_ile.samplers import regularize_cov as _regularize_cov
# ONE definition, imported rather than re-typed: a second copy of the
# "is this evidence trustworthy" rule is a copy that drifts.
from RIFT.likelihood.jax_ile.samplers import _finalize_evidence
from RIFT.likelihood.jax_ile.core import _GATHERERS as _JAX_GATHERERS, JAX_INTERP_DEFAULT
from RIFT.likelihood.jax_ile.anglemarg import (ANGLE_MARG_DEFAULT, ANGLE_MARG_LEGACY,
                                               ANGLE_MARG_CHOICES,
                                               ANGLE_MARG_CROSSOVER_AMPLITUDE)
from RIFT.likelihood.jax_ile.direct_marginalization_policy import (
    POLICY_CHOICES as DIRECT_MARG_POLICY_CHOICES,
    POLICY_DEFAULT as DIRECT_MARG_POLICY_DEFAULT,
    PolicyConfig as DirectMargPolicyConfig,
    BoundedMultipeakConfig, validate_bounded_multipeak_config,
    RESERVE_SCHEME_CHOICES as DIRECT_MARG_RESERVE_SCHEME_CHOICES,
    reserve_pair as _direct_marg_reserve_pair,
    summarize_policy_ledger as _summarize_policy_ledger,
    predict_reserve_pair, format_reserve_pair, RESERVE_SCHEME_EXECUTABLE)
_JAX_GATHERER_NAMES = tuple(_JAX_GATHERERS)
from RIFT.likelihood.jax_ile.wrapper import (
    JAXExtrinsicLikelihood, JAXDistanceMarginalizedLikelihood,
    BOUNDED_MULTIPEAK_LOG_ZERO,
)

MSUN = lal.MSUN_SI
PC = lal.PC_SI

# Angular parameter order for the distance-marginalized (default) path.
ANG_NAMES = ("ra", "dec", "psi", "incl", "phiref")
# Full extrinsic order for the fixed-distance path.
FULL_NAMES = ("ra", "dec", "psi", "incl", "phiref", "distMpc")


# ---------------------------------------------------------------------------
# Drop-in ILE argument compatibility
# ---------------------------------------------------------------------------
# Every option of integrate_likelihood_extrinsic_batchmode is accepted so this
# driver can be substituted into an existing production command line.  Options
# we implement are defined explicitly in build_parser(); the rest are registered
# here as accepted-but-ignored (with the correct arity so parsing succeeds), and
# a small set that would silently change the *science* if ignored is failed on.

# ILE's default cap for --fairdraw-extrinsic-output (integrate_likelihood_
# extrinsic_batchmode: --fairdraw-extrinsic-output-n-max default=5).
_FAIRDRAW_N_MAX_DEFAULT = 5
# Below this many effective samples the exported cloud is not a usable posterior
# sample however it is drawn.  ONE definition: fairdraw_indices warns at it, the
# tempering guard refuses at it, and --allow-degenerate-tempering's help quotes
# it, so the three cannot drift apart.
_USABLE_EXPORT_ESS = 200
# Default export-ESS budget for --auto-adapt-weight-exponent.  Named so the
# parser default and the "was this actually passed?" check are ONE value.
_TARGET_EXPORT_ESS_FRAC_DEFAULT = 0.9
# Modes whose sampler actually implements the anneal.  samplers.flowmc_sample_phimarg
# takes temper_adapt; samplers.flowmc_sample (plain `flowmc`) does not.
_ADAPT_ADAPT_MODES = frozenset((
    "flowmc-phimarg", "flowmc-phipsimarg", "flowmc-dpsimarg"))

# Modes whose sampler reports a TEMPERED state plus a genuine importance weight
# (post_weight = L^(1-inv_T)); only these honour --adapt-weight-exponent, and
# only these consume post_weight at export.
_TEMPERED_MODES = frozenset((
    "flowmc", "flowmc-phimarg", "flowmc-phipsimarg", "flowmc-dpsimarg"))
# Modes whose --save-samples export is a reweighted FAIR DRAW.  Only here do the
# fair-draw count options do anything; elsewhere the export is the sampler's own
# chain and the count flags are inert (and must be reported as ignored).
_FAIRDRAW_MODES = _TEMPERED_MODES | frozenset(("prior-mc", "laplace-is"))
# Modes whose blind full-sky prior draws go through eval_lnL, which stops the
# run on one uncertified bandlimited row.  The flowMC family and multistart-nuts
# pilot through the samplers' own draw and are not on this list.
_EVAL_LNL_STOP_MODES = frozenset(("prior-mc", "laplace-is", "map", "nuts"))
# Largest shift a wrong-sky draw makes to one detector's arrival time: the
# window must contain it or a blind row's arrival peak sits beyond the edge,
# which no certificate can converge on (measured 35/256 rows at 20 ms, 0/256 at
# 50 ms; DESIGN_jax_bandlimited_distmarg.md, "The endpoint certificate").
_BANDLIMITED_FULLSKY_HALF_WINDOW_MIN = 2.0 * lal.REARTH_SI / lal.C_SI
# Modes whose posterior is already phi_ref-marginalised ANALYTICALLY (a grid
# sum baked into the likelihood, not a sampled axis): --phase-marginalization
# has nothing left to do on these and is reported IGNORED rather than refused
# (RO'S 2026-09-08: missing knobs are fine to no-op for compatibility, but
# must warn).
_PHASE_ANALYTIC_MODES = frozenset((
    "flowmc-phimarg", "flowmc-phipsimarg", "flowmc-dpsimarg", "nuts-phimarg"))

# Boolean (zero-argument) ILE options (action=store_true/false).
_ILE_BOOL_OPTS = {
    "--calibration-conjugate-phase", "--calibration-global-norm",
    "--check-good-enough", "--zero-likelihood", "--random-event",
    "--soft-fail-event-range", "--fmin-template-correct-for-lmax",
    "--internal-use-gwpy", "--nr-lookup", "--nr-hybrid-use", "--rom-use-basis",
    "--rom-integrate-intrinsic", "--nr-perturbative-extraction",
    "--nr-perturbative-extraction-full", "--nr-use-provided-strain",
    "--no-memory", "--use-gwsignal", "--use-external-EOB", "--maximize-only",
    "--dump-lnL-time-series", "--time-marginalization",
    "--resample-time-marginalization", "--distance-marginalization",
    "--calibration-fused-kernel", "--calibration-export-posterior",
    "--extrinsic-proposal-adapt", "--vectorized", "--gpu", "--force-gpu-only",
    "--force-xpy", "--save-samples", "--save-samples-process-params",
    "--internal-hard-fail-on-error", "--internal-soft-fail-on-cuda-error",
    "--internal-make-empty-file-on-error", "--internal-waveform-fd-no-condition",
    "--verbose", "--save-eccentricity", "--fairdraw-extrinsic-output",
    "--convergence-tests-on", "--no-adapt", "--force-adapt-all",
    "--force-reset-all", "--no-adapt-distance", "--no-adapt-after-first",
    "--adapt-adapt", "--adapt-log", "--d-prior-redshift",
    "--declination-cosine-sampler", "--inclination-cosine-sampler",
    "--internal-rotate-phase", "--internal-sky-network-coordinates",
    "--internal-sky-network-coordinates-raw", "--auto-logarithm-offset",
    "--pin-distance-to-sim", "--export-eos-index",
    "--export-marginal-distance-grid", "--adapt-intrinsic",
    "--reject-collapsed-live-volume", "--rotation-slow", "--freqresponse",
    "--internal-waveform-fd-L-frame", "--save-EOB-parameters",
    "--save-hyperbolic", "--force-hyperbolic-22", "--save-meanPerAno",
    "--internal-gmm-correlate-all", "--internal-gmm-adaptive-components",
    "--portfolio-varaha-never-freeze", "--portfolio-varaha-can-freeze",
    "--portfolio-adaptive-alloc", "--sampler-sequential-warmstart",
    "--sampler-l0-rescue-accept-truncated", "--sampler-anisotropic-bins",
    "--internal-reparam-dl-incl", "--internal-use-lnL",
    "--distance-slice-all-fresh",
    "--distance-slice-randomize",
    "--psi-marginalization",
}
# ILE options taking repeated values (action=append).
_ILE_APPEND_OPTS = {
    "--channel-name", "--psd-file", "--fmin-ifo", "--nr-lookup-group",
    "--parameter", "--parameter-range", "--sampler-portfolio",
    "--sampler-portfolio-args",
}
# The full ILE option set (so anything not implemented is still accepted).
_ILE_ALL_OPTS = {
    "--adapt-adapt", "--adapt-floor-level", "--adapt-intrinsic", "--adapt-log",
    "--adapt-weight-exponent", "--amp-order", "--approximant",
    "--auto-logarithm-offset", "--cache-file", "--calibration-burn-in-neff",
    "--calibration-burn-in-nmax", "--calibration-dump-responsibilities",
    "--calibration-conjugate-phase", "--calibration-global-norm",
    "--calibration-mc-error-extrinsic", "--calibration-n-realizations-max",
    "--calibration-neff-cal-target",
    "--calibration-envelope-directory", "--calibration-export-posterior",
    "--calibration-fused-kernel", "--calibration-n-realizations",
    "--calibration-pilot-extrinsic", "--calibration-proposal-breadcrumb",
    "--calibration-spline-count", "--channel-name", "--check-good-enough",
    "--coinc-xml", "--convergence-tests-on", "--data-end-time",
    "--data-integration-window-half", "--data-start-time",
    "--declination-cosine-sampler", "--deff-lambda", "--distance-marginalization",
    "--distance-marginalization-lookup-table", "--distance-slice-all-fresh",
    "--distance-slice-chunk", "--distance-slice-randomize",
    "--distance-slice-skip-threshold", "--distance-slice-wing-delta-lnL",
    "--distance-slice-wing-neff", "--distance-slice-wing-nmax", "--d-max",
    "--d-min", "--d-prior", "--d-prior-redshift", "--eff-lambda", "--e-freq",
    "--dump-lnL-time-series", "--event", "--event-time",
    "--export-distance-slices", "--export-eos-index",
    "--export-marginal-distance-grid", "--extrinsic-proposal-adapt",
    "--extrinsic-proposal-breadcrumb", "--extrinsic-proposal-output",
    "--extrinsic-proposal-field", "--extrinsic-proposal-field-cover-frac",
    "--extrinsic-proposal-field-inflate",
    "--fairdraw-extrinsic-output", "--fairdraw-extrinsic-output-n-max", "--fmax",
    "--fmin-ifo", "--fmin-template", "--fmin-template-correct-for-lmax",
    "--force-adapt-all", "--force-gpu-only", "--force-hyperbolic-22",
    "--force-reset-all", "--force-xpy", "--freqresponse",
    "--freqresponse-arm-length", "--freqresponse-qmax", "--gpu",
    "--check-slowrot-pmax", "--check-finite-size-Qmax",
    "--choose-slowrot-pmax", "--choose-slowrot-Qmax",
    "--response-order-snr", "--response-order-lnL-tol",
    "--response-order-sky-samples", "--response-order-p-reference",
    "--response-order-Q-reference", "--response-order-max-bank-gib",
    "--inclination-cosine-sampler", "--internal-gmm-adaptive-components",
    "--internal-gmm-correlate-all", "--internal-gmm-defensive-frac",
    "--internal-gmm-inflate", "--internal-gmm-max-components",
    "--internal-gmm-phase-components", "--internal-gmm-sky-components",
    "--internal-data-storage-window-half", "--internal-hard-fail-on-error",
    "--internal-make-empty-file-on-error", "--internal-precompute-ignore-threshold",
    "--internal-rotate-phase", "--internal-sky-network-coordinates",
    "--internal-sky-network-coordinates-raw", "--internal-soft-fail-on-cuda-error",
    "--internal-reparam-dl-incl", "--internal-use-gwpy", "--internal-use-lnL",
    "--internal-waveform-extra-kwargs",
    "--internal-waveform-extra-lalsuite-args", "--internal-waveform-fd-no-condition",
    "--internal-waveform-fd-L-frame", "--internal-waveform-taper",
    "--interpolate-time", "--inv-spec-trunc-time", "--l-max",
    "--limit-declination", "--limit-distance", "--limit-inclination", "--limit-psi",
    "--limit-right-ascension", "--manual-logarithm-offset", "--mass1", "--mass2",
    "--mc-error-ess-trigger", "--mc-error-khat-trigger", "--mc-error-replicas",
    "--mc-error-sigma-trigger",
    "--maximize-only", "--n-chunk", "--n-distance-slice-core", "--n-eff",
    "--n-events-to-analyze", "--n-fairdraw-extrinsic-samples", "--n-max",
    "--n-distance-slice-wing", "--no-adapt", "--no-adapt-after-first",
    "--no-adapt-distance", "--no-memory", "--nf-flow-load", "--nf-flow-save",
    "--nr-group", "--nr-hybrid-method", "--nr-hybrid-use", "--nr-index",
    "--nr-lookup", "--nr-lookup-group", "--nr-param", "--nr-params",
    "--nr-perturbative-extraction", "--nr-perturbative-extraction-full",
    "--nr-use-provided-strain", "--output-file", "--output-format",
    "--parameter", "--parameter-range", "--pin-distance-to-sim",
    "--right-ascension", "--declination", "--psi", "--distance",
    "--phi-orb", "--t-ref", "--inclination",
    "--portfolio-adaptive-alloc", "--portfolio-alloc-exponent",
    "--portfolio-freeze-wt", "--portfolio-grace-iters",
    "--portfolio-probe-period", "--portfolio-quality-signal",
    "--portfolio-revive-period", "--portfolio-varaha-can-freeze",
    "--portfolio-varaha-max-frac", "--portfolio-varaha-min-frac",
    "--portfolio-varaha-never-freeze", "--portfolio-weight-clip", "--psd-file",
    "--psd-window-shape", "--psi-marginalization", "--random-event", "--reference-freq",
    "--resample-time-marginalization", "--restricted-mode-list-file",
    "--rom-group", "--rom-integrate-intrinsic", "--rom-limit-basis-size-to",
    "--reject-collapsed-live-volume", "--rom-param", "--rom-use-basis",
    "--rotation-n-harmonics", "--rotation-p-max", "--rotation-slow",
    "--sampler-anisotropic-bins", "--sampler-l0-rescue-accept-truncated",
    "--sampler-l0-rescue-puff-factor", "--sampler-l0-rescue-puff-scale",
    "--sampler-l0-rescue-puff-width-frac", "--sampler-l0-rescue-reject-dlnZ",
    "--sampler-load-state", "--sampler-method", "--sampler-portfolio",
    "--sampler-portfolio-args", "--sampler-save-state",
    "--sampler-sequential-warmstart", "--sampler-sequential-warmstart-cover-frac",
    "--sampler-sequential-warmstart-deltalnL", "--sampler-warmstart-cover-frac",
    "--sampler-warmstart-inflate", "--sampler-warmstart-retry-neff",
    "--sampler-warmstart-samples", "--sampler-xpy", "--save-EOB-parameters",
    "--save-P", "--save-deltalnL", "--save-eccentricity", "--save-hyperbolic",
    "--save-meanPerAno",
    "--save-samples", "--save-samples-process-params", "--seed", "--sim-grid",
    "--sim-xml", "--skymap-file", "--soft-fail-event-range", "--spin1z",
    "--spin2z", "--srate", "--srate-internal",
    "--srate-resample-time-marginalization",
    "--supplementary-likelihood-factor-code",
    "--supplementary-likelihood-factor-function",
    "--supplementary-likelihood-factor-ini", "--time-marginalization",
    "--time-posterior-export", "--use-external-EOB", "--use-gwsignal",
    "--use-gwsignal-lmax-nyquist", "--vectorized", "--verbose",
    "--window-shape", "--zero-likelihood",
}


def _register_compat_options(optp):
    """Add every ILE option we do not already define, as accepted-but-ignored."""
    import optparse
    known = set()
    for opt in optp._get_all_options():
        known.update(opt._long_opts)
    g = OptionGroup(optp, "ILE compatibility (accepted; ignored unless implemented)")
    for name in sorted(_ILE_ALL_OPTS):
        if name in known:
            continue
        if name in _ILE_BOOL_OPTS:
            g.add_option(name, action="store_true", default=False,
                         help=optparse.SUPPRESS_HELP)
        elif name in _ILE_APPEND_OPTS:
            g.add_option(name, action="append", default=[],
                         help=optparse.SUPPRESS_HELP)
        else:
            g.add_option(name, default=None, help=optparse.SUPPRESS_HELP)
    optp.add_option_group(g)


def _dest(name):
    return name.lstrip("-").replace("-", "_")


def check_critical_and_report(opts, optp):
    """Fail on science-changing options we do not implement; report the rest."""
    def is_set(name):
        v = getattr(opts, _dest(name), None)
        return bool(v) if name in _ILE_BOOL_OPTS else (v is not None and v != [])

    fatal = []
    # --distance-gh-nodes / JAX_ILE_DISTMARG_GH resolution (RO'S 2026-09-08,
    # revised same day per adversarial review: the option must default to
    # None, not 0, or an explicit ``--distance-gh-nodes 0`` is indistinguish-
    # able from "not passed" and a nonzero JAX_ILE_DISTMARG_GH silently wins
    # -- the exact opposite of the documented "CLI wins").  The per-sample
    # Gauss-Hermite distance quadrature must be reachable by an ILE argument,
    # not only the environment variable core.py reads at import.  Resolution:
    # the CLI flag, when GIVEN (including an explicit 0), always wins over
    # the environment variable; the environment variable is consulted only
    # when the CLI flag was not given at all.  The two are REFUSED, not
    # silently reconciled, whenever BOTH are given and differ -- including an
    # explicit CLI 0 against a nonzero env, which is a conflict like any
    # other: silently keeping either value would hide the disagreement.
    # Only the RESOLUTION happens here (needed below by the
    # distance-grid-scheme combination check); the actual global mutation and
    # banner are deferred until after the fatal gate below, so a refused
    # command line never mutates core's module-level GH state.
    _gh_cli = getattr(opts, "distance_gh_nodes", None)
    _gh_env_raw = os.environ.get("JAX_ILE_DISTMARG_GH")
    _gh_env = int(_gh_env_raw) if _gh_env_raw not in (None, "") else None
    if _gh_cli is not None and _gh_cli < 0:
        fatal.append("--distance-gh-nodes must be >= 0, got %d" % _gh_cli)
    elif _gh_cli is not None and _gh_env is not None and _gh_cli != _gh_env:
        fatal.append(
            "--distance-gh-nodes %d conflicts with JAX_ILE_DISTMARG_GH=%d in "
            "the environment; unset the environment variable or pass the "
            "same value on the command line" % (_gh_cli, _gh_env))
    if _gh_cli is not None:
        _gh_resolved = _gh_cli
    elif _gh_env is not None:
        _gh_resolved = _gh_env
    else:
        _gh_resolved = 0
    if getattr(opts, "calibration_n_realizations", None) not in (None, 1) \
            or is_set("--calibration-export-posterior") \
            or is_set("--calibration-envelope-directory"):
        fatal.append("calibration marginalization (--calibration-*) is not implemented")
    if is_set("--rom-use-basis") or is_set("--rom-group") or is_set("--rom-param"):
        fatal.append("ROM-basis waveforms (--rom-*) are not implemented")
    if any(is_set(x) for x in ("--supplementary-likelihood-factor-code",
                               "--supplementary-likelihood-factor-function",
                               "--supplementary-likelihood-factor-ini")):
        fatal.append("supplementary likelihood factors are not implemented")
    if is_set("--nr-lookup") or is_set("--nr-params") or is_set("--nr-index"):
        fatal.append("NR-waveform templates (--nr-*) are not implemented")
    if is_set("--zero-likelihood"):
        fatal.append("--zero-likelihood is not implemented")
    if is_set("--maximize-only"):
        fatal.append("--maximize-only is not implemented (this driver integrates)")
    _is_n = int(getattr(opts, "smc_is_samples", 60000))
    if _is_n < 0:
        fatal.append("--smc-is-samples must be >= 0; a negative count reaches "
                     "the SMC proposal draw, whose exception handler would "
                     "silently fall back to the raw SMC evidence")
    _sampler_method = getattr(opts, "sampler_method", None)
    _jax_av_active = _sampler_method in ("AV", "portfolio")
    if _jax_av_active:
        _dp = str(getattr(opts, "d_prior", None) or "euclidean").strip().lower()
        if _dp not in ("euclidean", "volumetric", "pseudo_cosmo"):
            fatal.append("--sampler-method %s supports distance priors "
                         "Euclidean/volumetric and pseudo_cosmo, not %s"
                         % (_sampler_method, _dp))
        try:
            resolve_av_angular_limits(opts)
        except SystemExit as exc:
            fatal.append(str(exc).strip())
        if getattr(opts, "mode", None) in ("map", "nuts", "multistart-nuts",
                                            "nuts-phimarg"):
            fatal.append("--sampler-method %s cannot override the chain/optimizer "
                         "mode --mode %s; use laplace-is or a flowmc-* mode to "
                         "select the desired JAX likelihood geometry"
                         % (_sampler_method, getattr(opts, "mode", None)))
        if int(getattr(opts, "n_eff", None) or 1000) <= 0:
            fatal.append("--n-eff must be positive for --sampler-method %s"
                         % _sampler_method)
        if int(getattr(opts, "jax_av_seed_pilot", 4000)) <= 0:
            fatal.append("--jax-av-seed-pilot must be positive")
        if int(getattr(opts, "jax_av_seed_modes", 4)) <= 0:
            fatal.append("--jax-av-seed-modes must be positive")
        _seed_points = getattr(opts, "jax_av_seed_points", None)
        if _seed_points is not None and int(_seed_points) <= 0:
            fatal.append("--jax-av-seed-points must be positive")
        _eval_chunk = getattr(opts, "jax_av_eval_chunk", None)
        if _eval_chunk is not None and int(_eval_chunk) <= 0:
            fatal.append("--jax-av-eval-chunk must be positive")
        if float(getattr(opts, "jax_av_sky_inflate", 2.0)) <= 0:
            fatal.append("--jax-av-sky-inflate must be positive")
        _prior_frac = float(getattr(opts, "jax_av_seed_prior_frac", 0.1))
        if not (0.0 <= _prior_frac <= 1.0):
            fatal.append("--jax-av-seed-prior-frac must lie in [0,1]")
        _members = [s.strip().upper()
                    for item in (getattr(opts, "sampler_portfolio", None) or [])
                    for s in str(item).split(",") if s.strip()]
        if _sampler_method == "AV" and _members:
            fatal.append("--sampler-portfolio is inert with --sampler-method AV")
        if any(name not in ("AV", "GMM") for name in _members):
            fatal.append("the JAX portfolio supports only AV and GMM members; got %s"
                         % ",".join(_members))
        if getattr(opts, "jax_av_seed", "none") == "none":
            for _name in ("--jax-av-seed-pilot", "--jax-av-seed-modes",
                          "--jax-av-seed-points", "--jax-av-sky-inflate",
                          "--jax-av-seed-prior-frac"):
                if was_supplied(opts, _name):
                    fatal.append("%s is inert with --jax-av-seed none" % _name)
        for _name in ("--auto-adapt-weight-exponent", "--target-export-ess-frac",
                      "--allow-degenerate-tempering", "--adapt-adapt",
                      "--temper-init", "--temper-ess-frac",
                      "--temper-max-stages", "--temper-max-dbeta",
                      "--fisher-precondition", "--fisher-is-samples",
                      "--smc-puffball", "--smc-walkers", "--smc-move-steps",
                      "--smc-is-samples", "--smc-puff-scale"):
            if was_supplied(opts, _name):
                fatal.append("%s configures the mode's gradient/SMC sampler and "
                             "is inert under --sampler-method %s"
                             % (_name, _sampler_method))
    else:
        for _name in ("--sampler-portfolio", "--sampler-anisotropic-bins",
                      "--n-eff"):
            if was_supplied(opts, _name):
                fatal.append("%s is inert unless --sampler-method AV/portfolio"
                             % _name)
        for _name in ("--jax-av-seed", "--jax-av-seed-pilot",
                      "--jax-av-seed-modes", "--jax-av-seed-points",
                      "--jax-av-sky-inflate", "--jax-av-seed-prior-frac",
                      "--jax-av-eval-chunk"):
            if was_supplied(opts, _name):
                fatal.append("%s requires --sampler-method AV/portfolio" % _name)
        # ``was_supplied`` deliberately fails open when an in-process caller
        # has no argv provenance.  The non-default value itself is conclusive.
        if getattr(opts, "jax_av_seed", "none") != "none" \
                and not any("--jax-av-seed requires" in item for item in fatal):
            fatal.append("--jax-av-seed requires --sampler-method AV/portfolio")
    if is_set("--psi-marginalization"):
        fatal.append("--psi-marginalization (analytic polarization-angle "
                     "marginalization) is not implemented for psi-sampling modes "
                     "on this driver; use --mode flowmc-phipsimarg or "
                     "flowmc-dpsimarg, which marginalize psi by construction")
    # Static multipeak options must be effective and valid before precompute.
    _bounded = getattr(opts, "angle_marg_scheme", None) == "multipeak-jax"
    if _bounded:
        if getattr(opts, "mode", None) != "flowmc-phipsimarg":
            fatal.append("multipeak-jax requires --mode flowmc-phipsimarg")
        try:
            bounded_multipeak_config_from_options(opts)
        except (TypeError, ValueError) as exc:
            fatal.append("multipeak-jax: %s" % exc)
    else:
        for field in (*BoundedMultipeakConfig._fields, "decline_action"):
            flag = "--multipeak-jax-" + field.replace("_", "-")
            if was_supplied(opts, flag):
                fatal.append("%s requires --angle-marg-scheme multipeak-jax" % flag)
    # Cross-axis policy scope, at parse time (external review of #278, P1):
    # only --mode flowmc-phipsimarg reads the policy, so a request anywhere
    # else would otherwise complete on the ordinary likelihood, and the
    # policy's knobs are inert unless the policy is on.
    _policy = getattr(opts, "direct_marginalization_policy",
                      DIRECT_MARG_POLICY_DEFAULT)
    _policy_knobs = ("--direct-marginalization-time-guard",
                     "--direct-marginalization-reserve-time-refine",
                     "--direct-marginalization-reserve-time-refine-max",
                     "--direct-marginalization-error-budget-nats",
                     "--direct-marginalization-batch-rows",
                     "--direct-marginalization-policy-probe-rows",
                     "--direct-marginalization-policy-probe-only",
                     "--direct-marginalization-max-modes",
                     "--direct-marginalization-enriched-max-modes",
                     "--direct-marginalization-base-oversample",
                     "--direct-marginalization-enriched-oversample",
                     "--direct-marginalization-max-starts",
                     "--direct-marginalization-max-time-nodes",
                     "--direct-marginalization-convergence-tol-nats",
                     "--direct-marginalization-time-guard-tol-nats")
    if _policy != "off":
        if getattr(opts, "mode", None) != "flowmc-phipsimarg":
            fatal.append("--direct-marginalization-policy %s applies only to "
                         "--mode flowmc-phipsimarg; --mode %s would run the "
                         "ordinary likelihood and silently ignore the request"
                         % (_policy, getattr(opts, "mode", None)))
        _g = int(getattr(opts, "direct_marginalization_time_guard", 16))
        if _g < 2:
            fatal.append("--direct-marginalization-time-guard must be >= 2 "
                         "(the two-guard comparison needs a half guard)")
        _f = int(getattr(opts, "direct_marginalization_reserve_time_refine", 4))
        if _f < 2 or _f % 2:
            fatal.append("--direct-marginalization-reserve-time-refine must be "
                         "an even integer >= 2 (the check rule is the "
                         "half-refined rule)")
        _fm = int(getattr(opts, "direct_marginalization_reserve_time_refine_max",
                          32))
        if _fm < _f or _fm % 2:
            fatal.append("--direct-marginalization-reserve-time-refine-max must "
                         "be an even integer >= the refine factor")
        _b = float(getattr(opts, "direct_marginalization_error_budget_nats",
                           1.0e-3))
        if not (np.isfinite(_b) and _b > 0.0):
            fatal.append("--direct-marginalization-error-budget-nats must be "
                         "finite and positive")
        # --d-prior is accepted by the parser but not forwarded to the JAX
        # wrapper (which is always volumetric); refuse a non-volumetric request
        # under the policy here rather than let the wrapper's refusal be
        # unreachable (wiring review, item 6).
        _dp = getattr(opts, "d_prior", None)
        if _dp not in (None, "", "euclidean", "volumetric"):
            fatal.append("--direct-marginalization-policy %s derives its "
                         "measure from the volumetric distance prior; "
                         "--d-prior %s is not supported by the composite"
                         % (_policy, _dp))
        _br = int(getattr(opts, "direct_marginalization_batch_rows", 1))
        if _br < 0:
            fatal.append("--direct-marginalization-batch-rows must be >= 0 "
                         "(1 = row at a time, 0 = one full batch)")
        _pr = int(getattr(opts, "direct_marginalization_policy_probe_rows", 0))
        if _pr < 0:
            fatal.append("--direct-marginalization-policy-probe-rows must be "
                         ">= 0 (0 disables the pre-sampling probe)")
        if getattr(opts, "direct_marginalization_policy_probe_only", False) \
                and _pr <= 0:
            fatal.append("--direct-marginalization-policy-probe-only needs "
                         "--direct-marginalization-policy-probe-rows > 0; "
                         "otherwise it would exit having measured nothing")
        for _k, _d in (("--direct-marginalization-max-modes",
                        "direct_marginalization_max_modes"),
                       ("--direct-marginalization-enriched-max-modes",
                        "direct_marginalization_enriched_max_modes"),
                       ("--direct-marginalization-base-oversample",
                        "direct_marginalization_base_oversample"),
                       ("--direct-marginalization-enriched-oversample",
                        "direct_marginalization_enriched_oversample"),
                       ("--direct-marginalization-max-starts",
                        "direct_marginalization_max_starts")):
            if int(getattr(opts, _d, 1)) < 1:
                fatal.append("%s must be >= 1" % _k)
        for _k, _d in (("--direct-marginalization-convergence-tol-nats",
                        "direct_marginalization_convergence_tol_nats"),
                       ("--direct-marginalization-time-guard-tol-nats",
                        "direct_marginalization_time_guard_tol_nats")):
            _v = float(getattr(opts, _d, 1.0e-3))
            if not (np.isfinite(_v) and _v > 0.0):
                fatal.append("%s must be finite and positive" % _k)
        if (int(getattr(opts, "direct_marginalization_enriched_max_modes", 8))
                < int(getattr(opts, "direct_marginalization_max_modes", 4))):
            fatal.append("--direct-marginalization-enriched-max-modes must be "
                         ">= --direct-marginalization-max-modes: the enriched "
                         "plan has to be able to nest the base plan, and a "
                         "narrower one declines on mode nesting every row")
    else:
        for _k in _policy_knobs:
            if was_supplied(opts, _k):
                fatal.append("%s is inert without --direct-marginalization-"
                             "policy auto; pass the policy or drop the option"
                             % _k)
    # Distance-grid option combinations that the wrapper would reject anyway --
    # caught HERE, at parse time, so the user is not made to sit through a full
    # precompute first (F8 of external review).
    dgs = getattr(opts, "distance_grid_scheme", "uniform")
    if dgs != "uniform":
        if getattr(opts, "mode", None) != "flowmc-phipsimarg":
            fatal.append("--distance-grid-scheme %s applies only to --mode "
                         "flowmc-phipsimarg (it is validated only for the dense "
                         "angle-marginalization schemes)" % dgs)
        elif getattr(opts, "angle_marg_scheme", ANGLE_MARG_DEFAULT) == "grid":
            fatal.append(
                "--distance-grid-scheme %s requires --angle-marg-scheme "
                "exact/laplace/auto: the log-uniform grid is sized from the "
                "data-derived angle amplitude, which the 'grid' scheme you "
                "asked for does not compute" % dgs)
        if getattr(opts, "distance_grid_points", None) is not None:
            fatal.append("--distance-grid-points and --distance-grid-scheme %s "
                         "both set the distance node count; pass one or the "
                         "other" % dgs)
        _tol = getattr(opts, "distance_grid_tol", None)
        if _tol is not None and not (0.0 < float(_tol) < 2.0):
            # The option's VALUE, checked beside its combinations.  The range
            # is closed form -- c(tol) = pi*sqrt(2/ln(2/tol)) needs
            # 0 < tol < 2 -- so making the user sit through a full precompute
            # to be told otherwise is the same avoidable cost F8 removed for
            # the combinations.
            fatal.append("--distance-grid-tol must be in (0, 2): it is a "
                         "FRACTIONAL error on the distance integral (~nats on "
                         "lnL), not a node count; got %r.  Note the interval "
                         "is the DERIVATION's domain, not a range of sensible "
                         "values: c(tol) diverges as tol -> 2, so 1.999 asks "
                         "for a 2-node grid.  The shipped default is %g."
                         % (_tol, _jax_core_dist_tol_default()))
        if _gh_resolved > 0:
            fatal.append(
                "--distance-grid-scheme %s cannot be combined with "
                "--distance-gh-nodes/JAX_ILE_DISTMARG_GH: the per-sample "
                "Gauss-Hermite distance quadrature uses only the SUPPORT of "
                "the grid, so the option would be bit-identically inert "
                "while still being reported as active" % dgs)
    elif getattr(opts, "distance_grid_tol", None) is not None:
        fatal.append("--distance-grid-tol applies only to "
                     "--distance-grid-scheme loguniform; it would be silently "
                     "inert here")
    # NOT ``... or 1``: 0 is falsy, so that idiom would silently promote an invalid
    # ``--q-time-pregrid-factor 0`` to the default and report nothing.  (It did, until
    # test_jax_dropin_manifest... caught it.)
    _qf = getattr(opts, "q_time_pregrid_factor", 1)
    _qf = 1 if _qf is None else int(_qf)
    if _qf < 1:
        fatal.append("--q-time-pregrid-factor must be >= 1, got %d" % _qf)
    elif _qf != 1:
        if (getattr(opts, "rotation_slow", False)
                or getattr(opts, "freqresponse", False)):
            fatal.append("--q-time-pregrid-factor is not implemented for banded "
                         "rotation/frequency-response likelihoods")
        # NOT a blanket refusal any more (it was, before the JAX Q path grew a
        # pregrid).  What survives is the narrower, still-true refusal, and it
        # mirrors conventional ILE's (#261): the pregrid buys sub-sample accuracy
        # with a FOUR-tap cubic, so a stencil that cannot use a sub-sample position
        # ('nearest') or a different stencil chosen behind the user's back would
        # both make the two arms answer differently for the same command line --
        # which is exactly how they came to ship opposite stencil defaults (#233).
        #
        # 'nearest' is separately UNIMPLEMENTED, not merely suboptimal: see
        # core._q_sample_positions -- a refined-grid nearest gather no longer reads
        # the sample the banded post-phase reconstructs from the coarse index, so
        # the data term and the model norm would drift apart by up to half a coarse
        # bin.
        # The ``!= JAX_INTERP_DEFAULT`` clause is not redundant with was_supplied:
        # was_supplied() FAILS OPEN by design ("no record -> assume not supplied"),
        # which is the safe direction for the conflict checks it was written for and
        # the WRONG one here, where the consequence of guessing "not supplied" is
        # silently replacing a stencil the caller chose.  A caller with no supplied
        # record whose interp is not the module default has demonstrably chosen it,
        # so treat that as explicit too.
        _explicit_interp = (was_supplied(opts, "--interp")
                            or was_supplied(opts, "--interpolate-time")
                            or getattr(opts, "interp", None) != JAX_INTERP_DEFAULT)
        if _explicit_interp and getattr(opts, "interp", None) != "cubic":
            fatal.append(
                "--q-time-pregrid-factor %d uses four-tap cubic interpolation "
                "(--interp %r was requested); remove the explicit stencil option "
                "or set it to cubic" % (_qf, getattr(opts, "interp", None)))
        else:
            opts._q_pregrid_fallback_interp = getattr(opts, "interp", None)
            opts.interp = "cubic"
            print(" Q_lm pregrid: ENABLED factor=%d boundary=even-reflection "
                  "arrival_stencil=cubic integration_grid=unchanged" % _qf)
    _tq = getattr(opts, "time_marginalization_quadrature", "simpson")
    if getattr(opts, "rotation_slow", False) and _tq == "bandlimited":
        fatal.append("bandlimited time quadrature is not valid when slow rotation "
                     "makes <h|h> arrival-time dependent; use simpson")
    if ((getattr(opts, "rotation_slow", False) or getattr(opts, "freqresponse", False))
            and getattr(opts, "phase_marginalization", False)):
        fatal.append("phase marginalization is not implemented for banded "
                     "rotation/frequency-response likelihoods")
    if int(getattr(opts, "rotation_p_max", 0) or 0) < 0:
        fatal.append("--rotation-p-max must be >= 0")
    if int(getattr(opts, "rotation_n_harmonics", 2) or 2) < 2:
        fatal.append("--rotation-n-harmonics must be >= 2")
    if int(getattr(opts, "freqresponse_qmax", 4) or 4) < 0:
        fatal.append("--freqresponse-qmax must be >= 0")
    _rotation = bool(getattr(opts, "rotation_slow", False))
    _freqresponse = bool(getattr(opts, "freqresponse", False))
    _order_p = bool(getattr(opts, "check_slowrot_pmax", False)
                    or getattr(opts, "choose_slowrot_pmax", False))
    _order_q = bool(getattr(opts, "check_finite_size_Qmax", False)
                    or getattr(opts, "choose_slowrot_Qmax", False))
    if _order_p and not _rotation:
        fatal.append("pmax check/choice requires --rotation-slow")
    if _order_q and not _freqresponse:
        fatal.append("Qmax check/choice requires --freqresponse")
    if (_order_p or _order_q) and getattr(opts, "response_order_snr", None) is None:
        fatal.append("response-order check/choice requires --response-order-snr")
    if (getattr(opts, "check_slowrot_pmax", False)
            and getattr(opts, "choose_slowrot_pmax", False)):
        fatal.append("choose either --check-slowrot-pmax or --choose-slowrot-pmax")
    if (getattr(opts, "check_finite_size_Qmax", False)
            and getattr(opts, "choose_slowrot_Qmax", False)):
        fatal.append("choose either --check-finite-size-Qmax or --choose-slowrot-Qmax")
    if ((_order_p and _order_q)
            and ((getattr(opts, "check_slowrot_pmax", False)
                  or getattr(opts, "check_finite_size_Qmax", False))
                 and (getattr(opts, "choose_slowrot_pmax", False)
                      or getattr(opts, "choose_slowrot_Qmax", False)))):
        fatal.append("do not mix check and choose controls in one compound response scan")
    if float(getattr(opts, "response_order_lnL_tol", 0.1)) <= 0:
        fatal.append("--response-order-lnL-tol must be positive")
    if int(getattr(opts, "response_order_sky_samples", 128)) < 8:
        fatal.append("--response-order-sky-samples must be at least 8")
    if (int(getattr(opts, "response_order_p_reference", 2)) < 0
            or int(getattr(opts, "response_order_Q_reference", 8)) < 0):
        fatal.append("response-order reference orders must be nonnegative")
    if (not np.isfinite(float(getattr(opts, "response_order_max_bank_gib", 4.0)))
            or float(getattr(opts, "response_order_max_bank_gib", 4.0)) <= 0):
        fatal.append("--response-order-max-bank-gib must be finite and positive")
    if not _rotation:
        for _name in ("--rotation-n-harmonics", "--rotation-p-max"):
            if was_supplied(opts, _name):
                fatal.append("%s is inert without --rotation-slow" % _name)
    elif _freqresponse and was_supplied(opts, "--rotation-n-harmonics"):
        fatal.append("--rotation-n-harmonics does not set the compound bank width; "
                     "its exact support is fixed by Qmax and pmax")
    if not _freqresponse:
        for _name in ("--freqresponse-qmax", "--freqresponse-arm-length"):
            if was_supplied(opts, _name):
                fatal.append("%s is inert without --freqresponse" % _name)
    elif getattr(opts, "freqresponse_arm_length", None) not in (None, ""):
        try:
            _parse_freqresponse_arm_length(opts.freqresponse_arm_length)
        except (TypeError, ValueError) as exc:
            fatal.append("invalid --freqresponse-arm-length: %s" % exc)
    # Refused HERE, before precompute: with a full-sky prior the stop below is
    # not a chance event but the first chunk (see _EVAL_LNL_STOP_MODES).
    _hw = getattr(opts, "data_integration_window_half", None)
    if (_tq == "bandlimited" and (getattr(opts, "mode", None) in _EVAL_LNL_STOP_MODES
                                  or _jax_av_active)
            and _hw is not None
            and float(_hw) < _BANDLIMITED_FULLSKY_HALF_WINDOW_MIN):
        optp.error(
            "--time-marginalization-quadrature bandlimited with --mode %s "
            "evaluates full-sky prior draws and stops on one uncertified row, and "
            "--data-integration-window-half %.4g s is below 2 R_earth / c = %.4f s, "
            "the largest arrival shift a wrong-sky draw makes at one detector: a "
            "row whose arrival peak lies beyond the window edge cannot be "
            "certified, and the first chunk contains such rows (measured 35 of 256 "
            "blind rows at 0.02 s, 0 of 256 at 0.05 s).  Raise "
            "--data-integration-window-half to at least %.4f s plus the "
            "event-time uncertainty (0.05 s is the measured clean value), or use "
            "--time-marginalization-quadrature simpson."
            % (opts.mode, float(_hw), _BANDLIMITED_FULLSKY_HALF_WINDOW_MIN,
               _BANDLIMITED_FULLSKY_HALF_WINDOW_MIN))
    if fatal:
        optp.error("Cannot run as a faithful drop-in: " + "; ".join(fatal)
                   + ".  (These would silently change the result if ignored.)")

    # Nothing above raised: safe to APPLY the resolved distance-GH-nodes count
    # (deferred from the resolution comment above) and report it, once, here.
    from RIFT.likelihood.jax_ile.core import set_distmarg_gh_nodes
    set_distmarg_gh_nodes(_gh_resolved)
    opts._distance_gh_nodes_resolved = _gh_resolved
    print("  distance quadrature: per-sample Gauss-Hermite nodes=%d%s"
          % (_gh_resolved, " (legacy uniform grid)" if _gh_resolved == 0 else
             " (source: %s)" % ("--distance-gh-nodes" if _gh_cli is not None
                                else "JAX_ILE_DISTMARG_GH")))

    # Report accepted-but-ignored options that were actually passed.
    ignored = []
    implemented = {"--cache-file", "--channel-name", "--psd-file",
                   "--data-start-time", "--data-end-time", "--window-shape",
                   "--psd-window-shape", "--mass1", "--mass2", "--spin1z",
                   "--spin2z", "--eff-lambda", "--deff-lambda", "--approximant",
                   "--l-max", "--event-time", "--fmin-template",
                   "--reference-freq", "--fmax", "--srate",
                   "--srate-internal",
                   "--data-integration-window-half",
                   "--internal-data-storage-window-half", "--d-min", "--d-max",
                   "--limit-distance",
                   "--n-max", "--n-chunk", "--output-file",
                   "--event", "--save-samples", "--verbose", "--seed",
                   "--sim-xml", "--sim-grid", "--n-events-to-analyze",
                   "--random-event", "--distance-marginalization",
                   "--time-marginalization", "--time-marginalization-quadrature",
                   "--interpolate-time",
                   "--vectorized", "--use-gwsignal", "--rotation-slow",
                   "--internal-waveform-fd-L-frame",
                   "--internal-waveform-fd-no-condition",
                   "--internal-precompute-ignore-threshold", "--no-memory",
                   "--e-freq", "--fmin-ifo",
                   "--rotation-n-harmonics", "--rotation-p-max",
                   "--freqresponse", "--freqresponse-qmax",
                   "--freqresponse-arm-length", "--check-slowrot-pmax",
                   "--check-finite-size-Qmax", "--choose-slowrot-pmax",
                   "--choose-slowrot-Qmax", "--response-order-snr",
                   "--response-order-lnL-tol", "--response-order-sky-samples",
                   "--response-order-p-reference", "--response-order-Q-reference",
                   "--response-order-max-bank-gib"}
    # These are implemented PER MODE.  Listing them unconditionally would claim
    # they act under --mode laplace-is (the default), nuts, map, multistart-nuts
    # and nuts-phimarg, where they are inert -- exactly the silent no-op this
    # driver's compat layer exists to prevent.
    mode = getattr(opts, "mode", None)
    if getattr(opts, "sampler_method", None) in ("AV", "portfolio"):
        implemented |= {"--sampler-method", "--n-eff",
                        "--d-prior",
                        "--sampler-anisotropic-bins",
                        "--limit-right-ascension", "--limit-declination",
                        "--limit-psi", "--limit-inclination"}
        if getattr(opts, "sampler_method", None) == "portfolio":
            implemented.add("--sampler-portfolio")
    _tempered_mode_active = (mode in _TEMPERED_MODES
                             and getattr(opts, "sampler_method", None)
                             not in ("AV", "portfolio"))
    if _tempered_mode_active:
        # static tempering exponent (samplers.flowmc_sample*: inv_T = 1/temper)
        implemented.add("--adapt-weight-exponent")
    if mode in _FAIRDRAW_MODES:
        implemented |= {"--fairdraw-extrinsic-output",
                        "--fairdraw-extrinsic-output-n-max",
                        "--n-fairdraw-extrinsic-samples"}
    # --d-prior is handled separately below with a substantive message (the
    # JAX driver's distance prior, not merely "not yet implemented"), so it is
    # excluded from the generic bag here rather than lumped into it.
    for name in sorted(_ILE_ALL_OPTS - implemented - {"--d-prior"}):
        if is_set(name):
            ignored.append(name)
    # JAX-NATIVE tempering flags.  These are not in _ILE_ALL_OPTS (they have no ILE
    # counterpart), so the loop above cannot see them -- and they act ONLY on the
    # tempered modes.  Report them explicitly rather than let a user's chooser
    # request evaporate under --mode laplace-is, which is the driver default.
    if not _tempered_mode_active:
        inert = [n for n in ("--auto-adapt-weight-exponent",
                             "--allow-degenerate-tempering")
                 if getattr(opts, _dest(n), False)]
        if _target_ess_was_given(opts):
            inert.append("--target-export-ess-frac")
        if inert:
            print("Note: %s only act on the tempered modes (%s); --mode %s ignores "
                  "them." % (" ".join(sorted(set(inert))),
                             " ".join(sorted(_TEMPERED_MODES)), mode))
    # AV/portfolio forwards both the volumetric and pseudo-cosmological priors.
    # Other samplers remain volumetric and must visibly report a non-volumetric
    # compatibility request as ignored.
    _dp = getattr(opts, "d_prior", None)
    if (getattr(opts, "sampler_method", None) not in ("AV", "portfolio")
            and _dp not in (None, "") and str(_dp).strip().lower() not in (
                "euclidean", "volumetric")):
        print("Note: --d-prior %r is accepted but IGNORED by the JAX driver: "
              "it always integrates distance against the volumetric d^2 "
              "prior on [--d-min, --d-max]." % (_dp,))
    # --phase-marginalization: forwarded on the 6-D/5-D likelihoods, but the
    # four phi_ref-marginalised modes already integrate phase out analytically
    # (a grid sum baked into the likelihood, not a sampled axis) -- outside
    # _ILE_ALL_OPTS, so the generic loop above cannot see it either.
    if getattr(opts, "phase_marginalization", False) and mode in _PHASE_ANALYTIC_MODES:
        print("Note: --phase-marginalization is accepted but IGNORED for "
              "--mode %s (phase is already marginalized analytically)." % mode)
    # --sky-coordinates network: wired only for --mode multistart-nuts (the
    # other modes' samplers assume the equatorial (ra, sin(dec)) parameterization
    # end to end -- reparameterizing them is more than plumbing, see the PR).
    _skyc = getattr(opts, "sky_coordinates", "equatorial")
    if _skyc != "equatorial" and mode != "multistart-nuts":
        print("Note: --sky-coordinates %s is accepted but IGNORED for --mode "
              "%s (only --mode multistart-nuts samples the sky in network-"
              "frame coordinates; other modes use the default equatorial "
              "parameterization)." % (_skyc, mode))
    if ignored:
        print("Note: the following ILE options are accepted but IGNORED by the "
              "JAX driver (not yet implemented; behavior may differ from ILE):")
        print("   " + " ".join(ignored))


# ---------------------------------------------------------------------------
# CLI
# ---------------------------------------------------------------------------
DISTANCE_GRID_POINTS_DEFAULT = 256


def _jax_core_dist_tol_default():
    """Default distance-grid tolerance, read from the module that owns it."""
    from RIFT.likelihood.jax_ile.core import DIST_GRID_TOL_DEFAULT
    return DIST_GRID_TOL_DEFAULT


def build_parser():
    optp = OptionParser(usage="%prog [options]", description=__doc__)

    g = OptionGroup(optp, "JAX compilation cache")
    g.add_option("--jax-cache-dir", default=None,
                 help="Persistent cache root. RIFT adds a JAX/JAXLIB/backend/"
                      "GPU compatibility namespace (default: $RIFT_JAX_CACHE_ROOT "
                      "or $XDG_CACHE_HOME/rift/jax).")
    g.add_option("--no-jax-persistent-cache", action="store_true", default=False,
                 help="Disable cross-process JAX compilation caching for this run.")
    g.add_option("--jax-cache-bundle", default=None,
                 help="Validate and import a warmed rift_jax_cache bundle before "
                      "constructing any ILE JIT (also $RIFT_JAX_CACHE_BUNDLE).")
    g.add_option("--jax-cache-profile", default=None,
                 help="Require --jax-cache-bundle to declare this warmup profile.")
    optp.add_option_group(g)

    g = OptionGroup(optp, "Data input (frame mode)")
    g.add_option("--cache-file", default=None)
    g.add_option("--channel-name", action="append", default=[],
                 help="instrument=channel. Repeatable.")
    g.add_option("--psd-file", action="append", default=[],
                 help="instrument=psd.xml.gz. Repeatable.")
    g.add_option("--data-start-time", type=float, default=None)
    g.add_option("--data-end-time", type=float, default=None)
    g.add_option("--window-shape", type=float, default=0.0)
    optp.add_option_group(g)

    g = OptionGroup(optp, "Data input (injection / self-test mode)")
    g.add_option("--inj-mode", action="store_true", default=False,
                 help="Synthesize zero-noise data from the template at the "
                      "truth sky location (no frames; for testing).")
    g.add_option("--inj-ra", type=float, default=1.2)
    g.add_option("--inj-dec", type=float, default=-0.4)
    g.add_option("--inj-psi", type=float, default=0.7)
    g.add_option("--inj-incl", type=float, default=0.9)
    g.add_option("--inj-phiref", type=float, default=2.1)
    g.add_option("--inj-distance", type=float, default=600.0)
    g.add_option("--inj-detectors", default="H1,L1,V1")
    g.add_option("--inj-deltaF", type=float, default=0.25)
    optp.add_option_group(g)

    g = OptionGroup(optp, "Intrinsic template parameters")
    g.add_option("--mass1", type=float, default=None, help="m1 (solar masses)")
    g.add_option("--mass2", type=float, default=None, help="m2 (solar masses)")
    g.add_option("--spin1z", type=float, default=0.0)
    g.add_option("--spin2z", type=float, default=0.0)
    g.add_option("--eff-lambda", type=float, default=None,
                 help="Tidal Lambda-tilde; converted to lambda1,lambda2.")
    g.add_option("--deff-lambda", type=float, default=None,
                 help="Tidal delta-Lambda-tilde.")
    g.add_option("--approximant", default="IMRPhenomD")
    g.add_option("--l-max", type=int, default=2)
    g.add_option("--sim-xml", default=None,
                 help="ligolw XML of intrinsic params (ChooseWaveformParams array). "
                      "Selects event(s) via --event / --n-events-to-analyze.")
    g.add_option("--sim-grid", default=None,
                 help="ASCII grid file (names=True header) of intrinsic params.")
    g.add_option("--event", type=int, default=0,
                 help="Index of the first event to analyze from --sim-xml/--sim-grid.")
    g.add_option("--n-events-to-analyze", type=int, default=1,
                 help="Number of consecutive events to analyze (batch loop).")
    g.add_option("--random-event", action="store_true", default=False,
                 help="Draw events at random from the table instead of in order.")
    optp.add_option_group(g)

    g = OptionGroup(optp, "Conditioning / windows / frequencies")
    g.add_option("--event-time", type=float, default=None,
                 help="GPS geocenter event time (fiducial epoch).")
    g.add_option("--fmin-template", type=float, default=30.0)
    g.add_option("--reference-freq", type=float, default=30.0)
    g.add_option("--fmax", type=float, default=1000.0)
    g.add_option("--srate", type=float, default=4096.0)
    g.add_option("--data-integration-window-half", type=float, default=75e-3)
    g.add_option("--internal-data-storage-window-half", type=float, default=0.15)
    optp.add_option_group(g)

    g = OptionGroup(optp, "Time-dependent detector response")
    g.add_option("--rotation-slow", action="store_true", default=False,
                 help="Use the sidereal slow-rotation likelihood. May be combined "
                      "with --freqresponse.")
    g.add_option("--rotation-n-harmonics", type=int, default=2,
                 help="Minimum sidereal harmonic half-width (default 2).")
    g.add_option("--rotation-p-max", type=int, default=0,
                 help="Maximum propagation-delay derivative order (default 0).")
    g.add_option("--freqresponse", action="store_true", default=False,
                 help="Use finite-arm frequency response. Combining this with "
                      "--rotation-slow constructs a compound banded likelihood.")
    g.add_option("--freqresponse-qmax", type=int, default=4,
                 help="Highest finite-arm projection power (default 4).")
    g.add_option("--freqresponse-arm-length", default=None,
                 help="Arm length in metres, globally or DET=value comma list.")
    g.add_option("--check-slowrot-pmax", action="store_true", default=False,
                 help="Estimate the required pmax and warn if --rotation-p-max is too small.")
    g.add_option("--check-finite-size-Qmax", "--check-finite-size-qmax",
                 dest="check_finite_size_Qmax", action="store_true", default=False,
                 help="Estimate the required Qmax and warn if --freqresponse-qmax is too small.")
    g.add_option("--choose-slowrot-pmax", action="store_true", default=False,
                 help="Choose the least pmax satisfying the response error budget.")
    g.add_option("--choose-slowrot-Qmax", "--choose-finite-size-Qmax",
                 dest="choose_slowrot_Qmax", action="store_true", default=False,
                 help="Choose the least finite-size Qmax satisfying the response error budget.")
    g.add_option("--response-order-snr", type=float, default=None,
                 help="Target network SNR for response-order checks/choice (required when active).")
    g.add_option("--response-order-lnL-tol", type=float, default=0.1,
                 help="Allowed worst scanned Asimov likelihood loss (default 0.1).")
    g.add_option("--response-order-sky-samples", type=int, default=128,
                 help="Deterministic full-prior angular design size (default 128).")
    g.add_option("--response-order-p-reference", type=int, default=2,
                 help="Highest p used by the diagnostic reference bank (default 2).")
    g.add_option("--response-order-Q-reference", "--response-order-q-reference",
                 dest="response_order_Q_reference", type=int, default=8,
                 help="Highest Q used by the diagnostic reference bank (default 8).")
    g.add_option("--response-order-max-bank-gib", type=float, default=4.0,
                 help="Refuse a diagnostic whose dense U,V planning estimate exceeds "
                      "this many GiB (default 4).")
    optp.add_option_group(g)

    g = OptionGroup(optp, "Extrinsic exploration / sampling")
    g.add_option("--mode", default="laplace-is",
                 choices=["prior-mc", "laplace-is", "map", "nuts",
                          "multistart-nuts", "flowmc", "flowmc-phimarg",
                          "flowmc-phipsimarg", "flowmc-dpsimarg", "nuts-phimarg"],
                 help="Extrinsic exploration: prior-mc | laplace-is | map | "
                      "nuts (single-chain) | multistart-nuts (mode-covering) | "
                      "flowmc (normalizing-flow, 5-D) | "
                      "flowmc-phimarg (phi_ref-marginalised flowMC, 4-D — removes "
                      "the psi/phi_ref degeneracy ridge) | "
                      "flowmc-dpsimarg (distance+psi-marginalised flowMC, 4-D "
                      "(ra,dec,phiref,incl) — psi-only scan, cheaper than phimarg) | "
                      "nuts-phimarg (Fisher-whitened multi-start NUTS on the 4-D "
                      "phi_ref-marginalised posterior — whitening keeps the NUTS "
                      "step size O(1) at any SNR; recommended for high SNR). "
                      "The last five need --distance-marginalization.")
    g.add_option("--num-warmup", type=int, default=500,
                 help="NUTS warmup iterations (--mode nuts).")
    g.add_option("--num-samples", type=int, default=2000,
                 help="NUTS posterior samples per chain (--mode nuts).")
    g.add_option("--num-chains", type=int, default=1)
    # DEFAULT IS OFF (no re-use).  Measured: across an 8-event batch at two seeds,
    # re-using the trained flow contracts the extrinsic posterior monotonically in
    # slot index -- psi to ~40% of its no-re-use width by slot 7, on BOTH seeds,
    # with slot 0 (where no re-use has happened yet) sitting at ~1.0 as a control.
    # It buys nothing to pay for that with: 1589 s mean wall with re-use vs 1567 s
    # without, a difference smaller than the seed-to-seed spread and of flipping
    # sign.  See analyses/jax_extrinsic_tempering/ in the paper repo.
    #
    # --no-flow-reuse is kept (now a no-op restating the default) so existing
    # command lines and scripts that pass it keep working; --flow-reuse is the way
    # back to the old behaviour.  Both write the same dest, so the LAST one on the
    # command line wins.
    g.add_option("--no-flow-reuse", action="store_true", dest="no_flow_reuse",
                 default=True,
                 help="Re-train the flow from scratch each event (the DEFAULT). "
                      "Kept for compatibility: it now restates the default.")
    g.add_option("--flow-reuse", action="store_false", dest="no_flow_reuse",
                 help="Bootstrap the trained flow across --n-events-to-analyze "
                      "instead of re-training per event.  NOT recommended for any "
                      "run whose extrinsic SAMPLES are used: it contracts the "
                      "posterior in later slots (measured ~40% of the no-re-use "
                      "psi width by slot 7 of 8) for no measurable wall-time "
                      "saving.  Its evidence is less affected than its samples.")
    # ILE-compatible semantics: store_true, default OFF.  When ON, distance is
    # marginalized analytically (5-D angular problem, well conditioned); when
    # OFF, distance is sampled explicitly (6-D).  Required for --mode nuts.
    g.add_option("--distance-marginalization", action="store_true", default=False,
                 help="Marginalize distance analytically (regulates the "
                      "factored-likelihood amplitude degeneracy; required for nuts).")
    g.add_option("--n-max", type=int, default=300000)
    g.add_option("--n-chunk", type=int, default=8000)
    g.add_option("--n-eff", type=int, default=None,
                 help="Target effective sample count for --sampler-method AV/portfolio "
                      "(default 1000 on those backends).")
    g.add_option("--sampler-method", default=None,
                 help="Optional non-AD integration backend: AV or portfolio.  The "
                      "existing --mode still selects the JAX likelihood geometry; "
                      "other legacy values remain accepted as compatibility no-ops.")
    g.add_option("--sampler-portfolio", action="append", default=[],
                 help="Portfolio members (repeat or comma-separate; JAX path supports "
                      "AV and GMM, default AV,GMM).")
    g.add_option("--sampler-anisotropic-bins", action="store_true", default=False,
                 help="Allocate AV bins preferentially along compressed coordinates.")
    g.add_option("--jax-av-seed", type="choice", choices=("none", "fisher-sky"),
                 default="none",
                 help="Optional JAX initializer for AV/portfolio: hill-climb several "
                      "modes, use their Fisher curvature on sky, and draw all other "
                      "coordinates from the physical prior (default none).")
    g.add_option("--jax-av-eval-chunk", type=int, default=None,
                 help="Fixed JAX likelihood batch inside each AV coverage chunk "
                      "(default min(--n-chunk, the JAX memory-aware default)).")
    g.add_option("--jax-av-seed-pilot", type=int, default=4000,
                 help="Prior likelihood evaluations used to find hill-climb starts.")
    g.add_option("--jax-av-seed-modes", type=int, default=4,
                 help="Maximum separated sky modes retained by the Fisher-sky seed.")
    g.add_option("--jax-av-seed-points", type=int, default=None,
                 help="Seed-cloud size (default --n-chunk).")
    g.add_option("--jax-av-sky-inflate", type=float, default=2.0,
                 help="Standard-deviation inflation of the local Fisher sky proposal.")
    g.add_option("--jax-av-seed-prior-frac", type=float, default=0.1,
                 help="Fraction of the seed cloud drawn from the full physical prior. "
                      "The portfolio GMM, not this finite cloud, guarantees full support.")
    g.add_option("--d-min", type=float, default=1.0, help="Min distance (Mpc).")
    g.add_option("--d-max", type=float, default=10000.0, help="Max distance (Mpc).")

    g.add_option("--distance-grid-points", type=int, default=None,
                 help="Node count for the DEFAULT uniform-in-distance grid "
                      "(default %d).  --distance-grid-scheme loguniform "
                      "derives its own node count from the data, so passing "
                      "both is REFUSED rather than silently ignoring this one."
                      % DISTANCE_GRID_POINTS_DEFAULT)
    g.add_option("--distance-grid-scheme", type="choice",
                 choices=("uniform", "loguniform"), default="uniform",
                 help="Distance quadrature for --mode flowmc-phipsimarg with "
                      "--angle-marg-scheme exact/laplace/auto.  'uniform' "
                      "(DEFAULT) is the historical fixed --distance-grid-points "
                      "grid over the whole prior range, so existing command "
                      "lines reproduce existing runs.  'loguniform' places the "
                      "SAME kind of static grid uniformly in ln d, with a node "
                      "count DERIVED from the run's own data-derived angle "
                      "amplitude so that every per-sample distance peak "
                      "(relative width 1/rho) is resolved anywhere in "
                      "[--d-min, --d-max].  See --distance-grid-tol and "
                      "RIFT/likelihood/jax_ile/DESIGN_jax_distance_quadrature.md.")
    g.add_option("--distance-grid-tol", type=float, default=None,
                 help="Target FRACTIONAL error of the distance quadrature "
                      "(~nats on lnL) for --distance-grid-scheme loguniform; "
                      "the node count follows from it in closed form.  "
                      "Default %g.  Only valid with that scheme."
                      % _jax_core_dist_tol_default())
    g.add_option("--distance-gh-nodes", type=int, default=None,
                 help="Node count for the per-sample Gauss-Hermite-style "
                      "distance quadrature (RIFT.likelihood.jax_ile.core."
                      "make_distance_gh / _distmarg_gh_logL): nodes centred "
                      "PER SAMPLE on the Gaussian peak of the distance "
                      "integrand, resolving it to machine precision at any "
                      "SNR with a few dozen nodes, instead of the legacy "
                      "fixed grid (see --distance-grid-points) which "
                      "under-resolves that peak at high SNR.  Not passing "
                      "this flag keeps the legacy uniform grid (current "
                      "behaviour, no change).  Equivalent to the environment "
                      "variable JAX_ILE_DISTMARG_GH, which is still honoured "
                      "for compatibility, but ONLY when this flag is not "
                      "given: when the flag IS given -- including an "
                      "explicit 0 -- it wins over the environment variable, "
                      "and a DIFFERENT nonzero JAX_ILE_DISTMARG_GH is "
                      "REFUSED rather than silently picked (an explicit "
                      "--distance-gh-nodes 0 against a nonzero "
                      "JAX_ILE_DISTMARG_GH is a conflict too, and is refused "
                      "the same way: silently choosing 0 or the environment "
                      "value would both hide the disagreement). "
                      "Not compatible with --distance-grid-scheme loguniform "
                      "(the per-sample quadrature reads only the SUPPORT of "
                      "that grid, so combining them would be silently inert "
                      "while reporting as active -- see --distance-grid-tol).")
    g.add_option("--limit-distance", default=None,
                 help="Restrict distance SAMPLING (and, when distance is "
                      "marginalized, the distance QUADRATURE) to 'LO,HI' in Mpc, "
                      "WITHOUT changing the prior or its normalization: the prior "
                      "stays normalized over the full [--d-min,--d-max], so the "
                      "reported lnZ needs no correction and stays comparable to a "
                      "full-range run and to the batchmode ILE -- comparable to the "
                      "SAMPLER's own systematic, not to machine precision: measured "
                      "end to end on real data at rho ~ 82, the box truncates <= 0.003 "
                      "nats of evidence while the observed lnZ difference is +0.37 +- "
                      "0.11 nats the other way, being the full-range run's own bias "
                      "that the box removes (RIFT_roboto_paper "
                      "analyses/limit_distance_e2e/).  That benefit is a THRESHOLD "
                      "in amplitude, not a general improvement: the same "
                      "measurement at rho 41 gives +0.13 +- 0.08 nats, consistent "
                      "with zero, because the full-range sampler is already "
                      "healthy there.  Intended for high "
                      "amplitude, where the distance posterior narrows as 1/rho.  "
                      "Likelihood outside the box is simply not integrated, so keep "
                      "the box comfortably wider than the posterior.  Must lie "
                      "inside [--d-min,--d-max].")
    g.add_option("--phase-marginalization", action="store_true", default=False)
    g.add_option("--time-marginalization-quadrature", type="choice",
                 choices=("simpson", "bandlimited"), default="simpson",
                 help="Rule for the time integral: historical fixed-grid "
                      "Simpson or adaptive reflected-FFT refinement followed "
                      "by a converged trapezoid. No resolution knob is exposed; "
                      "the factor is derived and rechecked per row. "
                      "'bandlimited' is honoured by the 6-D modes and, under "
                      "--distance-marginalization, by every mode except the "
                      "phi/psi-marginalized ones, which refuse it and say so. "
                      "One uncertified row stops the run.  The modes that push "
                      "full-sky prior draws through that stop (prior-mc, "
                      "laplace-is, map, nuts) are refused at parse time unless "
                      "--data-integration-window-half is at least 2 R_earth / c "
                      "(0.0426 s), the largest arrival shift of a wrong-sky "
                      "draw; 0.05 s is the measured clean value.")
    g.add_option("--resample-time-marginalization", action="store_true",
                 default=False,
                 help="Conventional ILE option; currently unsupported by JAX ILE.")
    g.add_option("--srate-resample-time-marginalization", type="int", default=None,
                 help="Conventional ILE option; currently unsupported by JAX ILE.")
    g.add_option("--q-time-pregrid-factor", type="int", default=1,
                 help="OPT-IN. Refine the STORED rholm buffers onto a factor-x finer "
                      "time grid ONCE, before sampling, by even reflection plus a "
                      "band-limited FFT interpolation, and evaluate detector arrival "
                      "times on it with the four-tap cubic stencil.  The integration "
                      "cadence -- deltaT, tvals and the Simpson weights -- is "
                      "UNCHANGED: this refines how Q is INTERPOLATED, not what the "
                      "likelihood integrates over.  Default 1 is the historical "
                      "behaviour and is bit-identical.  8 is the validated value on "
                      "both arms (RIFT PR #261 for conventional ILE); larger factors "
                      "are accepted but buy little, because past ~8 the residual is "
                      "the reflection boundary condition rather than the stencil "
                      "step.  Selects --interp cubic and REFUSES an explicit "
                      "different stencil, exactly as conventional ILE does.  Costs "
                      "factor-x device memory for Q (only Q -- not the per-sample "
                      "gather scratch, which is the large term).  Measured accuracy "
                      "and cost: RIFT/likelihood/DESIGN_q_window_stencil.md 9.7.")
    g.add_option("--n-phi", type=int, default=32,
                 help="phi_ref grid size for --mode flowmc-phimarg (default 32; "
                      "use 64-128 for l-max>=4 or production quality).")
    g.add_option("--n-psi", type=int, default=8,
                 help="psi (polarization) grid size for --mode flowmc-phipsimarg "
                      "(default 8; spin-2 -> exponential convergence, ~8 exact for "
                      "l_max=2).  Keep nphi*npsi small (~64) -- the (phi,psi) scan "
                      "cost scales with it.  CAVEAT: 'exact' above holds for "
                      "SAMPLING the trig polynomial lnL, not for the quadrature "
                      "of exp(lnL), whose peak width is ~1/SNR -- the grid "
                      "scheme's marginalization error grows without bound with "
                      "SNR (measured ~1e2 nats at SNR 40 for npsi=8).  See "
                      "--angle-marg-scheme for the fix.")
    # No type=: optparse infers type="choice" from choices=, and raises
    # "must not supply choices for type 'string'" if both are given.  --interp
    # above omits it for the same reason.
    g.add_option("--angle-marg-scheme", default=ANGLE_MARG_DEFAULT,
                 choices=sorted(ANGLE_MARG_CHOICES),
                 help="(phi_ref, psi) marginalization scheme for --mode "
                      "flowmc-phipsimarg.  DEFAULT '%s' since 2026-09-02; it "
                      "was '%s' before, and passing '%s' reproduces a "
                      "pre-2026-09-02 run exactly.  '%s': the historical "
                      "--n-phi x --n-psi quadrature -- it averages exp(lnL), "
                      "whose peak width is ~1/SNR, on those nodes, so its "
                      "error grows WITHOUT BOUND with SNR; at rho 652 it "
                      "displaced a recovered sky position by 0.53 deg and "
                      "ranked that artifact ABOVE the injection.  Kept for "
                      "reproducing archived runs, not for new ones.  'exact': "
                      "Fourier-coefficient bootstrap + dense reconstruction; "
                      "expensive likelihood calls fixed by MODE CONTENT, never "
                      "by SNR (cost ~SNR^2).  'laplace': analytic Laplace in "
                      "psi + dense phi, error O(1/SNR^2), cheaper than 'exact' "
                      "above rho ~200 -- but it cannot use the per-sample "
                      "adaptive distance quadrature, so on the default uniform "
                      "distance grid it was measured 43 nats from 'exact' at "
                      "rho 163.  'auto': select exact/laplace from the run's "
                      "own data-derived amplitude.  The scheme that actually "
                      "ran is printed.  See RIFT.likelihood.jax_ile.anglemarg."
                      % (ANGLE_MARG_DEFAULT, ANGLE_MARG_LEGACY,
                         ANGLE_MARG_LEGACY, ANGLE_MARG_LEGACY))
    g.add_option("--multipeak-jax-decline-action", type="choice",
                 choices=("drop", "refuse"), default="drop",
                 help="Handle bounded multipeak declines (default drop): "
                      "drop gives proposals zero numerical weight and labels "
                      "accepted-region evidence; refuse aborts on declined "
                      "log_likelihood batches only. Scalar MAP/Fisher/MALA "
                      "declines return NaN without raising or latching. "
                      "Important declined mass requires changing the config; "
                      "neither action runs a reserve.")
    # Each cap is a construction-time scalar.  The complete configuration is
    # recorded on the likelihood, and no reserve is available at any setting.
    for field, default in BoundedMultipeakConfig()._asdict().items():
        g.add_option("--multipeak-jax-" + field.replace("_", "-"),
                     type="int" if isinstance(default, int) else "float",
                     default=default,
                     help=("Static multipeak-jax %s (default %s); requires "
                           "--angle-marg-scheme multipeak-jax."
                           % (field.replace("_", " "), default)))
    g.add_option("--direct-marginalization-policy",
                 default=DIRECT_MARG_POLICY_DEFAULT,
                 choices=sorted(DIRECT_MARG_POLICY_CHOICES),
                 help="OPT-IN cross-axis policy for --mode flowmc-phipsimarg "
                      "(default '%s').  'auto': per likelihood evaluation, "
                      "build the bounded U,V,Q start portfolios, attempt the "
                      "four-axis peak-local integral over (t, phi_ref, psi, D), "
                      "and keep it only if every acceptance diagnostic passes "
                      "(capacity, finite stationary modes, nested geometry, "
                      "base/enriched agreement, nested quadrature, two-guard "
                      "time, omitted-time mass, total error score under "
                      "--direct-marginalization-error-budget-nats); otherwise "
                      "run the band-limited exact-angle reserve.  No SNR "
                      "threshold is coded.  Requires --angle-marg-scheme exact "
                      "(its reserve), the simpson time rule (its check rule), "
                      "the volumetric distance prior and a uniform distance "
                      "grid; anything else is refused, not ignored.  Value-only: "
                      "gradient parity is not yet validated (see "
                      "DESIGN_direct_marginalization_policy.md), so keep this "
                      "off for production until that ladder is recorded.  The "
                      "branch statistics of the exported samples are printed and "
                      "labelled in the output headers.  A row the controller "
                      "cannot warrant (after doubling the reserve rule up to "
                      "16x) is nan and the run is NOT published.  Applies only "
                      "to --mode flowmc-phipsimarg; any other mode is refused."
                      % DIRECT_MARG_POLICY_DEFAULT)
    g.add_option("--direct-marginalization-time-guard", type=int,
                 default=DirectMargPolicyConfig().time_guard,
                 help="Primitive-only support samples added at each end of the "
                      "time window for the policy's guarded reconstruction "
                      f"(default {DirectMargPolicyConfig().time_guard}; must "
                      "be >= 2).  Both the local integral and "
                      "the reserve are re-evaluated at half this guard and must "
                      "agree within the error budget; a too-small guard declines "
                      "to the reserve rather than biasing.")
    g.add_option("--direct-marginalization-reserve-time-refine", type=int,
                 default=DirectMargPolicyConfig().reserve_time_refine,
                 help="Refinement factor of the reserve's band-limited time rule "
                      "over the native cadence (default "
                      f"{DirectMargPolicyConfig().reserve_time_refine}; even, "
                      ">= 2).  The "
                      "half-refined rule is the check rule the reserve must "
                      "agree with, so the warrant is a convergence statement "
                      "about the refined rules; the native Simpson rule's own "
                      "error is what the refinement removes.")
    g.add_option("--direct-marginalization-reserve-scheme", type="choice",
                 choices=list(DIRECT_MARG_RESERVE_SCHEME_CHOICES),
                 default="exact",
                 help="The policy's reserve, a PAIR (angular kernel, time rule): "
                      "exact = exact angles on the whole-window refined rule "
                      "(default, unchanged); laplace = psi-Laplace angles on "
                      "that rule; peaklocal = psi-Laplace angles on a fixed-"
                      "count time rule sized per row from the PREDICTED peak "
                      "width 1/(2 pi rho sigma_f) (rho from the row's table, "
                      "sigma_f the rms frequency of the stored Q), with fine "
                      "blocks on one commensurate lattice around the time "
                      "maxima located on the primitive itself plus a coarse "
                      "scan of the window; its node count and "
                      "cost do not grow with rho, and the ledger prints the "
                      "prediction behind every row.  auto = choose the pair "
                      "from the precomputed inputs before any row runs "
                      "(predict_reserve_pair).  The whole-window refine "
                      "options are inert under peaklocal; a failed warrant "
                      "there doubles the fine lattice at fixed span, at most "
                      "--direct-marginalization-peaklocal-escalations times, "
                      "and is reported.  Measured: "
                      "DESIGN_direct_marginalization_policy.md, "
                      "\"Peak-local time reserve\".")
    g.add_option("--direct-marginalization-peaklocal-escalations", type=int,
                 default=int(DirectMargPolicyConfig().reserve_peaklocal_escalations),
                 help="Under the peaklocal reserve schemes, how many times a "
                      "failed resolution warrant may double the fine lattice at "
                      "fixed span before the row is unusable (default: the "
                      "PolicyConfig value).  The "
                      "first tier is sized from the predicted width, so an "
                      "escalation is a finding about the prediction, not a loop; "
                      "the ledger counts them.")
    g.add_option("--direct-marginalization-reserve-time-refine-max", type=int,
                 default=DirectMargPolicyConfig().reserve_time_refine_max,
                 help="Ceiling of the reserve rule's escalation on a failed "
                      "warrant (default "
                      f"{DirectMargPolicyConfig().reserve_time_refine_max}; "
                      "even, >= the refine factor).  It "
                      "also bounds REVERSE-MODE GRADIENT MEMORY: the dense "
                      "reserve's backward pass costs ~2.6 GiB per unit of "
                      "refinement per evaluated row on a 614-sample window "
                      "(measured, rho 163), and lax.cond reserves memory for the "
                      "largest tier whether or not it runs; 32 asked for 85 GiB. "
                      "LOWERING IT IS NOT FREE.  It does not change any value "
                      "that gets computed, but it changes which rows get a "
                      "value at all: a declined row gets one reserve attempt "
                      "per tier, and a row still unwarranted at the ceiling "
                      "returns nan, which this driver refuses to publish.  "
                      "Measured on full-sky prior draws at rho 163, refine-max "
                      "4 left 6 of 8 rows nan and the run could not complete; "
                      "at the default 32 the same rows escalated and all were "
                      "usable.  So on a small card the run fails rather than "
                      "returning a slightly worse number.  Lower it only with "
                      "the nan count in the ledger in front of you.")
    g.add_option("--direct-marginalization-error-budget-nats", type=float,
                 default=DirectMargPolicyConfig().total_value_error_budget_nats,
                 help="Value-error allowance, in nats, for the RESERVE's "
                      "time-resolution warrant (default "
                      f"{DirectMargPolicyConfig().total_value_error_budget_nats:g}"
                      "; 1e-3 shipped).  "
                      "It does NOT reach the local accept/decline decision, "
                      "despite what this help used to say: that decision is the "
                      "four-term boolean at all_axis_peaklocal.py:862 and carries "
                      "no tolerance, so acceptance measured identical (48/64 at "
                      "rho 652) across 1e-3, 1e-2 and 1e-1.  What it buys is "
                      "reserve cost: at rho 41 on identical rows, 1e-3 -> 1e-2 "
                      "took escalations from 2 to 0 and the run from 321.1 s to "
                      "146.9 s, a 2.19x speedup, moving lnL by 4.8e-12 nats.  "
                      "Raised to 1e-2 by RO on 2026-09-08 (evening).  Measured at "
                      "rho 41; the speedup is not yet confirmed at high SNR.")
    g.add_option("--direct-marginalization-max-time-nodes", type=int,
                 default=DirectMargPolicyConfig().max_time_nodes,
                 help="Time-node capacity of the policy's LOCAL plan (default "
                      f"{DirectMargPolicyConfig().max_time_nodes}"
                      "; 64 shipped, and was reachable from no flag).  A row "
                      "whose live time-node "
                      "count exceeds this declines to the exact reserve, which "
                      "costs ~4.7 h at rho 652 against 2.3 s for an accepted "
                      "row, so the capacity is sized to avoid the fallback "
                      "rather than to bound the plan.  Measured need at "
                      "rho 652 on full-sky prior draws: median 72, p90 217, "
                      "max 614; caps of 64/128/256/1024 hold 44/73/91/100% of "
                      "rows.  The library default was 64 and was unreachable "
                      "from any flag.  NOT value-neutral: it changes which "
                      "rows take the local branch.")
    g.add_option("--direct-marginalization-batch-rows", type=int, default=1,
                 help="Rows the policy's controller executes together under "
                      "one vmap (default 1: row at a time; 0: one full batch). "
                      "VALUES are identical at every size -- lnL, every branch "
                      "decision and the gradients -- and a test pins that.  "
                      "COST IS A REGRESSION, not neutral: under vmap the "
                      "accept/reserve cond becomes a select, so a locally "
                      "ACCEPTED row also executes the dense reserve it would "
                      "otherwise skip, and every escalation tier runs for "
                      "every row.  The penalty therefore scales with the "
                      "locally accepted fraction.  Device workspace grows as "
                      "0.046 + 0.385 B GiB, so a 24 GiB card holds 32 rows and "
                      "not 64.  Leave it at 1 unless a measurement on YOUR "
                      "rows says otherwise; see "
                      "DESIGN_direct_marginalization_policy.md.  "
                      "WITHDRAWN: an earlier version of this text said COST "
                      "ONLY and quoted 96.3 s per row at 1 against 99.9 at 8. "
                      "Those rows all declined (accepted_local 0), so they "
                      "priced the decline path, and 37-50% of them were nan "
                      "under a reserve-refine ceiling of 4.  Do not cite them.")
    g.add_option("--direct-marginalization-policy-probe-rows", type=int,
                 default=0,
                 help="Evaluate the policy ledger on this many prior draws "
                      "BEFORE sampling and print the accept/decline counts "
                      "(default 0 = off).  The end-of-run note reports the same "
                      "counts on the exported rows, but only a run that reaches "
                      "the export can be classified that way, and a declining "
                      "run is exactly the one that does not get there: each "
                      "declined row pays three exact reserve evaluations and "
                      "rows execute one at a time under jax.lax.map.  The probe "
                      "answers 'local or reserve' in one compile plus this many "
                      "rows.  The compile is not extra work for a policy run -- "
                      "the end-of-run note jits the identical function.")
    g.add_option("--direct-marginalization-policy-probe-only",
                 action="store_true", default=False,
                 help="Exit after the probe above, before any sampling.  For "
                      "surveying operating points: it reports which acceptance "
                      "diagnostic is declining without paying for a posterior.  "
                      "Writes no samples and no evidence.")
    # Start-portfolio and mode-plan sizing.  Every default below is the value
    # PolicyConfig already carried, so an unchanged command line is unchanged.
    # They are exposed because the four-axis controller's own accepting
    # operating point (analyses/va_sequence_20260902/RESULTS_20260907_aap268_ladder.md,
    # 7 of 8 accepts at rho 326.15) reaches the SAME function through a richer
    # portfolio than PolicyConfig's defaults build, and the difference was
    # unreachable from the command line -- so the policy could only ever be
    # observed declining.
    g.add_option("--direct-marginalization-max-modes", type=int,
                 default=DirectMargPolicyConfig().max_modes,
                 help="Modes retained in the policy's BASE plan (default "
                      f"{DirectMargPolicyConfig().max_modes}).")
    g.add_option("--direct-marginalization-enriched-max-modes", type=int,
                 default=DirectMargPolicyConfig().enriched_max_modes,
                 help="Modes retained in the policy's ENRICHED plan (default "
                      f"{DirectMargPolicyConfig().enriched_max_modes}).  "
                      "Acceptance compares the two plans, so an enriched "
                      "plan that cannot resolve the extra structure declines.")
    g.add_option("--direct-marginalization-base-oversample", type=int,
                 default=DirectMargPolicyConfig().base_oversample,
                 help="Angular oversampling of the base start portfolio "
                      f"(default {DirectMargPolicyConfig().base_oversample}).")
    g.add_option("--direct-marginalization-enriched-oversample", type=int,
                 default=DirectMargPolicyConfig().enriched_oversample,
                 help="Angular oversampling of the enriched start portfolio "
                      f"(default {DirectMargPolicyConfig().enriched_oversample}).")
    g.add_option("--direct-marginalization-max-starts", type=int,
                 default=DirectMargPolicyConfig().base_max_starts,
                 help="Cap on ranked optimizer starts per plan (default "
                      f"{DirectMargPolicyConfig().base_max_starts}).  "
                      "TWO THINGS THIS IS NOT.  (1) It is not the whole of "
                      "decline_capacity.  That ledger key reports the "
                      "conjunction at all_axis_peaklocal.py:862 -- "
                      "norm_nonnegative, a certified time cover, the time-node "
                      "capacity, AND this cap -- so a row is charged to "
                      "'capacity' for four different reasons.  Raising this "
                      "ALONE recovers ~16% and then plateaus -- at rho 652 on "
                      "43 such rows, 36 failed the TIME-NODE capacity and only "
                      "7 were count-driven.  That plateau describes the PAIR, "
                      "not this knob: it held while the time-node capacity was "
                      "pinned at the library default of 64.  Size this with "
                      "--direct-marginalization-max-time-nodes, which now "
                      "exists.  Measured acceptance at rho 652 as (time nodes, "
                      "starts): (64, 32) 28%, (256, 32) 31%, (256, 128) 75%, "
                      "(512, 256) 77%.  (2) It is not a "
                      "cost knob: it selects which modes are retained, so it "
                      "moves the integral.  Measured on identical rows at rho "
                      "163, 32 -> 128 shifted already-accepted values by 4.8e-4 "
                      "and 3.5e-3 nats, at and above the 1e-3 error budget.  "
                      "An A/B on this may not treat its value difference as "
                      "noise.  32 -> 128 was approved by RO on 2026-09-08 "
                      "(evening) with that shift on the record, because a "
                      "declined row runs the exact reserve at hours per row at "
                      "rho 652 against ~2.3 s for an accepted one.  It is "
                      "sized WITH --direct-marginalization-max-time-nodes; "
                      "either alone plateaus near 16%.")
    g.add_option("--direct-marginalization-convergence-tol-nats", type=float,
                 default=DirectMargPolicyConfig().convergence_tol_nats,
                 help="Agreement required between the base and enriched plans, "
                      "and between each plan's own two quadrature orders, in "
                      "nats (default 1e-3).  This is the tolerance the "
                      "'enrichment' and 'quadrature' declines are measured "
                      "against; it is SEPARATE from "
                      "--direct-marginalization-error-budget-nats, which is the "
                      "total value allowance.")
    g.add_option("--direct-marginalization-time-guard-tol-nats", type=float,
                 default=DirectMargPolicyConfig().time_guard_tol_nats,
                 help="Agreement required between the full-guard and half-guard "
                      "evaluations, in nats (default 1e-3).  This is the "
                      "tolerance the two-guard time warrant is measured "
                      "against, for the local branch and the reserve alike.")
    # flowMC tuning (modes flowmc / flowmc-phimarg).  Defaults match
    # samplers.flowmc_sample*; exposed so pipeline Makefiles can tune them.
    g.add_option("--n-training-loops", type=int, default=4,
                 help="flowMC training loops (flow fit refinements).")
    g.add_option("--n-production-loops", type=int, default=4,
                 help="flowMC production loops (post-training sampling).")
    g.add_option("--n-epochs", type=int, default=10,
                 help="flowMC NF training epochs per training loop.")
    g.add_option("--n-local-steps", type=int, default=20,
                 help="flowMC local (MALA) steps per loop.")
    g.add_option("--n-global-steps", type=int, default=20,
                 help="flowMC global (flow-proposal) steps per loop.")
    g.add_option("--n-prior-pilot", type=int, default=8000,
                 help="Prior pilot draws used to initialize chains (and to "
                      "seed multi-start optimizers).")
    g.add_option("--adapt-weight-exponent", type=float, default=1.0,
                 help="Static likelihood-tempering exponent beta (a la old-RIFT): "
                      "flowMC samples the broadened target L^beta * pi and "
                      "reweights to the true posterior. beta=1 (default) is OFF; "
                      "beta<1 broadens (helps the flow find sharp high-SNR peaks). "
                      "Maps to sampler temper = 1/beta (modes flowmc, flowmc-phimarg).")
    g.add_option("--auto-adapt-weight-exponent", action="store_true", default=False,
                 help="EXPERIMENTAL, and not sufficiently validated: the ESS "
                      "model it inverts is calibrated at SNR ~= 23.8 only and is "
                      "optimistic at higher SNR, so the fraction it targets is "
                      "not guaranteed. No paper result uses it. "
                      "Choose --adapt-weight-exponent automatically from the "
                      "EXPORT budget instead of by hand: pick the smallest beta "
                      "whose reweighted export keeps --target-export-ess-frac of "
                      "the cloud.  Keyed on the SAMPLED DIMENSION, not on SNR -- "
                      "the non-JAX helper's SNR rule does not transfer to this "
                      "path (jax_ile/DESIGN_jax_tempering.md).  Tempered modes only.")
    g.add_option("--target-export-ess-frac", type=float,
                 default=_TARGET_EXPORT_ESS_FRAC_DEFAULT,
                 help="Fraction of the sampler cloud the reweighted --save-samples "
                      "export must retain, for --auto-adapt-weight-exponent "
                      "(default 0.9).  ESS/N = [beta(2-beta)]^(dim/2).")
    g.add_option("--allow-degenerate-tempering", action="store_true", default=False,
                 help="Silence the warning issued when the ESTIMATED export ESS "
                      "is below the %d-sample usability floor.  The estimate is "
                      "advisory (calibrated at one SNR), so this suppresses a "
                      "warning rather than overriding a refusal." % _USABLE_EXPORT_ESS)
    g.add_option("--adapt-adapt", action="store_true", default=False,
                 help="Adaptive likelihood tempering (flowmc-phimarg): anneal the "
                      "tempering exponent inv_T from --temper-init up to 1.0, "
                      "warm-starting the flow each stage, choosing each step so the "
                      "tempered-reweight ESS stays >= --temper-ess-frac. The "
                      "high-SNR path that avoids the static-tempering neff->1 "
                      "collapse. Overrides --adapt-weight-exponent when set.")
    g.add_option("--temper-init", type=float, default=0.02,
                 help="Initial (widest) tempering exponent inv_T for --adapt-adapt.")
    g.add_option("--temper-ess-frac", type=float, default=0.5,
                 help="Target tempered-reweight ESS fraction per anneal step (--adapt-adapt).")
    g.add_option("--temper-max-stages", type=int, default=16,
                 help="Max anneal stages before forcing the final inv_T=1 pass (--adapt-adapt).")
    g.add_option("--temper-max-dbeta", type=float, default=0.15,
                 help="Cap on the inv_T step per anneal stage, so the thermodynamic-"
                      "integration evidence trapezoid stays accurate (--adapt-adapt).")
    g.add_option("--fisher-precondition", action="store_true", default=False,
                 help="flowmc-phi(psi)marg: whiten the sample space around the MAP "
                      "using the observed Fisher (theta=MAP+F^{-1/2} y) so the "
                      "(1/SNR)-narrow high-SNR posterior is O(1) for the flow. "
                      "Restarts chains at the MAP. Recommended at SNR>=320 to avoid "
                      "the sky-sample collapse; harmless (~no-op) at low SNR.")
    g.add_option("--smc-puffball", action="store_true", default=False,
                 help="flowmc-phi(psi)marg: use adaptive SMC with a cloud-covariance "
                      "'puffball' random-walk move instead of the flow. Robust at any "
                      "SNR (sample->puffball->sample, like RIFT-AV/nested sampling); "
                      "does not collapse on sharp high-SNR peaks. Recommended at "
                      "SNR>=320 for reliable sky samples.")
    g.add_option("--smc-walkers", type=int, default=2000,
                 help="SMC puffball: number of walkers in the cloud (default 2000).")
    g.add_option("--smc-move-steps", type=int, default=10,
                 help="SMC puffball: random-walk Metropolis moves per temperature rung.")
    g.add_option("--smc-is-samples", type=int, default=60000,
                 help="Draws in the SMC cloud-Gaussian importance-sampling "
                      "evidence stage (default 60000, the sampler's own value). "
                      "This stage is a FIXED cost paid after the ladder "
                      "finishes, and on a likelihood whose rows execute one at "
                      "a time -- the direct-marginalization policy -- it "
                      "dominates the whole run.  Lowering it trades the IS "
                      "evidence for the raw SMC evidence: the ESS gate then "
                      "declines the IS value and the reported logZ falls back "
                      "to the SMC estimator, which the run already prints.")
    g.add_option("--smc-puff-scale", type=float, default=1.0,
                 help="SMC puffball: proposal scale x sqrt(cloud covariance) (~1.0).")
    g.add_option("--fisher-is-samples", type=int, default=0,
                 help="flowmc-phi(psi)marg high-SNR FALLBACK: after annealing, draw "
                      "this many sky samples by Fisher-whitened importance sampling "
                      "about the (Newton-polished) MAP and reweight by the true lnL "
                      "(no flow training -> immune to the SNR>=640 NF-collapse). "
                      "Overrides the exported samples; TI evidence is unchanged. "
                      "0=off; ~40000 recommended at SNR>=640. Implies --fisher-precondition.")
    # Choices come from the gatherer registry rather than a literal list: a stencil added to
    # _GATHERERS and not here would be unreachable from the command line, which is how 'cubic'
    # spent a release being implemented but unselectable.
    g.add_option("--interp", default=JAX_INTERP_DEFAULT, choices=sorted(_JAX_GATHERER_NAMES),
                 help="Interpolation of the precomputed rholm timeseries in arrival time. "
                      "DEFAULT CHANGED 2026-08-26 from 'linear' to '" + JAX_INTERP_DEFAULT + "': "
                      "linear undershoots the rholm peak by an amount that depends on where it "
                      "falls between samples, which biases the recovered arrival time and hence "
                      "the sky, and at high SNR -- the only regime this driver is used in -- it "
                      "is worse than 'nearest'. THIS CHANGES RESULTS for anyone who did not pass "
                      "--interp; pass '--interp linear' to reproduce a pre-2026-08-26 run. "
                      "'cubic' and 'sinc' mirror the production factored_likelihood stencils of "
                      "the same names; which of the two is more accurate depends on how "
                      "oversampled the rholm timeseries is (on fmin and srate as well as mass) "
                      "and there is no automatic rule -- see "
                      "RIFT/likelihood/DESIGN_q_window_stencil.md. COST: 'sinc' is ~1.33x "
                      "'cubic' in GPU runtime and needs ~12x its XLA scratch (~512 MB against "
                      "~40 MB at the default --n-chunk 8000; measured, see §9.5).")
    g.add_option("--sky-coordinates", default="equatorial",
                 choices=["equatorial", "network"],
                 help="Optional: 'network' samples the sky in the two-detector "
                      "baseline frame (folds the time-delay ring); multistart-nuts "
                      "only, falls back to equatorial if <2 detectors.")
    g.add_option("--proposal-inflate", type=float, default=3.0,
                 help="Scale on the adaptive Gaussian proposal covariance.")
    g.add_option("--seed", type=int, default=0)
    optp.add_option_group(g)

    g = OptionGroup(optp, "Output")
    g.add_option("--output-file", default=None)
    g.add_option("--save-samples", action="store_true", default=False)
    g.add_option("--fairdraw-extrinsic-output", action="store_true", default=False,
                 help="Bound the fair-draw export to "
                      "--fairdraw-extrinsic-output-n-max samples (as ILE does). "
                      "The export is ALWAYS a fair draw when the sampler "
                      "supplies importance weights; this only caps the count.")
    g.add_option("--fairdraw-extrinsic-output-n-max", type=int, default=None,
                 help="Cap on fair draws per evaluation when "
                      "--fairdraw-extrinsic-output is set (ILE default 5).  "
                      "Left as None when unset so the ignored-option report does "
                      "not claim the user passed it; resolved to 5 downstream.")
    g.add_option("--n-fairdraw-extrinsic-samples", type=int, default=None,
                 help="Export exactly this many fair draws, clamped only by the "
                      "rows available. Those already carry ILE's 1.5*ESS cap "
                      "wherever the export weights were non-uniform, applied "
                      "against the EXPORT weights; there is deliberately no "
                      "second clamp by the evidence estimator's neff, which "
                      "describes a different cloud. "
                      "Overrides --fairdraw-extrinsic-output-n-max.")
    g.add_option("--verbose", action="store_true", default=False)
    optp.add_option_group(g)

    # Register every remaining ILE option for drop-in compatibility (accepted,
    # silently ignored unless implemented; science-critical ones fail below).
    _register_compat_options(optp)
    return optp


# ---------------------------------------------------------------------------
# Data loading
# ---------------------------------------------------------------------------
def make_template(opts, fiducial_epoch, deltaF, deltaT):
    P = lalsimutils.ChooseWaveformParams()
    P.m1 = opts.mass1 * MSUN
    P.m2 = opts.mass2 * MSUN
    P.s1z = opts.spin1z
    P.s2z = opts.spin2z
    if opts.eff_lambda is not None:
        P.lambda1, P.lambda2 = lalsimutils.tidal_lambda_from_tilde(
            opts.mass1, opts.mass2, opts.eff_lambda, opts.deff_lambda or 0.0)
    P.fmin = opts.fmin_template
    P.fref = opts.reference_freq
    P.deltaT = deltaT
    P.deltaF = deltaF
    P.fmax = 0.0
    P.approx = lalsim.GetApproximantFromString(opts.approximant)
    P.radec = True
    P.tref = fiducial_epoch
    P.dist = 1000.0 * 1e6 * PC
    return P


def _analysis_delta_t(opts):
    """Waveform/precompute cadence, distinct from the input-data cadence."""
    rate = (opts.srate if getattr(opts, "srate_internal", None) is None
            else opts.srate_internal)
    if rate is None or float(rate) <= 0:
        raise ValueError("--srate/--srate-internal must be positive")
    return 1.0 / float(rate)


def load_templates(opts, fiducial_epoch, deltaF, deltaT):
    """Return a list of intrinsic templates (ChooseWaveformParams).

    Mirrors the production driver's --sim-xml / --sim-grid / --mass1,--mass2
    handling, including the [event : event+n_events] slice and the
    extrinsic-zeroing / fiducial-distance conventions.  Extrinsic parameters
    (sky, distance, inclination, phase, psi) are placeholders -- they are the
    quantities this code samples/marginalizes.
    """
    def _finalize(P):
        P.radec = True
        P.fref = opts.reference_freq
        P.fmin = opts.fmin_template
        P.deltaT = deltaT
        P.deltaF = deltaF
        P.fmax = 0.0
        P.tref = fiducial_epoch
        P.dist = 1000.0 * 1e6 * PC  # fiducial template distance (= distMpcRef)
        # Extrinsic angles are sampled later, not baked into the base modes.
        # Match conventional ILE: some mode generators encode these angles in
        # h_lm, so retaining XML/grid values would apply them twice.
        P.phiref = 0.0
        P.psi = 0.0
        P.incl = 0.0
        if opts.approximant:
            P.approx = lalsim.GetApproximantFromString(opts.approximant)
        return P

    if opts.sim_xml:
        P_all = lalsimutils.xml_to_ChooseWaveformParams_array(str(opts.sim_xml))
        if opts.random_event:
            idx = np.random.choice(len(P_all),
                                   size=min(opts.n_events_to_analyze, len(P_all)),
                                   replace=False)
            P_list = [P_all[i] for i in idx]
        else:
            if len(P_all) <= opts.event:
                print(" Event index out of range; soft exit"); sys.exit(0)
            hi = min(len(P_all), opts.event + opts.n_events_to_analyze)
            P_list = list(P_all[opts.event:hi])
        return [_finalize(P) for P in P_list]

    if opts.sim_grid:
        grid = np.genfromtxt(opts.sim_grid, names=True)
        grid = np.atleast_1d(grid)
        names = list(grid.dtype.names)
        params = [n for n in names if n in set(lalsimutils.valid_params)]
        if len(grid) <= opts.event:
            print(" Event index out of range; soft exit"); sys.exit(0)
        hi = min(len(grid), opts.event + opts.n_events_to_analyze)
        P_list = []
        for row in grid[opts.event:hi]:
            P = lalsimutils.ChooseWaveformParams()
            for name in params:
                val = float(row[name])
                if hasattr(P, name):
                    setattr(P, name, val)
                else:
                    P.assign_param(name, val)
            if P.m1 < 1e15:           # masses given in solar masses -> SI
                P.m1 *= MSUN; P.m2 *= MSUN
            P_list.append(_finalize(P))
        return P_list

    # single template from --mass1/--mass2
    if opts.mass1 is None or opts.mass2 is None:
        raise SystemExit("Provide --mass1/--mass2, or --sim-xml / --sim-grid.")
    return [make_template(opts, fiducial_epoch, deltaF, deltaT)]


def load_injection(opts, fiducial_epoch):
    detectors = opts.inj_detectors.split(",")
    deltaT = _analysis_delta_t(opts)
    P = make_template(opts, fiducial_epoch, opts.inj_deltaF, deltaT)
    P.phi, P.theta = opts.inj_ra, opts.inj_dec
    P.psi, P.incl, P.phiref = opts.inj_psi, opts.inj_incl, opts.inj_phiref
    P.dist = opts.inj_distance * 1e6 * PC
    data_dict, psd_dict = {}, {}
    for det in detectors:
        Pdet = P.copy(); Pdet.detector = det
        data_dict[det] = lalsimutils.non_herm_hoff(Pdet)
        psd_dict[det] = lalsim.SimNoisePSDaLIGOZeroDetHighPower
    P_template = make_template(opts, fiducial_epoch, opts.inj_deltaF, deltaT)
    return P_template, data_dict, psd_dict, detectors, True


def load_frames(opts, fiducial_epoch):
    deltaT = 1.0 / opts.srate
    deltaT_internal = (None if getattr(opts, "srate_internal", None) is None
                       else _analysis_delta_t(opts))
    data_dict, psd_dict = {}, {}
    fmin_ifo = {}
    for item in getattr(opts, "fmin_ifo", []) or []:
        inst, value = item.split("=", 1)
        fmin_ifo[inst] = float(value)
    for inst, chan in map(lambda c: c.split("="), opts.channel_name):
        if opts.verbose:
            print("Reading channel %s:%s from %s" % (inst, chan, opts.cache_file))
        data_dict[inst] = lalsimutils.frame_data_to_non_herm_hoff(
            opts.cache_file, inst + ":" + chan,
            start=opts.data_start_time, stop=opts.data_end_time,
            window_shape=opts.window_shape, deltaT=deltaT,
            deltaT_internal=deltaT_internal)
    for inst, psdf in map(lambda c: c.split("="), opts.psd_file):
        if opts.verbose:
            print("Reading PSD for %s from %s" % (inst, psdf))
        psd_dict[inst] = lalsimutils.get_psd_series_from_xmldoc(psdf, inst)
        psd_dict[inst] = lalsimutils.resample_psd_series(
            psd_dict[inst], data_dict[inst].deltaF)
        psd_window_shape = float(getattr(opts, "psd_window_shape", None) or 0.0)
        if psd_window_shape > 0 or opts.window_shape > 0:
            psd_factor = lalsimutils.psd_windowing_factor(
                psd_window_shape, len(psd_dict[inst].data.data))
            data_factor = lalsimutils.psd_windowing_factor(
                opts.window_shape, len(data_dict[inst].data.data))
            psd_dict[inst].data.data *= data_factor / psd_factor
        if inst in fmin_ifo:
            freqs = (psd_dict[inst].f0
                     + psd_dict[inst].deltaF
                     * np.arange(psd_dict[inst].data.length))
            psd_dict[inst].data.data[freqs < fmin_ifo[inst]] = 0
    detectors = list(data_dict.keys())
    return data_dict, psd_dict, detectors, False


# ---------------------------------------------------------------------------
# Priors (the angular block is shared; distance added only when NOT marginalized)
# ---------------------------------------------------------------------------
def resolve_distance_limit(opts):
    """(lo, hi) in Mpc that distance is actually SAMPLED / integrated over.

    Equals (--d-min, --d-max) unless --limit-distance narrowed it.  The prior
    normalization always stays on (--d-min, --d-max): see log_prior() and the
    d_prior_range= argument threaded into the marginalized likelihoods.
    """
    v = getattr(opts, "limit_distance", None)
    if not v:
        return float(opts.d_min), float(opts.d_max)
    from RIFT.integrators.mcsampler import distance_limit_range
    try:
        return distance_limit_range(v, opts.d_min, opts.d_max)
    except ValueError as exc:
        raise SystemExit(" --limit-distance: %s" % exc)


def resolve_av_angular_limits(opts):
    """Return validated AV-only angular sampling windows.

    The likelihood prior remains normalized on the full physical domain.  These
    limits therefore restrict the region integrated; they do not define a new,
    box-normalized prior.  Wraparound RA windows are intentionally refused until
    AV can represent their union as more than one hyperrectangle.
    """
    specs = (("limit_right_ascension", "ra", 0.0, 2.0 * np.pi),
             ("limit_declination", "dec", -0.5 * np.pi, 0.5 * np.pi),
             ("limit_psi", "psi", 0.0, np.pi),
             ("limit_inclination", "incl", 0.0, np.pi))
    out = {}
    for attr, name, physical_lo, physical_hi in specs:
        raw = getattr(opts, attr, None)
        if raw in (None, ""):
            continue
        try:
            parts = [float(item.strip()) for item in str(raw).split(",")]
        except ValueError:
            parts = []
        if len(parts) != 2 or not np.all(np.isfinite(parts)):
            raise SystemExit(" --limit-%s must be a finite LO,HI pair" %
                             attr[len("limit_"):].replace("_", "-"))
        lo, hi = parts
        if lo >= hi:
            raise SystemExit(" --limit-%s requires LO < HI (RA wraparound is "
                             "not supported)" %
                             attr[len("limit_"):].replace("_", "-"))
        if lo < physical_lo or hi > physical_hi:
            raise SystemExit(" --limit-%s must lie within [%g,%g]" %
                             (attr[len("limit_"):].replace("_", "-"),
                              physical_lo, physical_hi))
        out[name] = (lo, hi)
    return out


def av_distance_sampling_kwargs(like, d_lo, d_hi):
    """Apply distance limits only when distance is an sampled coordinate."""
    from RIFT.likelihood.jax_ile import samplers as _samplers
    if "distMpc" not in _samplers._av_param_order(like):
        return {}
    return {"sample_d_min": d_lo, "sample_d_max": d_hi}


def log_distance_box_correction(opts, with_distance):
    """ln[ prior mass on (d_min,d_max) / prior mass on the sampled box ]; 0.0 unless
    --limit-distance narrowed the box.

    sample_prior() draws distance from the prior RESTRICTED to the box, so its
    density is the prior divided by that mass.  The rule is per ESTIMATOR, not
    per function: an estimator that assumes "proposal == prior" must subtract
    this, and one that forms ln w = lnL + ln p - ln q explicitly must NOT,
    because log_prior() is already normalized over the physical [d_min,d_max].

    run_laplace_is contains ONE OF EACH, which is why this docstring no longer
    names functions.  Its adapted loop forms ln w explicitly and does not
    subtract; its prior pilot (kept as the collapse reference, #227) draws from
    sample_prior and therefore does.  Both then estimate the SAME integral on
    the SAME full-range normalization, which is the only reason the two are
    comparable and the 5-nat guard means anything.
    """
    if not with_distance:
        return 0.0
    lo, hi = resolve_distance_limit(opts)
    if lo == float(opts.d_min) and hi == float(opts.d_max):
        return 0.0
    return float(np.log((opts.d_max ** 3 - opts.d_min ** 3) / (hi ** 3 - lo ** 3)))


def sample_prior(n, opts, rng, with_distance):
    ra = rng.uniform(0.0, 2 * np.pi, n)
    dec = np.arcsin(rng.uniform(-1.0, 1.0, n))
    psi = rng.uniform(0.0, np.pi, n)
    incl = np.arccos(rng.uniform(-1.0, 1.0, n))
    phiref = rng.uniform(0.0, 2 * np.pi, n)
    cols = [ra, dec, psi, incl, phiref]
    if with_distance:
        u = rng.uniform(0.0, 1.0, n)
        # the SAMPLED range (== [d_min,d_max] unless --limit-distance)
        d_lo, d_hi = resolve_distance_limit(opts)
        dmin3, dmax3 = d_lo ** 3, d_hi ** 3
        cols.append((dmin3 + u * (dmax3 - dmin3)) ** (1.0 / 3.0))
    theta = np.stack(cols, axis=-1)
    return theta, log_prior(theta, opts, with_distance)


def log_prior(theta, opts, with_distance):
    ra, dec, psi, incl, phiref = [theta[..., i] for i in range(5)]
    inb = ((ra >= 0) & (ra <= 2 * np.pi) & (dec >= -np.pi / 2) & (dec <= np.pi / 2)
           & (psi >= 0) & (psi <= np.pi) & (incl >= 0) & (incl <= np.pi)
           & (phiref >= 0) & (phiref <= 2 * np.pi))
    with np.errstate(divide="ignore", invalid="ignore"):
        logp = (np.log(np.cos(dec)) - np.log(2.0)
                + np.log(np.sin(incl)) - np.log(2.0)
                - np.log(2 * np.pi) - np.log(np.pi) - np.log(2 * np.pi))
        if with_distance:
            dist = theta[..., 5]
            # SUPPORT is the sampled box (--limit-distance); the NORMALIZATION is
            # always over the physical [d_min,d_max], so lnZ keeps its full-range
            # scale and does not have to be corrected before comparison.
            d_lo, d_hi = resolve_distance_limit(opts)
            inb = inb & (dist >= d_lo) & (dist <= d_hi)
            dmin3, dmax3 = opts.d_min ** 3, opts.d_max ** 3
            logp = logp + np.log(3.0) + 2 * np.log(dist) - np.log(dmax3 - dmin3)
    return np.where(inb, logp, -np.inf)


# ---------------------------------------------------------------------------
# Batched lnL
# ---------------------------------------------------------------------------
def bounded_multipeak_config_from_options(opts):
    """Build the exact envelope executed by the wrapper from registered flags."""
    return validate_bounded_multipeak_config(BoundedMultipeakConfig(**{
        field: getattr(opts, "multipeak_jax_" + field)
        for field in BoundedMultipeakConfig._fields}))


def require_bounded_multipeak_rows(like, lnL):
    """Apply the requested decline action before publishing either artifact."""
    if getattr(like, "angle_marg_scheme", None) != "multipeak-jax":
        return
    bad = ~np.isfinite(np.asarray(lnL)) | (np.asarray(lnL) == BOUNDED_MULTIPEAK_LOG_ZERO)
    if getattr(like, "bounded_multipeak_decline_action", "refuse") == "drop":
        if np.all(bad):
            raise RuntimeError("multipeak-jax: no accepted output rows; change "
                               "the static envelope. No samples or evidence published.")
        return
    if np.any(bad) or getattr(like, "bounded_multipeak_declined", False):
        raise RuntimeError(
            "--angle-marg-scheme multipeak-jax: %d of %d rows could not be "
            "warranted within the static envelope; refusing to publish samples "
            "or evidence (including any earlier declined batch). Change the "
            "--multipeak-jax-* envelope, or select decline-action drop."
            % (int(np.sum(bad)), int(np.size(bad))))


def eval_lnL(like, theta, opts, with_distance):
    N = theta.shape[0]
    out = np.empty(N)
    # anglemarg schemes multiply the batch by (quad, dist_block, phi_chunk)
    # transients; --n-chunk 8000 is sized for the grid scheme.  Independent
    # slices, so the cap changes peak memory and nothing else.
    chunk = _angle_marg_eval_chunk(like, opts.n_chunk)
    for i in range(0, N, chunk):
        sl = slice(i, min(i + chunk, N))
        cols = [theta[sl, j] for j in range(theta.shape[1])]
        out[sl] = np.asarray(like.log_likelihood(*cols))
        if getattr(like, "bounded_multipeak_decline_action", "refuse") == "refuse":
            require_bounded_multipeak_rows(like, out[sl])
        if (getattr(like, "direct_marginalization_policy", "off") != "off"
                and np.any(np.isnan(out[sl]))):
            raise RuntimeError(
                "--direct-marginalization-policy %s: %d of %d evaluated rows "
                "could not be warranted (local gate declined AND the escalated "
                "band-limited reserve failed its guard/resolution check, or the "
                "norm table varies with time).  No coarse likelihood is "
                "substituted and the run is not published.  Raise "
                "--direct-marginalization-time-guard, raise the reserve "
                "refinement, or run with the policy off."
                % (like.direct_marginalization_policy,
                   int(np.sum(np.isnan(out[sl]))), int(sl.stop - sl.start)))
        if (getattr(like, "time_quadrature", "simpson") == "bandlimited"
                and np.any(np.isnan(out[sl]))):
            bad = np.where(np.isnan(out[sl]))[0]
            names = FULL_NAMES if with_distance else ANG_NAMES
            rows = "; ".join(
                ", ".join("%s=%.4f" % (nm, v) for nm, v in zip(names, theta[sl][j]))
                for j in bad[:3])
            # One uncertified row fails the run: no coarse likelihood is
            # substituted, and a row's contribution cannot be bounded without a
            # value.  With a full-sky prior the usual cause is a blind draw whose
            # detector arrival peak lies at or beyond the window edge (a wrong
            # sky shifts arrival by up to 2 R_earth/c ~ 43 ms), which the
            # trapezoid and guard certificates cannot converge on; the
            # integration window has to contain those shifts.  Measured rates:
            # DESIGN_jax_bandlimited_distmarg.md.
            raise RuntimeError(
                "--time-marginalization-quadrature bandlimited: %d of %d rows in "
                "this chunk failed a certificate (peak-width resolution, factor "
                "doubling, or guard agreement); no coarse likelihood is "
                "substituted and the run stops.  First failing rows: %s.  A "
                "peak at the window edge is the common cause under a full-sky "
                "prior; raise --data-integration-window-half to cover the "
                "detector arrival shifts (>= 2 R_earth / c = %.4f s plus the "
                "event-time uncertainty; 0.05 s measured clean), or raise the "
                "input/rholm sample rate."
                % (len(bad), int(sl.stop - sl.start), rows,
                   _BANDLIMITED_FULLSKY_HALF_WINDOW_MIN))
    return out


# ---------------------------------------------------------------------------
# Evidence helpers
# ---------------------------------------------------------------------------
def require_finite_evidence(logZ, neff, mode):
    """A NON-FINITE evidence is a FAILED event, not a result.

    The estimators return nan when the proposal never bracketed the peak (see
    ``_finalize_evidence``), and publishing that as an ILE row hands a
    downstream CIP fit a nan ``lnL`` for a template that simply was not
    integrated -- the same "leave no artifact behind" rule ``write_samples``
    already enforces for a cloud that admits no fair draw.  Raising here means
    ``--soft-fail-event-range`` skips to the next event, and without it the run
    exits nonzero: on the #227 configuration the shipped code exited 0.
    """
    if np.isfinite(logZ):
        return
    raise RuntimeError(
        "extrinsic integration produced a non-finite evidence (logZ=%r, "
        "neff=%.3g, mode=%s): the proposal did not bracket the likelihood peak, "
        "so no result row is written for this event.  Every --mode whose evidence "
        "comes from a single moment-matched Gaussian fitted to its own draws can "
        "fail this way on a narrow, high-SNR extrinsic posterior; try another "
        "--mode, and check its reported neff rather than only its lnZ."
        % (logZ, neff, mode))


# The prior pilot is an UNBIASED estimator of the same Z, not a bound on it, and
# the difference is load-bearing -- see prior_pilot_floor below.
PILOT_FLOOR_FP_RATE = 3.4e-4          # = exp(-8); see the operating curve below


def prior_pilot_floor(logZ_pilot, false_positive_rate=PILOT_FLOOR_FP_RATE):
    """A lower confidence bound on ln Z built from the prior pilot's estimate.

    WHY A BOUND AND NOT THE ESTIMATE.  The pilot is prior Monte Carlo: unbiased
    for Z, but with a heavy RIGHT tail once the target occupies a tiny fraction
    of the prior.  For a mode of prior mass ``m``, a single draw landing near the
    peak makes the estimate ~``L_max / n_pilot`` while the truth is ~``L_max m``,
    so it overshoots by ``1/(n_pilot m)`` -- more than ``T`` nats whenever
    ``n_pilot m < exp(-T)``.  Using the raw estimate as a floor therefore rejects
    a CORRECT adapted answer at exactly the rate that tail occurs.  That is not
    hypothetical here: at a synthetic width of 0.05 rad, ``n_pilot m = 1.6e-3``
    and the pilot was measured above the truth by up to +5.46 nats, P = 1.1e-3
    over 900 seeds.

    THE BOUND.  Markov is enough and needs nothing but unbiasedness and
    non-negativity, both of which hold: ``P(Zhat >= Z / rate) <= rate``.  So
    ``Zhat * rate`` is a lower confidence bound on Z at level ``1 - rate``, and
    the floor is ``ln Zhat + ln rate``.  The threshold is therefore a CHOSEN
    false-positive rate rather than a tuned constant, and it is distribution-free
    -- in particular it does NOT assume the pilot resolved anything.  It cannot:
    the pilot's own ESS is ~1 in every regime where this guard matters (measured
    median 1.0-1.3 for widths 0.03-0.08 rad), so "the pilot rests on one draw" is
    the normal state here, not an exceptional one.

    WHY exp(-8) AND NOT exp(-5).  Measured over 5400 synthetic runs (widths
    0.03-0.15 rad), sweeping the threshold T:

        T      Markov FP <=    FP measured   power on inaccurate runs
        5          6.7e-3            0/1180        0.651
        6          2.5e-3            0/1180        0.636
        8          3.4e-4            0/1180        0.606
       10          5.0e-5            0/1180        0.578

    Going from 5 to 8 costs 4.5 points of power and buys a 20x smaller worst-case
    false-positive rate; a false positive here FAILS THE EVENT and writes no row,
    so it is worth paying for.  No false positive was observed at any threshold:
    the largest pilot-minus-adapted gap on an accurate run was +1.654 nats, so 8
    clears the measured margin by 6.3 nats.  The bound, not the measurement, is
    what the guarantee rests on.
    """
    if not np.isfinite(logZ_pilot):
        return -np.inf
    return float(logZ_pilot) + float(np.log(false_positive_rate))


def evidence_from_logweights(logw):
    """(logZ, sigma/Z, neff) for Z = E[w] from log importance weights."""
    fin = np.isfinite(logw)
    logw = logw[fin]
    if logw.size == 0:
        return -np.inf, np.inf, 0.0
    m = np.max(logw)
    w = np.exp(logw - m)
    n = len(w)
    Zhat = np.mean(w)
    logZ = m + np.log(Zhat)
    sigma_over_Z = np.sqrt(np.var(w) / n) / Zhat
    neff = (np.sum(w) ** 2) / np.sum(w ** 2)
    return logZ, sigma_over_Z, neff


def _gaussian_logq(theta, mu, cov):
    diff = theta - mu[None, :]
    sol = np.linalg.solve(cov, diff.T).T
    return (-0.5 * np.einsum("ij,ij->i", diff, sol)
            - 0.5 * theta.shape[1] * np.log(2 * np.pi)
            - 0.5 * np.linalg.slogdet(cov)[1])


def _moment_match(theta, logL):
    w = np.exp(logL - np.max(logL))
    w = w / np.sum(w)
    mu = np.sum(w[:, None] * theta, axis=0)
    d = theta - mu[None, :]
    cov = (w[:, None, None] * d[:, :, None] * d[:, None, :]).sum(axis=0)
    cov += 1e-9 * np.eye(theta.shape[1]) * (np.trace(cov) / theta.shape[1] + 1e-12)
    return mu, cov


# ---------------------------------------------------------------------------
# Estimators
# ---------------------------------------------------------------------------
def run_prior_mc(like, opts, rng, dim, with_distance):
    theta, _ = sample_prior(opts.n_max, opts, rng, with_distance)
    lnL = eval_lnL(like, theta, opts, with_distance)
    # p_s == p (the proposal IS the prior), so ln w = lnL -- EXCEPT under
    # --limit-distance, where the proposal is the prior restricted to the box while
    # the prior keeps its full-range normalization; the constant below is exactly
    # 0.0 without that option, so the historical path is untouched.
    logw = lnL - log_distance_box_correction(opts, with_distance)
    logZ, sig, neff = evidence_from_logweights(logw)
    # These raw draws are PRIOR samples, not posterior ones: they must be
    # fair-drawn before export.
    return logZ, sig, neff, opts.n_max, theta, lnL, logw


def run_laplace_is(like, opts, rng, dim, with_distance, n_adapt=2):
    """Prior-seeded adaptive Gaussian importance sampling (Z = E_q[L p / q])."""
    n_pilot = min(opts.n_max // 4, 40000)
    theta_p, _ = sample_prior(n_pilot, opts, rng, with_distance)
    lnL_p = eval_lnL(like, theta_p, opts, with_distance)
    logL_p = lnL_p + log_prior(theta_p, opts, with_distance)
    mu, cov = _moment_match(theta_p, logL_p)
    # KEEP the pilot's own (prior-proposal) evidence estimate.  It is crude -- one
    # prior scan, often with an ESS of order 1 -- but it estimates the SAME integral
    # from a proposal that is guaranteed to cover the prior, and importance sampling
    # from a proposal that MISSES mass is biased low.  So an adapted estimate coming
    # out far BELOW this one is evidence that the adaptation walked away from the
    # peak: the failure mode that remains once #227's draw/density mismatch is fixed.
    # It is a REFERENCE, not a floor -- prior_pilot_floor() turns it into one, and
    # the distinction is the whole of external review's P1 on this PR.
    logZ_pilot, _, _ = evidence_from_logweights(
        lnL_p - log_distance_box_correction(opts, with_distance))

    per_round = max(opts.n_max // (n_adapt + 1), 1)
    all_theta, all_logw, all_lnL = [], [], []
    for r in range(n_adapt + 1):
        # ONE matrix for the draw AND the density.  Drawing from
        # ``cov_use + 1e-12*I`` and scoring under bare ``cov_use`` is issue #227:
        # once the adapted covariance falls below the absolute jitter the weights
        # are computed against a distribution that was never sampled, and this
        # mode returned lnZ = 5.8e9 on real O4 data while exiting 0.
        cov_use = _regularize_cov(cov * opts.proposal_inflate)
        Lc = np.linalg.cholesky(cov_use)
        z = rng.standard_normal((per_round, dim))
        theta = mu[None, :] + z @ Lc.T
        logq = _gaussian_logq(theta, mu, cov_use)
        logp = log_prior(theta, opts, with_distance)
        valid = np.isfinite(logp)
        lnL = np.full(per_round, -np.inf)
        if valid.any():
            lnL[valid] = eval_lnL(like, theta[valid], opts, with_distance)
        logw = np.where(valid, lnL + logp - logq, -np.inf)
        all_theta.append(theta); all_logw.append(logw); all_lnL.append(lnL)
        good = np.isfinite(logw)
        if r < n_adapt and good.sum() > 50:
            mu, cov = _moment_match(theta[good], logw[good])

    theta = np.concatenate(all_theta); logw = np.concatenate(all_logw)
    lnL = np.concatenate(all_lnL)
    logZ, sig, neff = evidence_from_logweights(logw)
    # log Z <= max lnL for a normalized prior, and a low neff means the proposal
    # never bracketed the peak.  Same rule the library samplers already apply;
    # this driver applied none, which is why #227 exited 0 on lnZ = 5.8e9.
    logZ, sig, neff = _finalize_evidence(
        logZ, sig, neff, float(np.max(lnL)) if np.isfinite(lnL).any() else np.nan)
    pilot_floor = prior_pilot_floor(logZ_pilot)
    if np.isfinite(logZ) and np.isfinite(pilot_floor) and logZ < pilot_floor:
        print("  [laplace-is] adapted proposal gives lnZ = %.3f, %.1f nats BELOW the "
              "prior pilot's Markov floor (%.3f, from a pilot estimate of %.3f at a "
              "%.1e false-positive rate): the adaptation moved off the peak, so this "
              "evidence is reported as unreliable."
              % (logZ, pilot_floor - logZ, pilot_floor, logZ_pilot,
                 PILOT_FLOOR_FP_RATE))
        logZ, sig = np.nan, np.nan
    # theta follows the GAUSSIAN PROPOSAL q, not the posterior; logw is what
    # turns it into one.  Returned so write_samples() can fair-draw.
    return logZ, sig, neff, len(theta), theta, lnL, logw


def run_nuts(like, opts, rng, with_distance):
    """Gradient-based NUTS over the distance-marginalized angular posterior.

    Demonstrates the AD payoff: NUTS uses the exact JAX gradient of the
    (smooth, distance-marginalized) log-likelihood to climb and sample the
    narrow extrinsic peak efficiently -- where brute-force MC stalls.  The
    posterior samples then seed a high-neff importance estimate of the evidence
    (a Gaussian proposal moment-matched to the NUTS draws), and are written out
    as the per-sample product (the natural seed for a flowMC / AV closeout).
    """
    if with_distance:
        raise SystemExit("--mode nuts requires distance marginalization (the "
                         "bare 6-D likelihood is degenerate); pass "
                         "--distance-marginalization.")
    import numpyro
    import numpyro.distributions as dist
    from numpyro.infer import MCMC, NUTS
    import jax.numpy as jnp

    two_pi = float(2 * np.pi)

    def model():
        ra = numpyro.sample("ra", dist.Uniform(0.0, two_pi))
        sin_dec = numpyro.sample("sin_dec", dist.Uniform(-1.0, 1.0))
        psi = numpyro.sample("psi", dist.Uniform(0.0, float(np.pi)))
        cos_incl = numpyro.sample("cos_incl", dist.Uniform(-1.0, 1.0))
        phiref = numpyro.sample("phiref", dist.Uniform(0.0, two_pi))
        dec = jnp.arcsin(sin_dec)
        incl = jnp.arccos(cos_incl)
        lnL = like._scalar(jnp.stack([ra, dec, psi, incl, phiref]))
        numpyro.factor("loglike", lnL)

    # Seed NUTS at the best prior draw so it climbs the GLOBAL mode rather than
    # one of the (many) secondary sky modes from the detector time-delay ring.
    seed_theta, _ = sample_prior(4000, opts, rng, with_distance=False)
    seed_lnL = eval_lnL(like, seed_theta, opts, with_distance=False)
    th0 = seed_theta[np.argmax(seed_lnL)]
    init_vals = {
        "ra": float(th0[0]), "sin_dec": float(np.sin(th0[1])),
        "psi": float(th0[2]), "cos_incl": float(np.cos(th0[3])),
        "phiref": float(th0[4]),
    }
    print("  NUTS seeded at best-of-%d prior draw (lnL=%.2f)"
          % (len(seed_theta), seed_lnL.max()))
    from numpyro.infer import init_to_value
    key = jax.random.PRNGKey(opts.seed)
    kernel = NUTS(model, target_accept_prob=0.8,
                  init_strategy=init_to_value(values=init_vals))
    mcmc = MCMC(kernel, num_warmup=opts.num_warmup,
                num_samples=opts.num_samples, num_chains=opts.num_chains,
                progress_bar=opts.verbose)
    mcmc.run(key)
    s = mcmc.get_samples()
    ra = np.asarray(s["ra"]); dec = np.arcsin(np.asarray(s["sin_dec"]))
    psi = np.asarray(s["psi"]); incl = np.arccos(np.asarray(s["cos_incl"]))
    phiref = np.asarray(s["phiref"])
    theta = np.stack([ra, dec, psi, incl, phiref], axis=-1)
    lnL = eval_lnL(like, theta, opts, with_distance=False)
    print("  NUTS drew %d samples; lnL in [%.2f, %.2f], peak at:"
          % (len(theta), lnL.min(), lnL.max()))
    for nm, v in zip(ANG_NAMES, theta[np.argmax(lnL)]):
        print("     %-8s = %.5f" % (nm, v))

    # NUTS-seeded importance estimate of the evidence (Gaussian proposal
    # moment-matched to the posterior draws -> high neff vs prior seeding).
    mu, cov = _moment_match(theta, np.zeros(len(theta)))  # posterior moments
    n_is = min(opts.n_max, 40000)
    cov_is = _regularize_cov(cov * opts.proposal_inflate)   # one matrix (#227)
    Lc = np.linalg.cholesky(cov_is)
    z = rng.standard_normal((n_is, 5))
    th_is = mu[None, :] + z @ Lc.T
    logq = _gaussian_logq(th_is, mu, cov_is)
    logp = log_prior(th_is, opts, with_distance=False)
    valid = np.isfinite(logp)
    lnL_is = np.full(n_is, -np.inf)
    if valid.any():
        lnL_is[valid] = eval_lnL(like, th_is[valid], opts, with_distance=False)
    logw = np.where(valid, lnL_is + logp - logq, -np.inf)
    logZ, sig, neff = evidence_from_logweights(logw)
    # Bound with the NUTS CHAIN's peak, not the IS cloud's.  Both are valid upper
    # bounds on ln Z, the chain's is the larger (it sits ON the peak), and using the
    # smaller one would fail a correct run whose IS cloud happened to be broad.  This
    # is also the argument samplers.py passes at its own _finalize_evidence sites.
    logZ, sig, neff = _finalize_evidence(
        logZ, sig, neff, float(np.max(lnL)) if np.isfinite(lnL).any() else np.nan)
    # theta/lnL are the NUTS chain (already targets the posterior); the IS cloud
    # th_is/logw is only the evidence estimator, so there is nothing to reweight.
    # `neff` below therefore describes th_is, NOT the exported chain -- it must
    # never reach the export count.  It used to: fairdraw_size clamped the export
    # by 1.5*neff, truncating this chain by an unrelated cloud's ESS on what is
    # the driver's DEFAULT mode.  write_samples no longer accepts a neff at all.
    return logZ, sig, neff, n_is, theta, lnL, None


def run_map(like, opts, rng, dim, with_distance):
    """Gradient-ascend the (well-conditioned, distance-marginalized) peak."""
    from scipy.optimize import minimize
    ang_bounds = [(0, 2 * np.pi), (-np.pi / 2 + 1e-3, np.pi / 2 - 1e-3),
                  (0, np.pi), (1e-3, np.pi - 1e-3), (0, 2 * np.pi)]
    bounds = ang_bounds + ([resolve_distance_limit(opts)] if with_distance else [])
    theta_seed, _ = sample_prior(4000, opts, rng, with_distance)
    lnL_seed = eval_lnL(like, theta_seed, opts, with_distance)
    x0 = theta_seed[np.argmax(lnL_seed)]

    def negf(x):
        v, g = like.value_and_grad(x)
        return -float(v), -np.asarray(g)

    res = minimize(negf, x0, jac=True, method="L-BFGS-B", bounds=bounds,
                   options={"maxiter": 400})
    return res.x, -res.fun


# ---------------------------------------------------------------------------
# Output
# ---------------------------------------------------------------------------
def dat_path(opts, out_index):
    """Path of the normal per-event ILE result row (evidence, sigma, neff)."""
    return opts.output_file + "_" + str(out_index) + "_" + ".dat"


def samples_path(opts, out_index):
    """Path of this event's exported extrinsic samples."""
    return opts.output_file + "_" + str(out_index) + "_samples.dat"


def direct_marginalization_policy_note(like, theta, n_max=256, chunk=32,
                                       return_values=False):
    """Label the branch statistics of the exported rows under the policy.

    Returns "" when the policy is off.  Otherwise evaluates the ledger on up to
    ``n_max`` evenly spaced exported rows and reports how many took the local
    branch, the warranted reserve, or came back unusable (reserve failed its
    own warrant; the value is retained and the run is labelled).  The note is
    written into the sample and evidence headers next to the angle-grid label.
    """
    policy = getattr(like, "direct_marginalization_policy", "off")
    ledger_fn = getattr(like, "_batched_ledger", None)
    if policy == "off" or ledger_fn is None:
        return ("", None, None) if return_values else ""
    theta = np.asarray(theta)
    if theta.ndim != 2 or theta.shape[0] == 0:
        _empty = "DIRECT-MARG-POLICY=%s rows=0" % policy
        return (_empty, None, None) if return_values else _empty
    n = int(theta.shape[0])
    idx = np.unique(np.linspace(0, n - 1, min(n, int(n_max))).astype(int))
    parts = []
    for start in range(0, idx.size, int(chunk)):
        sub = theta[idx[start:start + int(chunk)]]
        cols = [jnp.asarray(sub[:, j]) for j in range(sub.shape[1])]
        _, ledger = ledger_fn(*cols)
        parts.append({k: np.asarray(v) for k, v in ledger.items()})
    ledger = {k: np.concatenate([p[k] for p in parts]) for k in parts[0]}
    summary = _summarize_policy_ledger(ledger)
    declines = ",".join("%s:%d" % (k.replace("decline_", ""), v)
                        for k, v in sorted(summary["declines"].items()))
    note = ("DIRECT-MARG-POLICY=%s rows=%d local=%d reserve=%d "
            "warranted-reserve=%d escalations=%d unusable=%d reconciles=%d "
            "declines=[%s] max-local-error-score-nats=%.3g"
            % (policy, summary["rows"], summary["accepted_local"],
               summary["reserve_executed"], summary["reserve_warranted"],
               summary.get("reserve_escalations", 0),
               summary["unusable"], summary["reconciles"], declines,
               summary["max_local_error_score_nats"]))
    if return_values:
        # The per-row selected value, so two probes at ONE seed -- one whose
        # config accepts, one forced onto the reserve -- can be differenced on
        # exactly the same prior draws.  An accepted row never evaluates the
        # reserve, so this is the only way to price the local branch against it.
        return (note, np.asarray(ledger["lnL"], dtype=float),
                np.asarray(ledger["accepted_local"], dtype=bool))
    return note


def angle_grid_suspect_note(scheme=None):
    """Label describing the angle-grid amplitude check for this event.

    Returns one of three things, and the THIRD states the exact coverage:

      ""                              -- the grid schemes were not used
      "SUSPECT-ANGLE-GRID ..."        -- undersizing was DETECTED
      "ANGLE-GRID-CHECK=OUTPUT-CLOUD-PASS ..." -- checked cloud was adequate
      "ANGLE-GRID-CHECK=NOT-PERFORMED ..."     -- amp-sized scheme, NO check ran

    The pure JIT returns its amplitude metric as ordinary data, and the Python
    boundary synchronously accumulates the maximum over every pilot, reweight,
    and final production/output-cloud batch.  That deterministic coverage is
    enough to label the points used for the published flow evidence and sample
    export.  Transient flow-training-only proposals (which do not enter those
    artifacts) are deliberately NOT claimed.  This split also keeps host
    callbacks out of the expensive graph so JAX can persistently cache it.
    """
    if scheme == "multipeak-jax":
        return ("MULTIPEAK-JAX bounded_cost=true dense_reserve=false "
                "fixed_plan_autodiff_only=true derivative_warrant_certified=false")
    st = _anglemarg.amp_failsafe_state()
    if st.get("tripped"):
        return ("SUSPECT-ANGLE-GRID amp_failsafe=TRIPPED worst_amp=%.6g "
                "amp_sizing=%.6g scheme=%s"
                % (st.get("worst_amp", float("nan")),
                   st.get("amp_sizing", float("nan")), st.get("scheme")))
    # 'peak-local' belongs here for the same reason it belongs in the failsafe: its phi
    # axis is dense and amp-sized, so its artifacts are entitled to no more confidence
    # than the other two, and a scheme missing from this list would publish output with
    # NO standing label at all -- the silence a reader would read as verification.
    # 'phi-local' belongs here too, and its omission was found by external review rather
    # than by this comment being read.  It runs _runtime_amp_failsafe and evaluates an
    # amp-sized dense fallback on every row, so it is amplitude-sized in exactly the sense
    # this label is about; leaving it out published its artifacts with an EMPTY note --
    # the silence the paragraph above says a reader would take for verification.
    if scheme in ("exact", "laplace", "peak-local", "phi-local"):
        # A PASS may only be claimed when a batch was actually inspected.  The
        # recorder is wired ONLY for direct_marginalization_policy="off"
        # (wrapper.py sets self._amp_record under that condition), while
        # self.angle_marg_scheme still names the amp-sized reserve scheme.  So
        # `--angle-marg-scheme exact --direct-marginalization-policy auto`
        # reaches here with n_calls == 0 and amp_sizing None.  Without this
        # branch that formatted "%.6g" % None and raised TypeError, killing the
        # event AFTER integration and BEFORE either writer -- and had amp_sizing
        # merely been 0.0 it would instead have published OUTPUT-CLOUD-PASS
        # worst_amp=0, an affirmative adequacy claim backed by zero checks,
        # which is the exact false negative this label exists to prevent.
        if not st.get("n_calls"):
            return ("ANGLE-GRID-CHECK=NOT-PERFORMED scheme=%s (no batched "
                    "amp-sized evaluation was recorded for this event, so the "
                    "adequacy of the dense (phi,psi) grids is UNKNOWN -- this "
                    "is NOT a pass; rebuild at a larger amp_sizing if it "
                    "matters)" % (scheme,))
        return ("ANGLE-GRID-CHECK=OUTPUT-CLOUD-PASS worst_amp=%.6g "
                "amp_sizing=%.6g scheme=%s (deterministic over pilot/reweight/"
                "final output-cloud evaluations; transient training-only "
                "proposals not inspected)"
                % (st.get("worst_amp", float("nan")),
                   st.get("amp_sizing", float("nan")), st.get("scheme")))
    return ""


def write_dat(opts, P, out_index, event_id, logZ, sigma_over_Z, ntotal, neff,
              angle_note=""):
    if not opts.output_file:
        return
    m1, m2 = P.m1 / MSUN, P.m2 / MSUN
    fname = dat_path(opts, out_index)
    row = np.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y,
                     P.s2z, logZ, sigma_over_Z, ntotal, neff]])
    hdr = "event_id m1 m2 s1x s1y s1z s2x s2y s2z lnL sigma_lnL ntotal neff"
    # The EVIDENCE artifact must carry the label independently of sample export:
    # write_samples() early-returns without --save-samples/--output-file, so a
    # run with export disabled would otherwise publish a numeric evidence row
    # indistinguishable from a clean integration.
    if angle_note:
        hdr += "\n" + angle_note
    np.savetxt(fname, row, header=hdr)
    print("Wrote %s" % fname)


def tempered_cloud_size(opts, n_starts):
    """Rows flowMC hands to the exporter: chains x loops x (local + global) steps.

    Verified against a real run: 20 chains, 6 production loops, 20 local + 20
    global steps -> 4800, the ntotal the driver reported.
    """
    n_chains = max(n_starts, 20)
    return int(n_chains * opts.n_production_loops
               * (opts.n_local_steps + opts.n_global_steps))


def record_supplied_options(opts, argv, parser=None):
    """Record which long option TOKENS appeared on the command line.

    Everything below used to infer "did the user pass this?" from "does its value
    differ from the default?".  That is wrong whenever the user explicitly passes
    the default: `--auto-adapt-weight-exponent --adapt-weight-exponent 1.0` was
    silently ACCEPTED and the chooser then replaced the user's explicit untempered
    target, which is the opposite of the documented conflict behaviour.  Handles
    both `--opt value` and `--opt=value`.
    """
    supplied = set()
    for tok in (argv if argv is not None else sys.argv[1:]):
        if not (isinstance(tok, str) and tok.startswith("--")):
            continue
        name = tok.split("=", 1)[0]
        # CANONICALISE.  optparse accepts unambiguous long-option PREFIXES, so
        # `--adapt-weight-exp=1.0` parses and sets the value while recording a
        # token that no was_supplied() query matches -- auto then silently
        # overwrote the user's explicit exponent.  Ask the parser what the token
        # actually resolved to.  Reported by review on #186.
        if parser is not None:
            try:
                name = parser._match_long_opt(name)
            except Exception:
                pass          # unknown/ambiguous: record it verbatim
        supplied.add(name)
    opts._supplied_options = supplied
    return supplied


def was_supplied(opts, flag):
    """True when ``flag`` was named on the command line.

    Falls back to False when the record is absent (a caller that built an options
    object directly rather than parsing), which is the safe direction: it means
    "assume the user did not pass it" and so never fabricates a conflict.
    """
    return flag in getattr(opts, "_supplied_options", set())


def resolve_ile_interface_aliases(opts, optp):
    """Resolve conventional ILE spellings into JAX-native option values."""
    if getattr(opts, "interpolate_time", None) is not None:
        ile_interp = str(opts.interpolate_time).strip().lower()
        # Conventional ILE accepted this as a boolean before it accepted the
        # explicit stencil names.  Production args_ile.txt files containing
        # either spelling must remain swappable into this driver.
        if ile_interp in ("1", "true", "t", "yes", "y", "on"):
            ile_interp = "cubic"
        elif ile_interp in ("0", "false", "f", "no", "n", "off", "none"):
            ile_interp = "nearest"
        if ile_interp not in _JAX_GATHERER_NAMES:
            optp.error("--interpolate-time must be one of %s or a legacy boolean" %
                       ", ".join(sorted(_JAX_GATHERER_NAMES)))
        if was_supplied(opts, "--interp") and opts.interp != ile_interp:
            optp.error("--interp %r and --interpolate-time %r disagree" %
                       (opts.interp, ile_interp))
        opts.interp = ile_interp
    return opts


def _normalize_interpolate_time_argv(argv):
    """Give conventional ILE's historical bare flag its ``True`` value."""
    raw = list(sys.argv[1:] if argv is None else argv)
    out = []
    for i, token in enumerate(raw):
        out.append(token)
        if (token == "--interpolate-time" and
                (i + 1 == len(raw) or raw[i + 1].startswith("--"))):
            out.append("True")
    return out


def _target_ess_was_given(opts):
    """True when --target-export-ess-frac was named on the command line."""
    return was_supplied(opts, "--target-export-ess-frac")


def resolve_tempering_exponent(opts, n_dim, n_cloud):
    """Settle --adapt-weight-exponent for a tempered mode, and refuse a dead one.

    Two jobs, both of which exist because beta means something different here
    than it does in non-JAX RIFT (jax_ile/DESIGN_jax_tempering.md):

      * ``--auto-adapt-weight-exponent`` picks beta from the EXPORT budget via
        ``samplers.beta_for_export_ess``, i.e. from the sampled dimension.  There
        is no SNR term, deliberately: the reweight cost
        ``ESS/N = [beta(2-beta)]^(n_dim/2)`` is independent of lnLmax, so the
        non-JAX helper's SNR-keyed rule would be keying on the wrong variable.
      * whatever beta ends up in force, the predicted export ESS is REPORTED, and
        a beta that cannot support a usable export RAISES.  Printing a warning
        above a 199-row file is the silent degradation this exists to remove.

    ``--adapt-adapt`` anneals inv_T up to 1 and therefore exports at full ESS;
    it is left alone here.

    RETURNS the exponent to use and does NOT write it back to ``opts``.  It used
    to assign ``opts.adapt_weight_exponent``, which gave per-RUN state a per-EVENT
    meaning: ``analyze_one`` is called once per intrinsic template with the SAME
    opts, so on event 1 the chooser read its own event-0 output as a
    user-supplied exponent and aborted the batch with SystemExit.  Reproduced
    with ``--n-events-to-analyze 3``; ILE_extr.sub runs batches.
    """
    from RIFT.likelihood.jax_ile.samplers import (
        beta_for_export_ess, export_ess_estimate)

    if opts.adapt_adapt and opts.mode not in _ADAPT_ADAPT_MODES:
        # samplers.flowmc_sample (the plain 5-D `flowmc` mode) has NO temper_adapt
        # argument -- it takes `temper` and samples a STATIC beta.  Returning the
        # "annealed to beta=1, full ESS" answer here reported annealing that never
        # happened AND skipped the degenerate-export check, so
        # `--mode flowmc --adapt-adapt --adapt-weight-exponent 0.1` sampled at
        # beta=0.1 unguarded while claiming full ESS.  Reported by review on #186.
        raise SystemExit(
            "--adapt-adapt is not implemented for --mode %s: its sampler "
            "(samplers.flowmc_sample) takes a STATIC tempering exponent and has "
            "no annealing path, so the flag would be silently ignored while the "
            "run reported annealing.  Use one of %s, or drop --adapt-adapt and "
            "set --adapt-weight-exponent explicitly."
            % (opts.mode, " ".join(sorted(_ADAPT_ADAPT_MODES))))

    if opts.adapt_adapt:
        if opts.auto_adapt_weight_exponent:
            raise SystemExit(
                "--auto-adapt-weight-exponent and --adapt-adapt both set.  The "
                "anneal already finishes at beta=1 (full export ESS), so there "
                "is no static exponent for the chooser to pick.  Use one.")
        print("Tempering: --adapt-adapt (anneal inv_T -> 1); export is untempered, "
              "full ESS.")
        return float(opts.adapt_weight_exponent)

    # --smc-puffball routes these modes to samplers.smc_puffball_sample, whose
    # signature swallows `temper` in **_ignore: it runs its OWN adaptive SMC
    # ladder, so --adapt-weight-exponent has NO effect there.  Guarding it would
    # REFUSE a run over a number that does nothing.  Report the no-op instead; a
    # silent one is what this driver's compat layer exists to prevent.
    #
    # (That sampler now reports the inv_T its ladder actually reached and a
    # matching L^(1-inv_T) post_weight -- so its export weights are uniform only
    # because analyze_one REFUSES the event unless the ladder finished at
    # inv_T=1.  The exponent being inert here does not depend on that, but do not
    # re-derive "its weights are always uniform" from this branch.)
    if getattr(opts, "smc_puffball", False):
        if opts.auto_adapt_weight_exponent or was_supplied(opts, "--adapt-weight-exponent"):
            print("Note: --smc-puffball ignores --adapt-weight-exponent / "
                  "--auto-adapt-weight-exponent (smc_puffball_sample runs its own "
                  "SMC temperature ladder; the event is refused unless that "
                  "ladder finishes at inv_T=1).")
        return 1.0

    if not opts.auto_adapt_weight_exponent and _target_ess_was_given(opts):
        # Setting a target and no chooser does nothing at all.  Say so rather than
        # let the request evaporate -- the same reason the ILE compat layer reports
        # accepted-but-ignored options.
        print("Note: --target-export-ess-frac %g has no effect without "
              "--auto-adapt-weight-exponent." % opts.target_export_ess_frac)

    if opts.auto_adapt_weight_exponent:
        # An explicit --adapt-weight-exponent alongside --auto is a contradiction.
        # Silently overriding it would be the worst of both: the run reports a
        # chooser it did not obey the user about.
        if was_supplied(opts, "--adapt-weight-exponent"):
            raise SystemExit(
                "--auto-adapt-weight-exponent was given together with an explicit "
                "--adapt-weight-exponent %g.  The chooser would overwrite it.  Pass "
                "one or the other.  (Detected from the command line, not from the "
                "value: passing the default explicitly, --adapt-weight-exponent 1.0, "
                "is still an explicit choice of an untempered target.)"
                % float(opts.adapt_weight_exponent))
        beta = beta_for_export_ess(opts.target_export_ess_frac, n_dim)
        print("Tempering: AUTO beta=%.5f for %.0f%% export ESS in %d-D "
              "(ESS/N=[beta(2-beta)]^(dim/2); no SNR term -- see "
              "jax_ile/DESIGN_jax_tempering.md)"
              % (beta, 100.0 * opts.target_export_ess_frac, n_dim))
    else:
        beta = float(opts.adapt_weight_exponent)

    if beta > 1.0:
        # NOT harmless, and NOT "untempered": samplers.flowmc_sample* take
        # temper = 1/beta, so beta>1 gives inv_T>1 -- it SHARPENS the target past
        # the posterior, and the export reweight L^(1-beta) has
        # ESS/N = [beta(2-beta)]^(dim/2), which is 0 at beta=2 and undefined
        # beyond.  An earlier version of this branch tested `beta >= 1.0` and
        # printed "beta=1 (untempered target)" here, which was false for every
        # beta > 1.
        raise SystemExit(
            "--adapt-weight-exponent %g is greater than 1.  On this path the "
            "exponent is applied to the SAMPLED target (inv_T = %g), so a value "
            "above 1 samples a target sharper than the posterior and the export "
            "reweight L^(1-beta) diverges.  Use beta in (0, 1]." % (beta, beta))
    if beta <= 0.0:
        raise SystemExit(
            "--adapt-weight-exponent %g must be positive: beta=0 samples the "
            "prior and the export reweight carries the entire likelihood."
            % beta)
    if beta == 1.0:
        print("Tempering: beta=1 (untempered target); export ESS is the full cloud.")
        return beta
    # An ESTIMATE, calibrated at SNR ~= 23.8 -- NOT a bound.  This used to REFUSE
    # a run below the floor, which is a guarantee the number cannot support: the
    # SNR ladder measures 0.00823 at beta=0.1/d=4/SNR~67 against an estimate of
    # 0.0285.  Refusing on it would be false precision in one direction and false
    # confidence in the other, so it now WARNS.  Reported by review on #186.
    frac = export_ess_estimate(beta, n_dim)
    ess = frac * n_cloud
    print("Tempering: beta=%.5f in %d-D -> export ESS/N >= %.4f "
          "(measured-calibrated lower bound), "
          "ESS~%.0f of %d rows" % (beta, n_dim, frac, ess, n_cloud))
    if ess < _USABLE_EXPORT_ESS and not opts.allow_degenerate_tempering:
        print(
            "  WARNING: --adapt-weight-exponent %g gives an ESTIMATED export ESS "
            "of %.0f (< %d), so the reweighted --save-samples cloud is unlikely "
            "to be a usable posterior sample.\n"
            "  This is the trap the non-JAX helper's rule sets here: it picks "
            "beta from the SNR (beta=0.1 at SNR<=22.5, 0.1*(22.5/SNR)^2 above), "
            "which is correct where beta only shapes a PROPOSAL, but on this "
            "path beta is the exponent of the SAMPLED target.\n"
            "  The estimate is calibrated at SNR ~= 23.8 and is OPTIMISTIC at "
            "higher SNR (measured 0.00823 vs an estimated 0.0285 at beta=0.1, "
            "d=4, SNR ~= 67), so the true figure may be several times worse.  "
            "It is a warning rather than a refusal precisely because it cannot "
            "support a hard floor; check the ESS the export actually reports."
            % (beta, ess, _USABLE_EXPORT_ESS), file=sys.stderr)
    return beta


def fairdraw_indices(logw, rng):
    """Indices that turn a WEIGHTED cloud into an equal-weight one, or ``None``.

    Production ILE convention (``RIFT/integrators/mcsampler.py::integrate`` and
    the identical block in mcsamplerGPU / mcsamplerAdaptiveVolume /
    mcsamplerEnsemble / mcsamplerPortfolio): normalise ``w = L * p / p_s`` and
    multinomial-resample WITH replacement, so the exported rows are equal weight
    and carry no weight column.  Capped at ``1.5*ESS``: never claim more fair
    draws than the weights support, or the file looks like N independent draws
    while holding ~ESS distinct points.

    This function does REWEIGHTING ONLY.  The export count requested by
    ``--n-fairdraw-extrinsic-samples`` / ``--fairdraw-extrinsic-output-n-max``
    is applied by the caller, on every path -- including this one's ``None``
    returns -- so that the count contract has exactly one implementation and
    cannot be quietly skipped for some configurations.

    Returns ``(indices_or_None, note)``; the note always records the ESS.  A
    ``FAILED`` note means the cloud CANNOT be fair-drawn: ``write_samples``
    refuses to write an export at all in that case, rather than passing the
    unreweighted cloud off as one.
    """
    logw = np.asarray(logw, dtype=float)
    fin = np.isfinite(logw)
    # FAIL LOUDLY, NOT OPEN.  "cannot compute weights" and "weights are already
    # uniform" both used to return None, so a degenerate weight vector silently
    # wrote the UNCORRECTED cloud under a header that promises a fair draw.
    if fin.sum() < 2:
        return None, ("FAILED: %d of %d weights are finite ESS=n/a n_in=%d"
                      % (int(fin.sum()), len(logw), len(logw)))
    lw = logw[fin] - np.max(logw[fin])
    w = np.exp(lw)
    tot = w.sum()
    if not np.isfinite(tot) or tot <= 0:
        return None, ("FAILED: weight sum is %r (overflow or all-zero) "
                      "ESS=n/a n_in=%d" % (tot, len(logw)))
    w = w / tot
    neff = 1.0 / np.sum(w ** 2)
    idx_fin = np.where(fin)[0]
    if np.allclose(w, w[0]):
        # Already equal weight -- the DEFAULT for the flowMC modes, since
        # --adapt-weight-exponent defaults to 1.  Nothing to reweight; the
        # caller still applies any requested count.
        # NOT an effective sample size: the weights carry no information about
        # how correlated the underlying chain is, so reporting the row count as
        # "ESS" would be a number that looks measured and is not.
        return None, "none (weights uniform) ESS=n/a n_in=%d" % (len(logw),)
    n_out = int(max(1, min(int(np.ceil(1.5 * neff)), len(logw))))
    print("  fairdraw: %d weighted samples (ESS=%.1f) -> %d equal-weight draws"
          % (len(logw), neff, n_out))
    if neff < _USABLE_EXPORT_ESS:
        print("  fairdraw: WARNING ESS=%.1f -- the proposal barely covers this "
              "posterior; the exported cloud is NOT a usable posterior sample "
              "however it is drawn." % neff, file=sys.stderr)
    note = "reweighted ESS=%.1f n_in=%d" % (neff, len(logw))
    return idx_fin[rng.choice(len(idx_fin), size=n_out, replace=True, p=w)], note


def fairdraw_size(opts, n_have):
    """Requested number of fair draws, or ``None`` for "as many as the weights
    support".

    ``--n-fairdraw-extrinsic-samples`` is an exact request;
    ``--fairdraw-extrinsic-output-n-max`` caps the count per evaluation.  Both
    are clamped by ``n_have`` and by nothing else.

    WHY THERE IS NO ESS TERM HERE, though ILE appears to have one.
    ``mcsampler.integrate`` clamps with
    ``n_extr = min(n_extr, 1.5*eff_samp, 1.5*neff)``, but there ``eff_samp``,
    ``neff`` and the fair draw are all properties of the SAME importance-weight
    vector ``wt`` that the draw then samples from -- one weight stream, so the
    clamp is self-consistent.

    That does not hold on this path.  For every flowMC mode ``theta`` is the
    production chain while ``res["neff"]`` comes from a SEPARATE estimator: the
    moment-matched Gaussian importance cloud built for the evidence, or (under
    ``--adapt-adapt``) the annealing ladder's minimum inter-stage ESS.  Neither
    describes the exported rows.  Clamping the export by it truncated a
    perfectly good equal-weight chain to a couple of rows whenever that
    unrelated proposal happened to have low ESS.

    The ESS that DOES bound this export is the one ``fairdraw_indices``
    computes from the export weights -- and it has already applied it, as
    ``n_out = min(ceil(1.5*ESS), n_in)``.  So ``n_have`` already carries the
    1.5*ESS cap wherever weights were non-uniform, and where they were uniform
    there is no computed chain ESS to clamp by (an autocorrelation-based one
    would be a different, unimplemented quantity).  Reported by review on #180."""
    n_req = getattr(opts, "n_fairdraw_extrinsic_samples", None)
    if n_req is None and getattr(opts, "fairdraw_extrinsic_output", False):
        # ILE's default cap is 5; kept out of the parser so an unset flag is not
        # reported as one the user passed.
        n_req = getattr(opts, "fairdraw_extrinsic_output_n_max", None)
        if n_req is None:
            n_req = _FAIRDRAW_N_MAX_DEFAULT
    if n_req is None:
        return None
    n_asked = int(n_req)
    n_req = max(1, min(n_asked, n_have))
    if n_req != n_asked:
        # ILE prints "Fairdraw size : n" whenever it clamps; silence here meant a
        # laplace-is run quietly turned a request for 137 into 32.
        print("  Fairdraw size : %d (requested %d, clamped by the %d available "
              "rows; those already carry the 1.5*ESS cap from fairdraw_indices "
              "when the weights were non-uniform)" % (n_req, n_asked, n_have))
    return n_req


def _remove_stale_artifact(path, what="export"):
    """Delete a leftover file at a path this failed event refuses to publish.

    The run was going to overwrite ``path``; if it is left behind, a re-run in
    an existing output directory leaves the PREVIOUS run's file where the
    pipeline looks for this one's, which is the same "read a file that is not
    this event's" failure the refusal exists to prevent."""
    if not os.path.exists(path):
        return
    try:
        os.remove(path)
    except OSError as e:
        print("  *** could not remove the stale %s %s (%s) -- it is NOT "
              "this run's output. ***" % (what, path, e), file=sys.stderr)
    else:
        print("  removed the stale %s %s (this event produced no valid result)"
              % (what, path), file=sys.stderr)


def write_samples(opts, out_index, theta, lnL, with_distance, logw=None,
                  angle_note=""):
    """Write the exported extrinsic samples.

    ``logw`` are per-sample LOG IMPORTANCE WEIGHTS ``ln(L p / p_s)`` for the
    rows of ``theta`` (``None`` when the sampler already targets the posterior,
    e.g. an untempered MCMC chain).  When they are non-uniform the cloud is
    fair-drawn against them BEFORE writing, so the exported rows are equal
    weight -- the same contract production ILE's ``--fairdraw-extrinsic-output``
    provides, and the one every downstream consumer of these files assumes.

    If the fair draw CANNOT be performed (degenerate or unnormalizable weights)
    no file is written and ``RuntimeError`` is raised: the alternative is to ship
    the raw proposal/prior cloud under the product name that means "posterior
    draws".  ``analyze_one`` calls this BEFORE ``write_dat``, so that refusal
    also leaves the event without a normal ILE result row.
    """
    if not (opts.output_file and opts.save_samples) or theta is None:
        return
    sname = samples_path(opts, out_index)
    # The export RNG is derived here and NOWHERE ELSE.  It must never be the
    # generator that feeds the samplers/estimators: --save-samples is an OUTPUT
    # flag and consuming the science stream made it change the lnL/logZ of every
    # later event in a batch.  There is deliberately no rng parameter, so that
    # mistake cannot be reintroduced by a caller.
    rng = np.random.default_rng((opts.seed, out_index))
    # n_in is the cloud as the sampler handed it over, before any filtering or
    # counting, so the header records what went in as well as what came out.
    note = ("not applicable (sampler targets the posterior) ESS=n/a n_in=%d"
            % len(theta))
    if logw is not None and len(logw) != len(theta):
        # FAIL LOUDLY, NOT OPEN.  This used to be a silent `and len(...) ==` in
        # the guard below, so a mismatched weight vector fell through to the
        # "not applicable (sampler targets the posterior)" note -- writing an
        # UNREWEIGHTED tempered cloud under a header claiming none was needed.
        # A length mismatch is a sampler bug, not a configuration; say so.
        raise RuntimeError(
            "export weights and samples disagree in length (%d weights, %d "
            "rows) for --mode %s.  Refusing to write: the fair draw cannot be "
            "performed and labelling the raw cloud 'not applicable' would be a "
            "false provenance line."
            % (len(logw), len(theta), opts.mode))
    if logw is not None:
        idx, note = fairdraw_indices(logw, rng)
        if idx is not None:
            theta, lnL = theta[idx], np.asarray(lnL)[idx]
        elif note.startswith("FAILED"):
            # NO ARTIFACT WHEN THE FAIR DRAW FAILS.  Writing the raw cloud here
            # put proposal/prior draws into the standard, weightless
            # `*_samples.dat` product under the ONE name every consumer reads as
            # equal-weight posterior draws -- and none of them is obliged to read
            # the provenance line -- so the contract broke precisely on the
            # collapsed integrations where the difference is largest.  Refuse to
            # produce the file at all: the event fails loudly (and
            # --soft-fail-event-range still skips to the next one).
            #
            # The `.dat` result row goes too.  analyze_one validates the export
            # BEFORE write_dat, so this run has not written one -- but a file
            # left at that path by an earlier run would be collected as this
            # event's successful integration, which is the same stale-artifact
            # failure as for the samples file.
            _remove_stale_artifact(sname, "export")
            _remove_stale_artifact(dat_path(opts, out_index), "result")
            raise RuntimeError(
                "fair draw failed for output index %d (%s): the exported cloud "
                "would be the raw, UNREWEIGHTED sampler/proposal samples, which "
                "is not a posterior sample.  No %s written.  This integration "
                "collapsed -- fix the run (more samples, better proposal, "
                "sanity-check the likelihood) rather than reading the cloud."
                % (out_index, note, sname))
    # Drop non-finite lnL FIRST.  Doing it last meant the count was applied to
    # rows that were then discarded, and the provenance n_out counted them: a
    # header saying n_out=1000 above a 963-row file, and 137 above 129 -- wrong
    # exactly when the likelihood misbehaved, which is when someone reads it.
    lnL = np.asarray(lnL)
    good = np.isfinite(lnL)
    if not good.all():
        theta, lnL = theta[good], lnL[good]
    # THE count contract, applied once.  Gated on _FAIRDRAW_MODES -- the SAME set
    # check_critical_and_report uses to decide whether to report these flags as
    # ignored, so report and behaviour cannot disagree.  Applying it everywhere
    # meant `--mode nuts --fairdraw-extrinsic-output` printed "IGNORED" and then
    # silently wrote 5 rows instead of 300; --fairdraw-extrinsic-output is in
    # ILE_extr.sub, so that is a real production command line losing 60x of its
    # export under a banner saying the flag did nothing.
    # Rows are equal weight by this point, so subsample WITHOUT replacement -- a
    # random subset of an equal-weight cloud is still a fair draw (verified) and
    # manufactures no duplicates.
    if opts.mode in _FAIRDRAW_MODES:
        # NO evidence-neff argument, deliberately: this function has no parameter
        # that could carry one.  See fairdraw_size for why an ESS term here was
        # wrong on the flowMC modes.
        n_req = fairdraw_size(opts, len(theta))
        if n_req is not None and n_req < len(theta):
            n_before = len(theta)
            sub = rng.choice(n_before, size=int(n_req), replace=False)
            theta, lnL = theta[sub], lnL[sub]
            print("  fairdraw: exporting %d of %d rows (requested count)"
                  % (int(n_req), n_before))
    provenance = "fairdraw: %s n_out=%d" % (note, len(theta))
    if angle_note:
        provenance += " " + angle_note
    phase_marginalized = bool(getattr(opts, "phase_marginalization", False))
    if phase_marginalized:
        provenance += " phase=analytically-marginalized"

    ndim = theta.shape[1] if theta.ndim == 2 else len(theta)
    if with_distance:
        # 6-D: ra, dec, psi, incl, phiref, dist
        if phase_marginalized:
            cols = np.column_stack([theta[:, 0], theta[:, 1], theta[:, 5],
                                    theta[:, 3], theta[:, 2], lnL])
            hdr = "right_ascension declination distance inclination psi loglikelihood"
        else:
            cols = np.column_stack([theta[:, 0], theta[:, 1], theta[:, 5],
                                    theta[:, 3], theta[:, 2], theta[:, 4], lnL])
            hdr = "right_ascension declination distance inclination psi phi_orb loglikelihood"
    elif ndim == 4 and opts.mode == "flowmc-dpsimarg":
        # 4-D (flowmc-dpsimarg): theta = ra, dec, phiref, incl (psi marginalised,
        # phi_ref sampled).  Write ra, dec, incl, phi_orb.
        cols = np.column_stack([theta[:, 0], theta[:, 1],
                                theta[:, 3], theta[:, 2], lnL])
        hdr = "right_ascension declination inclination phi_orb loglikelihood"
    elif ndim == 4:
        # 4-D (flowmc-phimarg): ra, dec, psi, incl (phi_ref marginalised out)
        cols = np.column_stack([theta[:, 0], theta[:, 1],
                                theta[:, 3], theta[:, 2], lnL])
        hdr = "right_ascension declination inclination psi loglikelihood"
    elif ndim == 3:
        # 3-D (flowmc-phipsimarg): ra, dec, incl (phi_ref AND psi marginalised out)
        cols = np.column_stack([theta[:, 0], theta[:, 1],
                                theta[:, 2], lnL])
        hdr = "right_ascension declination inclination loglikelihood"
    else:
        # 5-D: ra, dec, psi, incl, phiref
        if phase_marginalized:
            cols = np.column_stack([theta[:, 0], theta[:, 1], theta[:, 3],
                                    theta[:, 2], lnL])
            hdr = "right_ascension declination inclination psi loglikelihood"
        else:
            cols = np.column_stack([theta[:, 0], theta[:, 1], theta[:, 3],
                                    theta[:, 2], theta[:, 4], lnL])
            hdr = "right_ascension declination inclination psi phi_orb loglikelihood"
    # Column line FIRST (unchanged, so `head -1` parsers keep working); the
    # provenance line follows, so the artifact records how it was produced --
    # notably the export ESS, which was previously written nowhere.
    np.savetxt(sname, cols, header=hdr + "\nmode=%s %s" % (opts.mode, provenance))
    print("Wrote %s (%d samples)" % (sname, len(cols)))


def _parse_freqresponse_arm_length(value):
    if value in (None, ""):
        return None
    text = str(value)
    if "=" not in text:
        return float(text)
    result = {}
    for item in text.split(","):
        fields = item.split("=", 1)
        if len(fields) != 2 or not fields[0].strip():
            raise ValueError("invalid --freqresponse-arm-length item %r" % item)
        result[fields[0].strip()] = float(fields[1])
    return result


def _waveform_precompute_kwargs(opts):
    """Return the waveform controls shared with production numpy ILE.

    These options act while constructing the mode time series, before JAX sees
    the packed likelihood data.  Omitting them therefore changes the numerical
    likelihood rather than merely selecting a JAX execution detail.
    """
    e_freq = getattr(opts, "e_freq", None)
    extra = {"fd_alignment_postevent_time": 2,
             "e_freq": 1 if e_freq is None else int(e_freq)}
    if getattr(opts, "internal_waveform_fd_L_frame", False):
        extra["fd_L_frame"] = True
    if getattr(opts, "internal_waveform_fd_no_condition", False):
        extra["no_condition"] = True
    return dict(
        use_gwsignal=bool(getattr(opts, "use_gwsignal", False)),
        use_gwsignal_approx=(opts.approximant
                             if getattr(opts, "use_gwsignal", False) else None),
        ignore_threshold=getattr(opts, "internal_precompute_ignore_threshold", None),
        no_memory=bool(getattr(opts, "no_memory", False)),
        extra_waveform_kwargs=extra)


def analyze_one(opts, P, data_dict, psd_dict, analyticPSD_Q, fiducial_epoch,
                rng, out_index, event_id, flow_state=None):
    """Build the JAX likelihood for one intrinsic template and run --mode.

    Returns ``(value, flow_state)`` where ``flow_state`` is the trained-flow
    state to bootstrap the next event (``--mode flowmc`` only; ``None`` else).
    """
    # Per-EVENT state: a batch run analyzes several events in one process, and
    # an undersizing on event 0 must not label event 1.
    _anglemarg.reset_amp_failsafe()
    tq = opts.time_marginalization_quadrature
    if tq == "bandlimited":
        required, g0, gcert = bandlimited_storage_requirement(
            P.deltaT, opts.data_integration_window_half)
        if opts.internal_data_storage_window_half < required:
            print("  widening rholm storage for bandlimited time support: "
                  "%.6g -> %.6g s" %
                  (opts.internal_data_storage_window_half, required))
            opts.internal_data_storage_window_half = required
        print("  precompute time support: guard_initial=%d guard_certificate=%d "
              "samples storage_half=%.6g s" %
              (g0, gcert, opts.internal_data_storage_window_half))
    print("Building JAX likelihood (production precompute + banded pack)...")
    _order_active = any((opts.check_slowrot_pmax, opts.check_finite_size_Qmax,
                         opts.choose_slowrot_pmax, opts.choose_slowrot_Qmax))
    order_control = None
    if _order_active:
        order_control = dict(
            target_snr=float(opts.response_order_snr),
            lnL_tolerance=float(opts.response_order_lnL_tol),
            n_samples=int(opts.response_order_sky_samples),
            p_reference=int(opts.response_order_p_reference),
            q_reference=int(opts.response_order_Q_reference),
            max_bank_gib=float(opts.response_order_max_bank_gib),
            check_p=bool(opts.check_slowrot_pmax),
            check_q=bool(opts.check_finite_size_Qmax),
            choose_p=bool(opts.choose_slowrot_pmax),
            choose_q=bool(opts.choose_slowrot_Qmax))
    waveform_kw = _waveform_precompute_kwargs(opts)
    if opts.rotation_slow and opts.freqresponse:
        arm_length = _parse_freqresponse_arm_length(opts.freqresponse_arm_length)
        like_data, extras = build_rotating_freqresponse_data_from_precompute(
            P.copy(), data_dict, psd_dict, fiducial_epoch,
            opts.data_integration_window_half, opts.l_max, opts.fmax,
            t_window=opts.internal_data_storage_window_half,
            Qmax=opts.freqresponse_qmax, L_arm=arm_length,
            p_max=opts.rotation_p_max, analyticPSD_Q=analyticPSD_Q,
            verbose=opts.verbose, order_control=order_control, **waveform_kw)
    elif opts.rotation_slow:
        nh = int(opts.rotation_n_harmonics)
        like_data, extras = build_rotation_data_from_precompute(
            P.copy(), data_dict, psd_dict, fiducial_epoch,
            opts.data_integration_window_half, opts.l_max, opts.fmax,
            t_window=opts.internal_data_storage_window_half,
            harmonics=tuple(range(-nh, nh + 1)), p_max=opts.rotation_p_max,
            analyticPSD_Q=analyticPSD_Q, verbose=opts.verbose,
            order_control=order_control, **waveform_kw)
    elif opts.freqresponse:
        arm_length = _parse_freqresponse_arm_length(opts.freqresponse_arm_length)
        like_data, extras = build_freqresponse_data_from_precompute(
            P.copy(), data_dict, psd_dict, fiducial_epoch,
            opts.data_integration_window_half, opts.l_max, opts.fmax,
            t_window=opts.internal_data_storage_window_half,
            Qmax=opts.freqresponse_qmax, L_arm=arm_length,
            analyticPSD_Q=analyticPSD_Q, verbose=opts.verbose,
            order_control=order_control, **waveform_kw)
    else:
        like_data, extras = build_data_from_precompute(
            P.copy(), data_dict, psd_dict, fiducial_epoch,
            opts.internal_data_storage_window_half, opts.data_integration_window_half,
            opts.l_max, opts.fmax,
            analyticPSD_Q=analyticPSD_Q, verbose=opts.verbose,
            q_time_pregrid_factor=(
                1 if getattr(opts, "q_time_pregrid_factor", 1) is None
                else int(opts.q_time_pregrid_factor)), **waveform_kw)
    if extras.get('order_report') is not None:
        opts.rotation_p_max = int(extras['order_report']['final_p'])
        opts.freqresponse_qmax = int(extras['order_report']['final_Q'])
        print("  response order used: p_max=%d Qmax=%d" %
              (opts.rotation_p_max, opts.freqresponse_qmax))
    print("  feature:", getattr(like_data, "feature", None),
          " modes:", like_data.lms,
          " guessed SNR:", extras.get("guess_snr", "not estimated"))
    if getattr(like_data, "feature", None) is not None:
        _d0 = like_data.detectors[like_data.detector_names[0]]
        _A = int(_d0["Q_bank"].shape[0])
        _bank_bytes = sum(int(like_data.detectors[d]["Q_bank"].nbytes
                              + like_data.detectors[d]["U_bank"].nbytes
                              + like_data.detectors[d]["V_bank"].nbytes)
                          for d in like_data.detector_names)
        print("  [banded] basis=%d ordered-pairs=%d device-bank=%.1f MB "
              "(%d detectors)" % (_A, _A * _A, _bank_bytes / 2.0**20,
                                    len(like_data.detector_names)))
    if like_data.q_time_pregrid_factor != 1:
        _d0 = like_data.detectors[like_data.detector_names[0]]
        _bytes = sum(int(like_data.detectors[d]["Q"].nbytes)
                     for d in like_data.detector_names)
        print("  [q-pregrid] factor %d: Q sampled at deltaT/%d, %d -> %d samples, "
              "%.1f MB of Q on device (%d detectors); integration cadence "
              "unchanged at deltaT=%.6g s"
              % (like_data.q_time_pregrid_factor, like_data.q_time_pregrid_factor,
                 _d0["npts_full_coarse"], _d0["Q"].shape[0], _bytes/2.0**20,
                 len(like_data.detector_names), like_data.deltaT))

    with_distance = not opts.distance_marginalization
    # --limit-distance: (d_lo,d_hi) is what is SAMPLED / quadratured; the prior is
    # always normalized over (opts.d_min, opts.d_max), which is what keeps the
    # reported lnZ on the same scale as a full-range run.
    d_lo, d_hi = resolve_distance_limit(opts)
    if (d_lo, d_hi) != (float(opts.d_min), float(opts.d_max)):
        print("  [limit] restricting distance sampling/quadrature to [%.4f, %.4f] Mpc "
              "(prior normalization UNCHANGED over [%.4f, %.4f] Mpc)"
              % (d_lo, d_hi, opts.d_min, opts.d_max))
    # Fail closed on a flag that would otherwise be silently inert.  The
    # log-uniform distance quadrature is wired -- and validated -- only for the
    # dense (phi,psi) schemes of --mode flowmc-phipsimarg.
    if (getattr(opts, "distance_grid_scheme", "uniform") != "uniform"
            and getattr(opts, "distance_grid_points", None) is not None):
        raise SystemExit(
            "--distance-grid-points and --distance-grid-scheme loguniform both "
            "set the distance node count.  The log-uniform grid derives its "
            "count from the data (see --distance-grid-tol); pass one or the "
            "other, not both.")
    # Resolved into a LOCAL, never written back onto opts: this driver has a
    # documented history of an option written back on event 0 being read as
    # event 1's choice in a batch loop.
    n_dist_grid = (DISTANCE_GRID_POINTS_DEFAULT
                   if opts.distance_grid_points is None
                   else int(opts.distance_grid_points))
    if (getattr(opts, "distance_grid_scheme", "uniform") != "uniform"
            or getattr(opts, "distance_grid_tol", None) is not None):
        if opts.mode != "flowmc-phipsimarg":
            raise SystemExit(
                "--distance-grid-scheme/--distance-grid-tol apply only to "
                "--mode flowmc-phipsimarg (they are validated only for the "
                "dense angle-marginalization schemes); got --mode %s."
                % opts.mode)
    if opts.mode in ("flowmc-phimarg", "nuts-phimarg"):
        # phi_ref-marginalised: requires distance marginalisation (baked in);
        # produces a 4-D (ra, dec, psi, incl) posterior.
        if not opts.distance_marginalization:
            raise SystemExit("--mode %s requires --distance-marginalization."
                             % opts.mode)
        from RIFT.likelihood.jax_ile.wrapper import JAXDistPhiMargLikelihood
        nphi = getattr(opts, "n_phi", 32)
        print("Distance + phi_ref marginalization: ON (grid=%d, nphi=%d, d in [%g,%g] Mpc)"
              % (n_dist_grid, nphi, d_lo, d_hi))
        like = JAXDistPhiMargLikelihood(
            like_data, d_lo, d_hi,
            nphi=nphi, n_grid=n_dist_grid,
            interp=opts.interp, guess_snr=extras["guess_snr"],
            time_quadrature=tq, d_prior_range=(opts.d_min, opts.d_max))
        if getattr(like, "dist_grid_info", {}).get("mode") == "adaptive":
            gi = like.dist_grid_info
            print("  distance grid: ADAPTIVE  d_peak=%.3g Mpc  sigma_d=%.3g Mpc  npts=%d"
                  % (gi["d_peak"], gi["sigma_d"], gi["n"]))
        with_distance = False   # distance is already marginalised inside like
        dim = 4
    elif opts.mode == "flowmc-phipsimarg":
        # distance + phi_ref + psi marginalised: 3-D (ra, dec, incl) posterior.
        if not opts.distance_marginalization:
            raise SystemExit("--mode flowmc-phipsimarg requires --distance-marginalization.")
        from RIFT.likelihood.jax_ile.wrapper import JAXDistPhiPsiMargLikelihood
        nphi = getattr(opts, "n_phi", 32)
        npsi = getattr(opts, "n_psi", 16)
        angle_marg = getattr(opts, "angle_marg_scheme", ANGLE_MARG_DEFAULT)
        print("Distance + phi_ref + psi marginalization: ON (grid=%d, nphi=%d, npsi=%d, d in [%g,%g] Mpc)"
              % (n_dist_grid, nphi, npsi, d_lo, d_hi))
        dist_grid = getattr(opts, "distance_grid_scheme", "uniform")
        dist_tol = getattr(opts, "distance_grid_tol", None)
        if dist_tol is None:
            dist_tol = _jax_core_dist_tol_default()
        elif dist_grid == "uniform":
            raise SystemExit(
                "--distance-grid-tol only applies to --distance-grid-scheme "
                "loguniform; it would be silently inert here.")
        policy = getattr(opts, "direct_marginalization_policy",
                         DIRECT_MARG_POLICY_DEFAULT)
        policy_config = None
        if policy != "off":
            _f0 = int(opts.direct_marginalization_reserve_time_refine)
            _scheme = getattr(opts, "direct_marginalization_reserve_scheme",
                              "exact")
            if _scheme != "auto":
                _angular, _time_rule = _direct_marg_reserve_pair(_scheme)
                print("  direct-marginalization reserve scheme: %s (angular "
                      "kernel %s, time rule %s)" % (_scheme, _angular, _time_rule))
            policy_config = DirectMargPolicyConfig(
                time_guard=int(opts.direct_marginalization_time_guard),
                reserve_scheme=_scheme,
                reserve_peaklocal_escalations=int(
                    opts.direct_marginalization_peaklocal_escalations),
                reserve_time_refine=_f0,
                reserve_time_refine_max=int(
                    opts.direct_marginalization_reserve_time_refine_max),
                total_value_error_budget_nats=float(
                    opts.direct_marginalization_error_budget_nats),
                reserve_batch_rows=int(
                    opts.direct_marginalization_batch_rows),
                max_modes=int(opts.direct_marginalization_max_modes),
                enriched_max_modes=int(
                    opts.direct_marginalization_enriched_max_modes),
                base_oversample=int(
                    opts.direct_marginalization_base_oversample),
                enriched_oversample=int(
                    opts.direct_marginalization_enriched_oversample),
                base_max_starts=int(opts.direct_marginalization_max_starts),
                max_time_nodes=int(
                    opts.direct_marginalization_max_time_nodes),
                convergence_tol_nats=float(
                    opts.direct_marginalization_convergence_tol_nats),
                time_guard_tol_nats=float(
                    opts.direct_marginalization_time_guard_tol_nats))
        try:
            like = JAXDistPhiPsiMargLikelihood(
                like_data, d_lo, d_hi, nphi=nphi, npsi=npsi,
                n_grid=n_dist_grid, interp=opts.interp,
                guess_snr=extras["guess_snr"], angle_marg=angle_marg,
                time_quadrature=tq, d_prior_range=(opts.d_min, opts.d_max),
                dist_grid=dist_grid, dist_grid_tol=dist_tol,
                direct_marginalization_policy=policy,
                policy_config=policy_config,
                bounded_multipeak_config=(
                    bounded_multipeak_config_from_options(opts)
                    if angle_marg == "multipeak-jax" else None),
                bounded_multipeak_decline_action=opts.multipeak_jax_decline_action)
        except ValueError as e:
            if policy == "off":
                raise
            raise SystemExit("--direct-marginalization-policy %s refused: %s"
                             % (policy, e))
        # ALWAYS report the resolved scheme (requested may be 'auto'; this
        # pipeline has a documented history of silently-inert flags).
        print("  angle-marg scheme: %s (requested %s): %s"
              % (like.angle_marg_scheme, angle_marg,
                 "; ".join("%s=%s" % kv for kv in
                           sorted(like.angle_marg_info.items())
                           if kv[0] not in ("scheme", "requested"))))
        print("  direct-marginalization policy: %s%s"
              % (like.direct_marginalization_policy,
                 "" if like.policy_info is None else ": " + "; ".join(
                     "%s=%s" % kv for kv in sorted(like.policy_info.items()))))
        print("  distance grid: %s"
              % "; ".join("%s=%s" % kv
                          for kv in sorted(getattr(like, "dist_grid_info", {}).items())))

        # THE PAIR IS CHOSEN AND ANNOUNCED BEFORE ANY ROW IS EVALUATED.
        # RO, 2026-09-08: rely on analysis and the known physics to pick the
        # (local, reserve) pair, rather than try-then-decline-then-refine.  A
        # run whose branch can only be read from a ledger at the end is a run
        # that spent its time in a method nobody chose.
        if like.direct_marginalization_policy != "off":
            _requested = getattr(
                opts, "direct_marginalization_reserve_scheme",
                DirectMargPolicyConfig().reserve_scheme)
            # WHAT IS ON THE MENU IS A PROPERTY OF THIS DATA, not a constant.
            # Offering laplace unconditionally would let the analysis choose a
            # reserve whose premise is absent here, and the run would then be
            # refused (or worse, carried) for a reason that has nothing to do
            # with the signal.  Two conditions, both about the DISTANCE
            # integral the reserve's table is contracted against:
            #
            #  1. The per-sample adaptive Gauss-Hermite quadrature must be on.
            #     On the legacy static uniform grid laplace is measured to cost
            #     43.2 nats at rho 163 -- that is what the GRID costs, not the
            #     scheme, and it is why the reserve may not use it there.  The
            #     loguniform static grid is NOT admitted either: it is sized
            #     from the angle amplitude and may well be adequate, but no
            #     measurement exists, and "probably fine" is the thing this
            #     selector was built to stop saying.
            #  2. The A0 == 0 / B1 == 0 identity the adaptive node placement is
            #     DERIVED from must hold on this data, measured on concrete
            #     tables.  Same predicate, same probe direction, and the same
            #     one the angle scheme is gated by -- one definition, in
            #     anglemarg.gh_laplace_supported_for_data.
            _available = ["exact"]
            _lap_why = None
            if int(getattr(opts, "_distance_gh_nodes_resolved", 0)) <= 0:
                _lap_why = ("the per-sample adaptive distance quadrature is "
                            "off (--distance-gh-nodes 0 / JAX_ILE_DISTMARG_GH "
                            "unset), and on a static distance grid the laplace "
                            "reserve is measured at 43.2 nats of error at "
                            "rho 163")
            else:
                # The wrapper has already measured this whenever the ANGLE
                # scheme could use the placement, and stored it.  Reuse that
                # answer rather than measuring a second time: two measurements
                # of one property can disagree, and the one the run reports
                # should be the one the run is gated on.
                _am_info = getattr(like, "angle_marg_info", None) or {}
                if "gh_laplace_ok" in _am_info:
                    _lap_ok, _lap_info = _am_info["gh_laplace_ok"], _am_info
                else:
                    _lap_ok, _lap_info = (
                        _anglemarg.gh_laplace_supported_for_data(
                            like_data,
                            getattr(like, "interp", None)
                            or JAX_INTERP_DEFAULT))
                if not _lap_ok:
                    _lap_why = _lap_info.get("gh_laplace_reason",
                                             "identity absent")
            if _lap_why is None:
                _available.append("laplace")
            else:
                print("  reserve roster: laplace NOT offered -- %s" % _lap_why)
            _pair, _pair_info = predict_reserve_pair(
                like_data, extras.get("guess_snr"),
                reserve_time_refine_max=int(
                    like.policy_config.reserve_time_refine_max),
                crossover_amplitude=ANGLE_MARG_CROSSOVER_AMPLITUDE,
                max_time_nodes=int(like.policy_config.max_time_nodes),
                requested=_requested, available=tuple(_available))
            _line = format_reserve_pair(_pair, _pair_info)
            print("  " + _line)
            sys.stderr.write("NOTE integrate_likelihood_extrinsic_jax: %s\n"
                             % _line)
            if _pair is None:
                # REFUSE, do not fall back.  Falling back to whole-window
                # refinement is the failure this analysis exists to prevent:
                # it would carry the rows without anyone choosing it.
                raise SystemExit(
                    "--direct-marginalization-reserve-scheme auto: %s  No "
                    "implemented reserve is adequate for this signal, so the "
                    "run is refused rather than silently falling back.  Pass "
                    "an explicit --direct-marginalization-reserve-scheme to "
                    "override, and record that you did."
                    % _pair_info.get("reason", ""))
            # THE PAIR THE ANALYSIS CHOSE MUST BE THE PAIR THAT RUNS.
            # Announcing it and leaving policy_config.reserve_scheme == "auto"
            # would make the printed line a claim about nothing: the composite
            # reads the config, not this print.  validate_policy_config admits
            # "auto" ONLY because it is resolved here, so this write is what
            # makes that admission true (found by the #304 session, whose
            # composite dispatches on the string).
            if _pair not in RESERVE_SCHEME_EXECUTABLE:
                raise SystemExit(
                    "the analysis selected the %r reserve (%s), and the "
                    "composite dispatches only %r.  Refused rather than "
                    "running a different reserve under the selected one's "
                    "name.  Wiring %r is RIFT PR #304."
                    % (_pair, _pair_info.get("reason", ""),
                       RESERVE_SCHEME_EXECUTABLE, _pair))
            # The wrapper's jitted closures capture the config at
            # construction, so writing the attribute alone would leave the
            # composite running the string 'auto' (which reserve_pair refuses):
            # rebuild the likelihood with the resolved pair.
            if _pair != like.policy_config.reserve_scheme:
                like = JAXDistPhiPsiMargLikelihood(
                    like_data, d_lo, d_hi, nphi=nphi, npsi=npsi,
                    n_grid=n_dist_grid, interp=opts.interp,
                    guess_snr=extras["guess_snr"], angle_marg=angle_marg,
                    time_quadrature=tq, d_prior_range=(opts.d_min, opts.d_max),
                    dist_grid=dist_grid, dist_grid_tol=dist_tol,
                    direct_marginalization_policy=policy,
                    policy_config=like.policy_config._replace(
                        reserve_scheme=_pair))
            # policy_info's "reserve_scheme" reported the resolved ANGLE
            # scheme, which is a different quantity that happened to share the
            # name while 'exact' was the only reserve.  Both are reported now,
            # under names that say which is which.
            like.policy_info["reserve_angle_scheme"] = \
                like.policy_info.get("reserve_scheme")
            like.policy_info["reserve_scheme"] = _pair
            print("  reserve pair RESOLVED: reserve_scheme=%s" % _pair)
        with_distance = False
        dim = 3
    elif opts.mode == "flowmc-dpsimarg":
        # distance + psi marginalised, phi_ref SAMPLED: 4-D (ra,dec,phiref,incl).
        if not opts.distance_marginalization:
            raise SystemExit("--mode flowmc-dpsimarg requires --distance-marginalization.")
        from RIFT.likelihood.jax_ile.wrapper import JAXDistPsiMargLikelihood
        npsi = getattr(opts, "n_psi", 8)
        print("Distance + psi marginalization (phi_ref sampled): ON "
              "(grid=%d, npsi=%d, d in [%g,%g] Mpc)"
              % (n_dist_grid, npsi, d_lo, d_hi))
        like = JAXDistPsiMargLikelihood(
            like_data, d_lo, d_hi, npsi=npsi,
            n_grid=n_dist_grid, interp=opts.interp,
            guess_snr=extras["guess_snr"], time_quadrature=tq,
            d_prior_range=(opts.d_min, opts.d_max))
        if getattr(like, "dist_grid_info", {}).get("mode") == "adaptive":
            gi = like.dist_grid_info
            print("  distance grid: ADAPTIVE  d_peak=%.3g Mpc  sigma_d=%.3g Mpc  npts=%d"
                  % (gi["d_peak"], gi["sigma_d"], gi["n"]))
        with_distance = False
        dim = 4
    elif opts.distance_marginalization:
        print("Distance marginalization: ON (grid=%d, d in [%g,%g] Mpc)"
              % (n_dist_grid, d_lo, d_hi))
        like = JAXDistanceMarginalizedLikelihood(
            like_data, d_lo, d_hi, n_grid=n_dist_grid,
            interp=opts.interp, phase_marginalization=opts.phase_marginalization,
            time_quadrature=tq, d_prior_range=(opts.d_min, opts.d_max))
        dim = 5
    else:
        like = JAXExtrinsicLikelihood(
            like_data, interp=opts.interp,
            phase_marginalization=opts.phase_marginalization,
            time_quadrature=tq)
        dim = 6

    if like.time_quadrature != tq:
        raise RuntimeError("constructed likelihood changed time quadrature: %r != %r"
                           % (like.time_quadrature, tq))
    print("  time-marginalization quadrature: %s (interpolate-time: %s)"
          % (like.time_quadrature, opts.interp))
    if like.time_quadrature == "bandlimited":
        print("  time guard: derived_initial=%d samples certified=%d samples; "
              "fine factor is curvature-derived per row with one doubling certificate"
              % (like.time_guard_initial, like.time_guard_certified))

    # In-flight branch classification, BEFORE any sampling.
    #
    # direct_marginalization_policy_note() below reports the same counts, but it
    # runs on the EXPORTED rows, so only a run that reaches the export can be
    # classified -- and a declining run is precisely the one that does not get
    # there.  A declined row executes the exact reserve (three evaluations: the
    # refined rule, its half-refined check, and the guard comparison) and rows
    # run one at a time under jax.lax.map, so declining costs enough that the
    # 32-row prior pilot itself does not finish.  "Did the controller accept"
    # was therefore unanswerable until the end of a run that never ends.
    #
    # The probe draws from the same prior the pilot draws from and calls the
    # same ledger, so it predicts the pilot rather than describing something
    # else.  Its compile is not extra work for a policy run: the end-of-run note
    # jits the identical _batched_ledger.
    _probe_rows = int(getattr(opts, "direct_marginalization_policy_probe_rows", 0))
    if (getattr(like, "direct_marginalization_policy", "off") != "off"
            and _probe_rows > 0):
        import time as _time
        from RIFT.likelihood.jax_ile.samplers import sample_prior_3 as _sp3
        # sample_prior_3 is correct here ONLY because the policy is refused
        # outside --mode flowmc-phipsimarg, whose parameter order is
        # (ra, dec, incl).  Check rather than assume: another order would hand
        # the ledger rows of the wrong width.
        _order = tuple(getattr(like, "ANGULAR_PARAM_ORDER", ()))
        if len(_order) != 3:
            raise RuntimeError(
                "policy probe expects the 3-parameter (ra, dec, incl) "
                "likelihood the policy is restricted to; this one is %r"
                % (_order,))
        _theta_probe = _sp3(_probe_rows, np.random.default_rng(opts.seed))
        _t0 = _time.time()
        _probe_note, _probe_vals, _probe_acc = (
            direct_marginalization_policy_note(
                like, _theta_probe, n_max=_probe_rows, chunk=_probe_rows,
                return_values=True))
        _dt = _time.time() - _t0
        # stderr AND stdout: the surveys read one or the other, and an
        # unclassifiable run is the failure this whole probe exists to remove.
        _msg = ("PROBE integrate_likelihood_extrinsic_jax: prior-draw policy "
                "ledger on %d rows in %.1f s (%.2f s/row incl. compile): %s"
                % (_probe_rows, _dt, _dt / max(_probe_rows, 1), _probe_note))
        if _probe_vals is not None:
            # lnL PAIRED WITH THE DISPOSITION, per row.  The counts alone
            # cannot say whether the rows that decline are rows that matter:
            # capacity declines concentrate far from truth, where the
            # likelihood is negligible, and a decline fraction quoted without
            # the likelihood of the declining rows overstates what it costs
            # the posterior.
            _msg += ("\nPROBE-VALUES seed=%d rows=%d lnL=[%s] local=[%s]"
                     % (int(opts.seed), _probe_rows,
                        ",".join("%.10g" % v for v in _probe_vals),
                        ",".join("1" if a else "0" for a in _probe_acc)))
        print(_msg)
        sys.stderr.write(_msg + "\n")
        sys.stdout.flush()
        sys.stderr.flush()
        if getattr(opts, "direct_marginalization_policy_probe_only", False):
            print("--direct-marginalization-policy-probe-only: exiting before "
                  "sampling; no samples and no evidence were written.")
            return None, None

    if opts.mode == "map":
        theta_map, lnL_map = run_map(like, opts, rng, dim, with_distance)
        fish = like.fisher(theta_map)
        names = FULL_NAMES if with_distance else ANG_NAMES
        print("peak lnL = %.5f" % lnL_map)
        for nm, v in zip(names, theta_map):
            print("   %-8s = %.6f" % (nm, v))
        print("Fisher diag:", np.array2string(np.diag(fish), precision=4))
        write_dat(opts, P, out_index, event_id, lnL_map, float("nan"), 0, float("nan"))
        return lnL_map, None

    out_flow_state = None
    if getattr(opts, "sampler_method", None) in ("AV", "portfolio"):
        from RIFT.likelihood.jax_ile import samplers as _samplers
        _portfolio_members = [
            name.strip() for item in (opts.sampler_portfolio or [])
            for name in str(item).split(",") if name.strip()]
        if not _portfolio_members:
            _portfolio_members = ["AV", "GMM"]
        _sample_bounds = resolve_av_angular_limits(opts)
        print("  JAX-%s: value-only integration, coverage chunk=%d, eval chunk=%s%s"
              % (opts.sampler_method, opts.n_chunk,
                 (str(opts.jax_av_eval_chunk)
                  if opts.jax_av_eval_chunk is not None else "auto"),
                 ("; members=" + ",".join(_portfolio_members))
                 if opts.sampler_method == "portfolio" else ""))
        _distance_sampling = av_distance_sampling_kwargs(like, d_lo, d_hi)
        res = _samplers.adaptive_volume_sample(
            like, opts.d_min, opts.d_max,
            sampler_method=opts.sampler_method,
            portfolio_members=_portfolio_members,
            nmax=opts.n_max, neff=(opts.n_eff or 1000),
            n_chunk=opts.n_chunk, eval_chunk=opts.jax_av_eval_chunk,
            seed=opts.seed, seed_method=opts.jax_av_seed,
            seed_pilot=opts.jax_av_seed_pilot,
            seed_modes=opts.jax_av_seed_modes,
            seed_points=opts.jax_av_seed_points,
            sky_inflate=opts.jax_av_sky_inflate,
            seed_prior_frac=opts.jax_av_seed_prior_frac,
            anisotropic_bins=opts.sampler_anisotropic_bins,
            verbose=opts.verbose, sample_bounds=_sample_bounds,
            distance_prior=opts.d_prior,
            **_distance_sampling)
        theta, lnL = res["theta"], res["lnL"]
        logZ, sig, neff = res["logZ"], res["sigma_over_Z"], res["neff"]
        ntot = res["n_eval"]
        logw_export = res.get("log_weight")
        print("  JAX-%s resolved fixed eval chunk=%d"
              % (opts.sampler_method, res["eval_chunk"]))
    elif opts.mode in ("multistart-nuts", "flowmc", "flowmc-phimarg",
                     "flowmc-phipsimarg", "flowmc-dpsimarg", "nuts-phimarg"):
        if opts.mode not in ("flowmc-phimarg", "flowmc-phipsimarg",
                             "flowmc-dpsimarg", "nuts-phimarg") and with_distance:
            raise SystemExit("--mode %s requires --distance-marginalization "
                             "(it samples the 5-D angular posterior)." % opts.mode)
        from RIFT.likelihood.jax_ile import samplers as _samplers
        n_starts = opts.num_chains if opts.num_chains and opts.num_chains > 1 else 8
        # Resolved per event and used directly.  Deliberately a LOCAL: writing it
        # back to opts made event 1 of a batch read event 0's choice.
        resolved_beta = float(opts.adapt_weight_exponent)
        if opts.mode in _TEMPERED_MODES:
            resolved_beta = resolve_tempering_exponent(
                opts, dim, tempered_cloud_size(opts, n_starts))
        if opts.mode == "nuts-phimarg":
            # Fisher-whitened multi-start NUTS on the 4-D phimarg posterior.
            res = _samplers.fisher_nuts_sample_phimarg(
                like, num_warmup=opts.num_warmup,
                num_samples=opts.num_samples,
                n_starts=max(n_starts, 8),
                n_prior_pilot=opts.n_prior_pilot,
                seed=opts.seed, verbose=opts.verbose)
        elif opts.mode == "multistart-nuts":
            res = _samplers.multistart_nuts(
                like, opts.d_min, opts.d_max, n_starts=n_starts,
                num_warmup=opts.num_warmup, num_samples=opts.num_samples,
                seed=opts.seed, sky_coords=opts.sky_coordinates,
                verbose=opts.verbose)
        elif (opts.mode in ("flowmc-phimarg", "flowmc-phipsimarg", "flowmc-dpsimarg")
              and getattr(opts, "smc_puffball", False)):
            # Robust high-SNR alternative to the flow: adaptive SMC with a
            # cloud-covariance "puffball" random-walk move (sample->puffball->sample,
            # a la RIFT-AV / nested sampling).  Does not collapse on sharp peaks.
            res = _samplers.smc_puffball_sample(
                like, opts.d_min, opts.d_max,
                n_walkers=opts.smc_walkers, n_move=opts.smc_move_steps,
                is_samples=opts.smc_is_samples,
                ess_frac=opts.temper_ess_frac, max_dbeta=opts.temper_max_dbeta,
                max_stages=max(opts.temper_max_stages, 80),
                puff_scale=opts.smc_puff_scale,
                seed=opts.seed, verbose=opts.verbose)
            # AN UNFINISHED LADDER IS NOT A POSTERIOR DRAW.  smc_puffball_sample
            # also stops at max_stages (and on a cloud that has lost its finite
            # lnL), leaving a cloud that targets L**inv_T * prior -- reachable
            # exactly in the sharp high-SNR regime this fallback exists for, where
            # the ESS rule picks very small temperature steps.  NEITHER artifact
            # may be published from that state: the export would ship a TEMPERED
            # cloud as an equal-weight posterior draw, and the `.dat` would carry
            # the SMC log Z(inv_T) rather than log Z.  Fail the event here, before
            # write_samples/write_dat, and clear anything an earlier run left at
            # those paths (--soft-fail-event-range still skips to the next event).
            _inv_T = res.get("inv_T")
            if _inv_T is None or not float(_inv_T) >= 1.0 - 1e-9:
                _remove_stale_artifact(samples_path(opts, out_index), "export")
                _remove_stale_artifact(dat_path(opts, out_index), "result")
                raise RuntimeError(
                    "--smc-puffball: the SMC temperature ladder reached inv_T=%s, "
                    "not 1, so the cloud still targets L**inv_T * prior and is "
                    "neither a posterior sample nor an evidence.  Nothing written "
                    "for output index %d.  Give the ladder room to finish (raise "
                    "--temper-max-stages / --smc-walkers, or lower "
                    "--temper-ess-frac) rather than exporting the tempered cloud."
                    % ("(not reported by the installed sampler)" if _inv_T is None
                       else "%.4g" % float(_inv_T), out_index))
        elif opts.mode in ("flowmc-phimarg", "flowmc-phipsimarg", "flowmc-dpsimarg"):
            # 4-D phi-marg (ra,dec,psi,incl), 3-D phi+psi-marg (ra,dec,incl), or
            # 4-D d+psi-marg (ra,dec,phiref,incl); flowmc_sample_phimarg is
            # dimension-agnostic (helpers chosen from ANGULAR_PARAM_ORDER).
            # Static tempering: beta = --adapt-weight-exponent -> temper = 1/beta.
            _beta = resolved_beta
            _temper = 1.0 / _beta if _beta > 0 else 1.0
            res = _samplers.flowmc_sample_phimarg(
                like, opts.d_min, opts.d_max,
                n_chains=max(n_starts, 20),
                n_local_steps=opts.n_local_steps,
                n_global_steps=opts.n_global_steps,
                n_training_loops=opts.n_training_loops,
                n_production_loops=opts.n_production_loops,
                n_epochs=opts.n_epochs,
                n_prior_pilot=opts.n_prior_pilot, temper=_temper,
                temper_adapt=opts.adapt_adapt, temper_init=opts.temper_init,
                temper_ess_frac=opts.temper_ess_frac,
                temper_max_stages=opts.temper_max_stages,
                temper_max_dbeta=opts.temper_max_dbeta,
                fisher_precondition=(opts.fisher_precondition
                                     or opts.fisher_is_samples > 0),
                fisher_is_samples=opts.fisher_is_samples,
                seed=opts.seed, reuse_state=flow_state, verbose=opts.verbose)
            out_flow_state = res.get("flow_state")
        else:
            _beta = resolved_beta
            _temper = 1.0 / _beta if _beta > 0 else 1.0
            res = _samplers.flowmc_sample(
                like, opts.d_min, opts.d_max, n_chains=max(n_starts, 20),
                n_local_steps=opts.n_local_steps,
                n_global_steps=opts.n_global_steps,
                n_training_loops=opts.n_training_loops,
                n_production_loops=opts.n_production_loops,
                n_epochs=opts.n_epochs,
                n_prior_pilot=opts.n_prior_pilot, temper=_temper,
                seed=opts.seed, reuse_state=flow_state, verbose=opts.verbose)
            out_flow_state = res.get("flow_state")  # bootstrap the next event
        theta, lnL = res["theta"], res["lnL"]
        logZ, sig, neff, ntot = res["logZ"], res["sigma_over_Z"], res["neff"], len(theta)
        # post_weight means DIFFERENT things per sampler, so only consume it
        # where it is an importance weight.
        #  * flowMC modes: L^(1-inv_T), the correction from the TEMPERED state
        #    actually sampled (exponent = --adapt-weight-exponent) to the exact
        #    posterior.  Uniform only at inv_T == 1.  This is a genuine w.
        #    Under --smc-puffball the same key carries the SMC ladder's
        #    L^(1-inv_T) correction, and the gate above has already refused the
        #    event unless the ladder finished -- so it is uniform by then.
        #  * multistart-nuts / nuts-phimarg: NOT an importance weight.
        #    samplers.py builds np.full(n_per[k], mass[k]/n_per[k]) -- a
        #    per-chain Laplace MODE-EVIDENCE weight, constant within a chain
        #    (multistart_nuts estimates it as peak_k + 0.5*logdet(sky cov),
        #    which its own comment flags as approximate).  Fair-drawing against
        #    it would re-weight whole chains by that estimator and duplicate
        #    ~half the rows, not correct any proposal mismatch.  It may well fix
        #    a real mode-mass bias, but that is unvalidated -- so these modes
        #    export their chains unreweighted until it is measured.
        _pw = res.get("post_weight") if opts.mode in _TEMPERED_MODES else None
        if _pw is not None and len(_pw) != len(theta):
            # Silently dropping the weights here made the export claim no
            # reweighting was needed; write_samples now refuses a mismatch, and
            # this raises at the source so the mode is named in the message.
            raise RuntimeError(
                "--mode %s returned %d post_weight entries for %d samples; the "
                "tempering correction cannot be applied."
                % (opts.mode, len(_pw), len(theta)))
        logw_export = (np.log(np.asarray(_pw, dtype=float))
                       if _pw is not None else None)
    elif opts.mode == "prior-mc":
        logZ, sig, neff, ntot, theta, lnL, logw_export = run_prior_mc(like, opts, rng, dim, with_distance)
    elif opts.mode == "nuts":
        logZ, sig, neff, ntot, theta, lnL, logw_export = run_nuts(like, opts, rng, with_distance)
    else:
        logZ, sig, neff, ntot, theta, lnL, logw_export = run_laplace_is(like, opts, rng, dim, with_distance)

    print("\n==== Result (event %d) ====" % event_id)
    print("  log evidence (lnL marginal over extrinsic) = %.5f" % logZ)
    print("  sigma_lnL = %.4g    neff = %.1f    ntotal = %d" % (sig, neff, ntot))
    require_finite_evidence(logZ, neff, opts.mode)
    # EXPORT FIRST, THEN PUBLISH THE RESULT ROW.  write_samples raises when the
    # cloud admits no fair draw, and that refusal means the integration itself
    # collapsed -- so the event must leave NO artifact behind.  Writing the
    # `.dat` first published a normal ILE result row (a finite, collapsed
    # evidence in the one-finite-weight case) that --soft-fail-event-range then
    # left in place for downstream collectors while the batch carried on.
    #
    # write_samples takes NO rng: it derives its own from (seed, out_index).
    # Passing the shared `rng` here -- which also feeds run_laplace_is /
    # run_prior_mc and the samplers -- made --save-samples, an OUTPUT flag,
    # change the lnL/logZ of every later event in the batch.
    # Report ONCE per event, here, BEFORE either artifact is written -- and
    # independently of --save-samples/--output-file, since write_samples()
    # early-returns when export is off and a run with export disabled would
    # otherwise get no warning and no persistent label at all.
    # Compute ONCE from the RESOLVED scheme and hand the same string to both
    # writers.  Recomputing inside each writer with no argument left `scheme`
    # None, so the standing output-cloud label never emitted and every artifact
    # stayed silent -- an inert guard, which is the exact failure mode this
    # label exists to prevent.
    _scheme = getattr(like, "angle_marg_scheme", None)
    require_bounded_multipeak_rows(like, lnL)
    _ev_note = angle_grid_suspect_note(_scheme)
    if _scheme == "multipeak-jax":
        _ev_note += " config=" + json.dumps(
            like.angle_marg_info["config"], sort_keys=True, separators=(",", ":"))
        _action = like.bounded_multipeak_decline_action
        _ev_note += " decline_action=" + _action
        if _action == "drop":
            _ev_note += " evidence_scope=accepted-region omitted_mass=unbounded"
            # Evidence has already counted every proposal, including finite
            # log-zero declines. Only the output cloud is filtered here.
            _keep = np.isfinite(lnL) & (np.asarray(lnL) != BOUNDED_MULTIPEAK_LOG_ZERO)
            _ev_note += " output_rows_dropped=%d" % int(np.sum(~_keep))
            theta, lnL = np.asarray(theta)[_keep], np.asarray(lnL)[_keep]
            if logw_export is not None:
                logw_export = np.asarray(logw_export)[_keep]
        _ev_note += (" audit_scope=log_likelihood-batches"
                     " audit_excludes=scalar-MAP-Fisher-MALA"
                     " refuse_latch_scope=log_likelihood-batches"
                     " scalar_refuse_decline=nan-without-raise-or-latch")
        _ev_note += " host_evaluation_audit=" + json.dumps(
            like.bounded_multipeak_audit, sort_keys=True, separators=(",", ":"))
    # Record the resolved distance-GH-nodes count in the same artifact header
    # line as the angle-marg/policy notes and (via write_samples' provenance
    # line) the mode/ESS record, so a reader of either artifact can see
    # whether the per-sample quadrature ran without re-deriving it from the
    # command line or the environment.
    _gh_n = int(getattr(opts, "_distance_gh_nodes_resolved", 0) or 0)
    _ev_note = (_ev_note + " gh_nodes=%d" % _gh_n).strip()
    # Cross-axis policy audit: which branch the exported rows actually took.
    # Evaluated on a subsample AFTER sampling because the sampler consumes the
    # value-only path; the ledger is the same computation with its record kept.
    if (getattr(like, "direct_marginalization_policy", "off") != "off"
            and not np.all(np.isfinite(np.asarray(lnL)))):
        raise SystemExit(
            "--direct-marginalization-policy %s: %d of %d exported rows have a "
            "non-finite likelihood because the controller could not warrant "
            "them; refusing to publish samples or evidence.  See the ledger "
            "note: %s" % (like.direct_marginalization_policy,
                          int(np.sum(~np.isfinite(np.asarray(lnL)))),
                          int(np.size(lnL)),
                          direct_marginalization_policy_note(
                              like, theta, n_max=32)))
    _policy_note = direct_marginalization_policy_note(like, theta)
    if _policy_note:
        sys.stderr.write("NOTE integrate_likelihood_extrinsic_jax: %s\n"
                         % _policy_note)
        _ev_note = (_ev_note + " " + _policy_note).strip()
    if _ev_note.startswith("SUSPECT-ANGLE-GRID"):
        sys.stderr.write(
            "WARNING integrate_likelihood_extrinsic_jax: angle-marginalization "
            "dense grids were UNDERSIZED at some evaluated points this event "
            "(%s).  Samples and evidence are LABELLED SUSPECT in their headers; "
            "rebuild with a larger amp_sizing before using them.  The run was "
            "NOT aborted and no points were discarded -- discarding would excise "
            "exactly the region the estimator missed.\n" % _ev_note)
    elif _ev_note:
        # Deterministic pass over the artifact-producing cloud, with scope
        # stated in the label (training-only proposals are not claimed).
        sys.stderr.write(
            "NOTE integrate_likelihood_extrinsic_jax: %s\n" % _ev_note)
    write_samples(opts, out_index, theta, lnL, with_distance, angle_note=_ev_note,
                  logw=logw_export)
    write_dat(opts, P, out_index, event_id, logZ, sig, ntot, neff,
              angle_note=_ev_note)
    return logZ, out_flow_state


# ---------------------------------------------------------------------------
# Main
# ---------------------------------------------------------------------------
def main(argv=None):
    if _JAX_CACHE_DIR is not None:
        print("JAX persistent compilation cache:", _JAX_CACHE_DIR)
    optp = build_parser()
    argv = _normalize_interpolate_time_argv(argv)
    opts, _ = optp.parse_args(argv)
    # BEFORE anything reads an option: which tokens did the user actually type?
    record_supplied_options(opts, argv, optp)
    resolve_ile_interface_aliases(opts, optp)
    check_critical_and_report(opts, optp)

    if opts.event_time is None:
        if opts.inj_mode:
            opts.event_time = 1126259462.0
        else:
            optp.error("--event-time is required (frame mode)")
    fiducial_epoch = opts.event_time
    rng = np.random.default_rng(opts.seed)
    deltaT = _analysis_delta_t(opts)

    # --- data (loaded once; shared across the intrinsic batch) ---
    if opts.inj_mode:
        if opts.mass1 is None or opts.mass2 is None:
            optp.error("--inj-mode requires --mass1 and --mass2")
        P_inj, data_dict, psd_dict, detectors, analyticPSD_Q = load_injection(opts, fiducial_epoch)
        deltaF = data_dict[detectors[0]].deltaF
        P_list = [P_inj]            # injection mode is single-event
    else:
        if not opts.cache_file or not opts.channel_name:
            optp.error("frame mode requires --cache-file and --channel-name (or --inj-mode)")
        data_dict, psd_dict, detectors, analyticPSD_Q = load_frames(opts, fiducial_epoch)
        deltaF = data_dict[detectors[0]].deltaF
        P_list = load_templates(opts, fiducial_epoch, deltaF, deltaT)

    print("Detectors:", detectors)
    print("Analyzing %d event(s) [event index base = %d]" % (len(P_list), opts.event))

    # --- batch loop over intrinsic templates ---
    # For --mode flowmc the trained normalizing flow is threaded from one event
    # to the next (partial-flow re-use), unless --no-flow-reuse is given.
    flow_state = None
    for out_index, P in enumerate(P_list):
        P.deltaT = deltaT
        P.deltaF = deltaF
        event_id = opts.event if (opts.sim_xml or opts.sim_grid or opts.inj_mode) else -1
        print("\n======== event %d (output index %d) : m1=%.3f m2=%.3f ========"
              % (event_id, out_index, P.m1 / MSUN, P.m2 / MSUN))
        try:
            _, new_flow_state = analyze_one(
                opts, P, data_dict, psd_dict, analyticPSD_Q,
                fiducial_epoch, rng, out_index, event_id,
                flow_state=(None if opts.no_flow_reuse else flow_state))
            if not opts.no_flow_reuse and new_flow_state is not None:
                flow_state = new_flow_state
        except Exception as e:
            if getattr(opts, "soft_fail_event_range", False):
                print("  event %d failed (soft-fail): %s" % (event_id, e))
                continue
            raise


if __name__ == "__main__":
    main()
