cosmo_val: pair-count shape noise in the OneCovariance covariance

docs/pair-count-noise.md

@@ -0,0 +1,75 @@
+# Pair-count shape noise
+
+`cosmo_val` workflows (whole-sample and tomographic) use measured effective pairs for OneCovariance's real-space Gaussian shape-noise term by default.
+Set `covariance.shape_noise: analytic` to recover the uniform-area approximation.
+Other workflows retain their analytic default.
+This changes neither blinding nor the mixed, sample-variance, non-Gaussian or SSC terms.
+
+## Cache and normalization
+
+`covariance_pair_counts` reads only the configured RA, Dec, weight and (for tomography) bin-ID columns of the same selected shear catalogue as `xi`.
+It counts on each requested covariance grid, using TreeCorr NN correlations with weights `w` and `w²`, spherical coordinates, logarithmic bins and `min_top=6`, matching `xi`.
+The cache is independent of cosmology and contains no shapes, signals or blind records.
+The reader projects fields directly from the lazy catalogue chunks, including compound HDF5 datasets, rather than materialising rows containing shape fields.
+It is a separate rule so covariance-only reruns need no shape measurement or access to a blind.
+Its `threads`, memory and runtime can be overridden just like `xi` (large catalogues may need more than the default 40 GB).
+For the full v2 catalogue, use `--set-resources covariance_pair_counts:mem_mb=128000 covariance_pair_counts:runtime=1440`, alongside the run's existing resource overrides; the two NN passes can peak above the ξ rule's memory usage.
+
+For unordered auto pairs in bin b, the NPZ stores:
+
+- `weight`: W = Σ w_i w_j.
+- `weight_squared`: Q = Σ (w_i w_j)².
+- `npairs`: the raw unordered count, plus exact bin edges and nominal centres.
+- `ordered_effective_pairs`: 2 W²/Q.
+
+With independent isotropic shapes of constant per-component dispersion s, the pure-noise variance of either ξ component is `2 s⁴ Q/W²`.
+The total two-component RMS is σ_e = √2 s, so the same variance is `σ_e⁴/(2 N_eff)`, where `N_eff = W²/Q` is **unordered**.
+Mixing the two dispersion conventions gives a factor-of-four error.
+
+## OneCovariance interface
+
+Verified against OneCovariance `c8755d6de43aadaa1cd4fe132ecf1c1d534e61f8`, the source used by the v2 run with the `ddd58fe0` container.
+The external checkout is selected by `tools.onecov_executable`; the image supplies its dependencies, not the source.
+`onecov/cov_input.py` reads `[tabulated inputs files]` keys `npair_directory` and `npair_mm_file` (first column θ in arcmin, last column counts).
+`onecov/cov_setup.py:get_npair` accepts a matching grid without rebinning and otherwise falls back to an unsafe-for-effective-counts sum.
+The workflow therefore writes **nominal geometric bin centres**, not measured mean separations, and counts separately on every covariance grid.
+Empty or non-finite pair bins fail instead of silently restoring analytic noise.
+
+OneCovariance's analytic auto count is `A_arcmin² n_eff² π(θ_hi²−θ_lo²)`, which counts **ordered** pairs.
+At `onecov/cov_theta_space.py` lines 1288–1297, the two Kronecker contractions and the `/0.5` yield `4 s⁴/N_ordered` for a single source bin.
+The two-column `.txt` thus supplies `N_ordered = 2 W²/Q`, giving exactly `2 s⁴ Q/W²`.
+No area rescaling or post-run whole-covariance scaling is applied.
+
+The normalization test checks the uniform-field expectation analytically, including non-unit equal weights and both dispersion conventions.
+A four-object positions-only catalogue checks TreeCorr's unordered convention against explicit pair enumeration.
+HDF5 tests forbid whole-dataset or shape-column reads for single, chunked and joined layouts.
+DAG tests verify the cache-to-ini dependency, analytic opt-out, dispersion boundary and tomographic file order and footprint inputs.
+
+## Tomographic pairs
+
+Each positive tomographic bin pair is counted separately, including cross correlations between disjoint bins.
+`shear.tomo_bin_col` names the bin-ID field; zero and negative IDs are excluded, matching ξ± selection.
+One-based indices in filenames refer to the sorted positive catalogue IDs, in the same order as the n(z) columns and `cov_th.tomo` survey parameters.
+The workflow writes `{version}_pairs_minsep=..._maxsep=..._nbins=..._tomo_bin=i-j.npz` and `.txt`, ordered `(1,1), (1,2), …, (2,2), …` in `npair_mm_file`.
+Each cache records the bin indices and the actual catalogue IDs.
+The area's and mask's source remains `cov_th.tomo` as for analytic noise.
+Geometry references must also have identical bin IDs when tomography is enabled.
+
+For an auto pair the table contains **2W²/Q**, with `ordered_effective_pairs` in the cache.
+For a cross pair it contains **W²/Q**, with `effective_cross_pairs` in the cache: TreeCorr counts each cross pair once, and only one of OneCovariance's two Kronecker contractions survives.
+Both yield the independent-shape variance `2 s_i² s_j² Q/W²`.
+Passing 2W²/Q for a cross would underestimate its noise by a factor of two.
+Pure shape noise has no covariance between different unordered bin pairs or between ξ+ and ξ−.
+
+The opt-in test `test_pair_counts_onecov.py` executes the external checkout's actual file reader and real-space noise kernel with unequal bin dispersions and explicit weighted-pair answers.
+Only the Fourier cosmological terms are stubbed to zero; it does not run the full cosmological integration.
+Run it inside the dependency container with `ONECOVARIANCE_PATH` set to the workflow's external source checkout.
+The tomographic probe is verified against candide checkout `70265be3a4bffcfded7b9251a3c4cbe8159d793f`.
+
+## Existing products
+
+Existing covariance files do not encode the noise prescription in their filenames.
+When migrating a completed run, use a fresh output tree or explicitly rerun `covariance_ini` and its descendants; do not interpret an old cached matrix as measured-pair noise.
+Do not sum fine-grid effective counts to obtain coarse-grid counts: sum W and Q separately before forming W²/Q, or count directly on that grid as this rule does.
+Pair noise includes the observed spatial selection and weight variation; it is not an independent mask-only simulation.
+Mask-dependent mixed and sample-variance corrections remain separate validation work.

papers/cosmo_val/config/config.yaml

@@ -114,6 +114,9 @@ var_method: jackknife
 
 covariance:
   default_masked: true
+  # Real-space shape noise, including tomography: measured W²/Q includes pair loss.
+  # Set analytic to restore OneCovariance's uniform-area approximation.
+  shape_noise: pair_counts
 
 cl:
   n_ell_bins: 32

src/sp_validation/pair_counts.py

@@ -0,0 +1,147 @@
+"""Non-cosmological weighted pair counts for real-space shape noise.
+
+TreeCorr auto counts are unordered. OneCovariance's xi_pm noise is
+4*s**4/N_ordered, where s is the per-component dispersion, so supply
+N_ordered = 2*W**2/Q for autos. Cross bins use W**2/Q and one
+Kronecker contraction. No shape columns are accessed.
+"""
+
+from pathlib import Path
+
+import numpy as np
+import treecorr
+
+from sp_validation.io import Catalogue
+
+
+def ordered_effective_pairs(weight, weight_squared):
+    """Convert unordered W=sum(w_i*w_j), Q=sum((w_i*w_j)**2) to 2 W²/Q.
+
+    Empty/non-finite bins cannot define a covariance; fail rather than
+    silently falling back to analytic counts.
+    """
+    weight = np.asarray(weight, dtype=float)
+    weight_squared = np.asarray(weight_squared, dtype=float)
+    if (
+        weight.shape != weight_squared.shape
+        or not np.all(np.isfinite(weight))
+        or not np.all(np.isfinite(weight_squared))
+        or np.any(weight <= 0)
+        or np.any(weight_squared <= 0)
+    ):
+        raise ValueError("Pair counts require matching, finite, positive W and Q")
+    return 2 * weight**2 / weight_squared
+
+
+def measure_pair_counts(
+    catalogue,
+    output,
+    *,
+    min_sep,
+    max_sep,
+    nbins,
+    ra_col="RA",
+    dec_col="Dec",
+    w_col="w",
+    num_threads=1,
+    tomo_bin_col=None,
+    bin_pair=None,
+):
+    """Cache W, Q and raw counts using positions and weights only.
+
+    Uses the xi rule's log binning, spherical metric and min_top=6.
+    The adjacent two-column .txt has nominal log-bin centres in arcmin,
+    never mean pair separations: OneCovariance must not rebin W²/Q.
+    ``bin_pair`` contains one-based indices into the sorted positive IDs in
+    ``tomo_bin_col``, matching the n(z) column order used by OneCovariance.
+    """
+    # Catalogue.read materialises whole compound HDF5 chunks before selecting
+    # columns. Index the public lazy chunk tables field-by-field instead, so
+    # shapes are never loaded, for either FITS or any supported HDF5 layout.
+    columns = (ra_col, dec_col, w_col) + ((tomo_bin_col,) if bin_pair else ())
+    with Catalogue(str(catalogue)) as cat:
+        data = np.empty(len(cat), dtype=cat.dtype(columns))
+        start = 0
+        for table in cat.chunks.values():
+            end = start + len(table)
+            for name in columns:
+                data[name][start:end] = table[name]
+            start = end
+    if bin_pair:
+        ids = np.unique(data[tomo_bin_col])
+        ids = ids[ids > 0]
+        # OneCovariance indexes the sorted positive bins from one, as n(z) does.
+        selected_ids = ids[np.asarray(bin_pair) - 1]
+        selections = [data[tomo_bin_col] == value for value in selected_ids]
+    else:
+        selections = [slice(None)] * 2
+    auto = bin_pair is None or bin_pair[0] == bin_pair[1]
+    weights = np.asarray(data[w_col], dtype=float)
+    if (
+        not np.all(np.isfinite(weights))
+        or np.any(weights < 0)
+        or not np.all(np.isfinite(data[ra_col]))
+        or not np.all(np.isfinite(data[dec_col]))
+    ):
+        raise ValueError("Positions and non-negative weights must be finite")
+    config = dict(
+        min_sep=float(min_sep),
+        max_sep=float(max_sep),
+        nbins=int(nbins),
+        sep_units="arcmin",
+        min_top=6,
+        num_threads=int(num_threads),
+    )
+    constant_weight = (
+        float(weights[0]) if len(weights) and np.all(weights == weights[0]) else None
+    )
+    counts = []
+    for power in (1, 2):
+        cats = [
+            treecorr.Catalog(
+                ra=data[ra_col][mask],
+                dec=data[dec_col][mask],
+                w=weights[mask] ** power,
+                ra_units="deg",
+                dec_units="deg",
+            )
+            for mask in selections[: 1 if auto else 2]
+        ]
+        nn = treecorr.NNCorrelation(config)
+        nn.process(*cats)
+        counts.append(nn.weight.copy())
+        print(f"Completed NN weights^{power}", flush=True)
+        if constant_weight is not None:
+            # For w_i=c, W=c² N and Q=c⁴ N=c² W on the same geometry.
+            # The second traversal carries no additional information.
+            counts.append(counts[0] * constant_weight**2)
+            print("Derived squared-weight counts for constant weights", flush=True)
+            break
+    W, Q = counts
+    # One contraction survives for disjoint cross bins, two for an auto.
+    # OneCovariance therefore needs W²/Q for crosses and 2 W²/Q for autos.
+    effective = ordered_effective_pairs(W, Q) / (1 if auto else 2)
+    edges = np.geomspace(float(min_sep), float(max_sep), int(nbins) + 1)
+    theta = np.sqrt(edges[:-1] * edges[1:])
+    output = Path(output)
+    output.parent.mkdir(parents=True, exist_ok=True)
+    np.savez(
+        output,
+        **(
+            {"bin_pair": np.asarray(bin_pair), "bin_ids": selected_ids}
+            if bin_pair
+            else {}
+        ),
+        theta=theta,
+        edges=edges,
+        npairs=nn.npairs,
+        weight=W,
+        weight_squared=Q,
+        **{("ordered_effective_pairs" if auto else "effective_cross_pairs"): effective},
+    )
+    np.savetxt(
+        output.with_suffix(".txt"),
+        np.column_stack((theta, effective)),
+        header="theta_arcmin "
+        + ("N_ordered_effective" if auto else "N_cross_effective"),
+    )