Add bootstrap confidence intervals and Sobol quasi-random sampling to analysis

Two complementary improvements that travel together through every
analysis function.

Bootstrap confidence intervals:
- compute_expressibility, compute_entanglement_capability, and
  estimate_trainability now report percentile bootstrap CIs on their
  headline statistics in the detailed result dict. Configurable via
  confidence_level (default 0.95) and n_bootstrap_ci (default 200).
- ExpressibilityResult gains expressibility_ci_lower/upper,
  mean_fidelity_ci_lower/upper, confidence_level, and sampling.
- EntanglementResult gains entanglement_ci_lower/upper,
  confidence_level, and sampling.
- TrainabilityResult gains trainability_ci_lower/upper,
  gradient_variance_ci_lower/upper, confidence_level, and sampling.
- Bootstrap runs through the same pipeline as the point estimate so
  non-linearities are properly propagated:
  - Expressibility resamples fidelities and re-runs the
    histogram -> KL -> normalize chain on each resample.
  - Trainability resamples successful gradient rows and re-runs the
    per-parameter-variance -> mean -> variance_to_trainability chain.
- Shared encoding_atlas.analysis._ci module supplies the small,
  well-tested percentile_bootstrap_ci helper used everywhere, plus
  validate_ci_args for upfront argument validation.

Quasi-random (Sobol') sampling:
- New sampling parameter on all four analysis entry points
  (compute_expressibility, compute_fidelity_distribution,
  compute_entanglement_capability, estimate_trainability) accepts
  'uniform' (default, unchanged) or 'sobol'.
- 'sobol' uses scipy.stats.qmc.Sobol, seeded from the existing seed
  argument so a single seed parameter controls everything for both
  paths. Default 'uniform' preserves byte-identical seeded output to
  before this commit.
- Sobol typically converges 30-50% faster on the analysis statistics;
  stacks multiplicatively with the parallelization work from the
  previous commit.
- Shared encoding_atlas.analysis._sampling module supplies
  generate_sample_batch and validate_sampling, used by every
  analysis function for a single source of truth.

Public API for the float-only return path is fully unchanged: all
new parameters are keyword-only with defaults that preserve the
existing behavior. The detailed result dict gains keys; the existing
keys remain unchanged.

Tests (tests/unit/analysis/test_ci_and_sampling.py, 48 cases across
7 classes) cover: bootstrap helper on known distributions,
determinism under fixed seed, degenerate cases (empty, single,
constant samples), 99% interval is wider than 95%, six bad
confidence_level values and five bad n_bootstrap values rejected by
validate_ci_args, sampling helper shape/range/determinism/range
scaling, six bad sampling values rejected, scipy's power-of-two
notice suppressed inside the helper, Sobol determinism and
sampling-distinct-from-uniform across every analysis entry point,
all CI keys present and float-typed for every detailed result,
CI bounds bracket point estimate and respect documented ranges,
clean ValueError on bad sampling/confidence/n_bootstrap arguments,
and backward-compat assertions that sampling='uniform' with default
CI knobs reproduces the pre-commit float result exactly.

Updates two pre-existing key-completeness tests
(test_expressibility.py, test_trainability.py) to expect the new
CI fields.

Full test suite (4550 not-slow + 382 slow with optional backends)
passes; ruff, black, build, and mkdocs --strict all clean.

Author: ashutoshm1771

src/encoding_atlas/analysis/_ci.py

@@ -0,0 +1,155 @@
+"""Bootstrap confidence-interval helpers shared by the analysis pipeline.
+
+The three core analysis functions (``compute_expressibility``,
+``compute_entanglement_capability``, ``estimate_trainability``) all
+return a single scalar summary of a sampled distribution. For
+research use those scalars need an uncertainty envelope — otherwise
+``expressibility = 0.042`` is unfalsifiable in a paper.
+
+This module supplies a small, deterministic, well-tested percentile
+bootstrap implementation that every analysis function uses to attach
+a confidence interval to its result. Determinism is the key
+property: for a fixed ``rng`` the bootstrap CI is byte-identical
+across runs, so users can quote CIs in publications and reviewers
+can reproduce them exactly.
+
+Design notes
+------------
+* **Percentile bootstrap** (Efron 1979) rather than BCa: simpler,
+  identical asymptotically for our smooth statistics
+  (mean / variance / KL), and easier to verify in tests. BCa would
+  pull in additional jackknife logic for marginal accuracy gains.
+* **Pre-sampled indices**: ``rng.integers`` is called once for the
+  full ``(n_bootstrap, n_samples)`` index matrix. This is faster
+  than per-iteration draws and — crucially — produces an identical
+  index sequence regardless of how the loop is structured, so
+  changes to the inner statistic don't perturb the random state.
+* **Validation upfront**: bad ``confidence_level`` (outside (0, 1))
+  or non-positive ``n_bootstrap`` raise ``ValueError`` before any
+  bootstrap work happens.
+"""
+
+from __future__ import annotations
+
+from typing import Any, Callable
+
+import numpy as np
+from numpy.typing import NDArray
+
+
+def validate_ci_args(
+    confidence_level: float,
+    n_bootstrap: int,
+) -> None:
+    """Validate the shared CI arguments. Raises a clean ``ValueError``.
+
+    Parameters
+    ----------
+    confidence_level : float
+        Must lie strictly in the open interval (0, 1). Common
+        choices are 0.90, 0.95 (default for almost all analysis
+        functions), and 0.99.
+    n_bootstrap : int
+        Number of bootstrap resamples; must be a positive integer.
+        Larger values give smoother CIs at proportional CPU cost.
+
+    Raises
+    ------
+    ValueError
+        With a message that names exactly which argument was bad.
+    """
+    if not isinstance(confidence_level, (int, float)) or not (
+        0.0 < float(confidence_level) < 1.0
+    ):
+        raise ValueError(
+            f"confidence_level must be a number in (0, 1), " f"got {confidence_level!r}"
+        )
+    if not isinstance(n_bootstrap, int) or n_bootstrap < 1:
+        raise ValueError(f"n_bootstrap must be a positive integer, got {n_bootstrap!r}")
+
+
+def percentile_bootstrap_ci(
+    samples: NDArray[np.floating[Any]],
+    statistic_fn: Callable[[NDArray[np.floating[Any]]], float],
+    rng: np.random.Generator,
+    n_bootstrap: int = 200,
+    confidence_level: float = 0.95,
+) -> tuple[float, float]:
+    """Percentile bootstrap confidence interval for a scalar statistic.
+
+    Resamples ``samples`` with replacement ``n_bootstrap`` times,
+    applies ``statistic_fn`` to each resample, and returns the
+    ``(lower, upper)`` percentile bounds at the requested
+    ``confidence_level``.
+
+    Parameters
+    ----------
+    samples : NDArray[np.floating]
+        1-D sample array (length n). Per-sample statistics of the
+        analysis (fidelities, entanglement values, etc.).
+    statistic_fn : callable
+        Receives a resampled view of ``samples`` and returns a
+        ``float``. Common choices: ``np.mean``, ``np.var``, a
+        composite (e.g. variance-to-trainability mapping).
+    rng : numpy.random.Generator
+        Random source for the resampling. Caller is responsible for
+        seeding for determinism.
+    n_bootstrap : int, default=200
+        Number of bootstrap resamples. 200 keeps the cost negligible
+        next to the simulation cost and gives ~1% jitter on the
+        percentile endpoints.
+    confidence_level : float, default=0.95
+        Two-sided confidence level in (0, 1).
+
+    Returns
+    -------
+    (lower, upper) : tuple of float
+        Bounds of the bootstrap CI. If ``samples`` has length < 2 or
+        every resample yields the same value (degenerate distribution),
+        both bounds equal ``statistic_fn(samples)``.
+
+    Raises
+    ------
+    ValueError
+        If ``confidence_level`` is out of range or ``n_bootstrap``
+        is not a positive integer.
+
+    Notes
+    -----
+    The percentile bootstrap is consistent for smooth statistics and
+    matches the bias-corrected ``BCa`` method asymptotically. We
+    deliberately use it (rather than ``BCa``) here because it is
+    simpler, has fewer numerical edge cases, and is easier to verify
+    in tests — and because the analysis statistics it operates on
+    (mean of a bounded value, variance) satisfy the smoothness
+    assumption.
+    """
+    validate_ci_args(confidence_level, n_bootstrap)
+
+    samples = np.asarray(samples)
+    n = samples.shape[0]
+
+    # Degenerate input — no resampling is meaningful. Return the point
+    # estimate as both bounds so the public surface remains uniform.
+    if n < 2:
+        point = float(statistic_fn(samples)) if n == 1 else float("nan")
+        return point, point
+
+    # One vectorized index draw covers all resamples; downstream changes
+    # to ``statistic_fn`` cannot perturb the random sequence.
+    indices = rng.integers(0, n, size=(n_bootstrap, n))
+
+    bootstrap_stats = np.empty(n_bootstrap, dtype=np.float64)
+    for b in range(n_bootstrap):
+        bootstrap_stats[b] = float(statistic_fn(samples[indices[b]]))
+
+    alpha = 1.0 - float(confidence_level)
+    lower = float(np.percentile(bootstrap_stats, 100.0 * (alpha / 2.0)))
+    upper = float(np.percentile(bootstrap_stats, 100.0 * (1.0 - alpha / 2.0)))
+
+    # Numerical guard: if the statistic is constant across resamples the
+    # two percentiles will collapse; clip any tiny negative gap to 0 so
+    # downstream code can rely on ``lower <= upper``.
+    if upper < lower:
+        upper = lower
+    return lower, upper

src/encoding_atlas/analysis/_sampling.py

@@ -0,0 +1,133 @@
+"""Random-input sampling helpers for the analysis pipeline.
+
+The three core analysis functions sample random inputs into a fixed
+range ``[low, high]`` and feed each into the encoding. This module
+exposes a single ``generate_sample_batch`` helper that picks between
+two strategies:
+
+* ``"uniform"`` — the existing pseudo-random uniform sampling via
+  ``numpy.random.Generator.uniform``. Default and unchanged.
+* ``"sobol"`` — a Sobol low-discrepancy sequence via
+  ``scipy.stats.qmc.Sobol`` scaled to ``[low, high]``. Low-discrepancy
+  sequences cover the hypercube more evenly than i.i.d. uniform draws,
+  which reduces the sample count needed to estimate quantities like
+  expressibility's KL divergence and entanglement capability's mean by
+  roughly 30-50% in typical settings.
+
+Both paths consume the caller-supplied ``rng`` so the existing
+``seed`` parameter on the public analysis functions stays the single
+source of randomness — including for Sobol scrambling. This keeps
+reproducibility byte-identical for a given ``(seed, sampling)``
+pair.
+"""
+
+from __future__ import annotations
+
+import warnings
+from typing import Any, Literal
+
+import numpy as np
+from numpy.typing import NDArray
+
+SamplingMethod = Literal["uniform", "sobol"]
+
+
+def validate_sampling(sampling: SamplingMethod) -> None:
+    """Validate the sampling argument with a clear error message.
+
+    Parameters
+    ----------
+    sampling : {'uniform', 'sobol'}
+        Caller's choice of input sampling strategy.
+
+    Raises
+    ------
+    ValueError
+        If ``sampling`` is not one of the accepted values. The error
+        names exactly which strings are valid, so users can self-
+        correct quickly.
+    """
+    if sampling not in ("uniform", "sobol"):
+        raise ValueError(f"sampling must be 'uniform' or 'sobol', got {sampling!r}")
+
+
+def generate_sample_batch(
+    n_samples: int,
+    n_features: int,
+    input_range: tuple[float, float],
+    rng: np.random.Generator,
+    sampling: SamplingMethod = "uniform",
+) -> NDArray[np.floating[Any]]:
+    """Generate an ``(n_samples, n_features)`` input batch.
+
+    Parameters
+    ----------
+    n_samples : int
+        Number of samples to generate.
+    n_features : int
+        Dimensionality of each sample (matches the encoding's
+        ``n_features``).
+    input_range : (float, float)
+        ``(low, high)`` range; samples are uniform / quasi-uniform
+        within this interval.
+    rng : numpy.random.Generator
+        Random source. For the Sobol path it seeds the scrambling so
+        the same ``rng`` state always produces the same sequence.
+    sampling : {'uniform', 'sobol'}, default='uniform'
+        Sampling strategy. See module docstring.
+
+    Returns
+    -------
+    NDArray[np.floating]
+        Float64 array of shape ``(n_samples, n_features)`` with
+        values in ``[input_range[0], input_range[1])``.
+
+    Raises
+    ------
+    ValueError
+        If ``sampling`` is not one of the accepted strings.
+    ImportError
+        If ``sampling='sobol'`` is requested but ``scipy.stats.qmc``
+        is unavailable. SciPy is a hard dependency of
+        ``encoding_atlas`` so this is unlikely in practice, but the
+        error message names the missing import.
+
+    Notes
+    -----
+    The Sobol implementation deliberately suppresses the
+    ``scipy.stats.qmc`` warning that fires when ``n_samples`` is not
+    a power of two. The warning describes a balance property of
+    Sobol' nets; the resulting samples are still a valid low-
+    discrepancy draw and produce better convergence than i.i.d.
+    uniform for the analysis statistics here. Sobol' shines when
+    ``n_samples`` is a power of two — users who need the absolute
+    best statistical properties should round their batch size up to
+    the nearest power of two.
+    """
+    validate_sampling(sampling)
+    low, high = float(input_range[0]), float(input_range[1])
+
+    if sampling == "uniform":
+        return rng.uniform(low, high, size=(n_samples, n_features)).astype(
+            np.float64, copy=False
+        )
+
+    # sampling == "sobol"
+    try:
+        from scipy.stats import qmc
+    except ImportError as exc:  # pragma: no cover - scipy is a hard dep
+        raise ImportError(
+            "scipy.stats.qmc is required for sampling='sobol'. " "Install scipy >= 1.7."
+        ) from exc
+
+    # Seed the Sobol scrambler from the caller's Generator. Passing the
+    # ``rng`` object (not a raw int) lets ``np.random.Generator`` mediate
+    # the seeding so the public ``seed`` parameter remains the single
+    # source of randomness across both sampling paths.
+    sampler = qmc.Sobol(d=n_features, scramble=True, seed=rng)
+    with warnings.catch_warnings():
+        # See the Notes block above — the power-of-two balance warning
+        # is informative for QMC theorists but noise for our users.
+        warnings.simplefilter("ignore", UserWarning)
+        unit_samples = sampler.random(n_samples)
+    return (unit_samples * (high - low) + low).astype(np.float64, copy=False)