Fix c-bias error estimation

cosmo_val/cat_config.yaml

@@ -1108,6 +1108,12 @@ SP_v1.4.6.3_uncal:
   psf:
     path: unions_shapepipe_psf_2024_v1.4.a.fits
     hdu: 1
+  star:
+    ra_col: RA
+    dec_col: Dec
+    e1_col: e1
+    e2_col: e2
+    path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits
   patch_number: 100
 SP_v1.4.6.3_uncal_w_iv:
   blind: none
@@ -1124,6 +1130,12 @@ SP_v1.4.6.3_uncal_w_iv:
   psf:
     path: unions_shapepipe_psf_2024_v1.4.a.fits
     hdu: 1
+  star:
+    ra_col: RA
+    dec_col: Dec
+    e1_col: e1
+    e2_col: e2
+    path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits
   patch_number: 100
 SP_v1.4.6.3_uncal_w_1:
   blind: none
@@ -1140,6 +1152,12 @@ SP_v1.4.6.3_uncal_w_1:
   psf:
     path: unions_shapepipe_psf_2024_v1.4.a.fits
     hdu: 1
+  star:
+    ra_col: RA
+    dec_col: Dec
+    e1_col: e1
+    e2_col: e2
+    path: /n17data/UNIONS/WL/v1.4.x/unions_shapepipe_star_2024_v1.4.a.fits
   patch_number: 100
 SP_v1.4.5_uncal:
   blind: none

src/sp_validation/cosmo_val/catalog_characterization.py

@@ -14,11 +14,14 @@
 import healpy as hp
 import matplotlib.pyplot as plt
 import numpy as np
+import treecorr
 from cs_util import plots as cs_plots
+from uncertainties import ufloat
 
 from ..io import open_entry
 from ..survey import (
     additive_bias,
+    additive_bias_errors,
     area_from_coords,
     effective_survey_stats,
     n_eff_density,
@@ -413,6 +416,105 @@ def calculate_additive_bias(self):
                 self._c1[ver], self._c2[ver] = additive_bias(e1, e2, w, R)
         self.print_done("Finished additive bias calculation.")
 
+    def calculate_additive_bias_errors(self, npatch=None):
+        """Shape-noise and jackknife errors of the additive bias.
+
+        The jackknife deletes one patch at a time; patches are the seeded
+        k-means patches of the 2PCF (``_patch_centers``), so it includes
+        spatially correlated contributions such as cosmic shear and PSF
+        residuals.
+
+        Parameters
+        ----------
+        npatch : int, optional
+            Number of jackknife patches; default ``self.npatch``. With
+            ``npatch <= 1`` only the shape-noise errors are computed.
+        """
+        npatch = int(npatch or self.npatch)
+        self.print_start(f"Calculating additive bias errors ({npatch} patches):")
+        self._c_err = {}
+        for ver in self.versions:
+            self.print_magenta(ver)
+            R = self.cc[ver]["shear"]["R"]
+            with self.results[ver].temporarily_read_data():
+                e1, e2, w = self._read_shear_cols(ver, "e1_col", "e2_col", "w_col")
+                patch = None
+                if npatch > 1:
+                    cols = {
+                        "ra": np.asarray(self.results[ver].dat_shear["RA"]),
+                        "dec": np.asarray(self.results[ver].dat_shear["Dec"]),
+                        "w": np.asarray(w),
+                    }
+                    patch = treecorr.Catalog(
+                        ra=cols["ra"],
+                        dec=cols["dec"],
+                        ra_units=self.treecorr_config["ra_units"],
+                        dec_units=self.treecorr_config["dec_units"],
+                        patch_centers=self._patch_centers(cols, npatch),
+                    ).patch
+                self._c_err[ver] = additive_bias_errors(e1, e2, w, R, patch=patch)
+        self.print_done("Finished additive bias errors.")
+
+    def print_additive_bias(self, npatch=None, out_base="c_non_tomographic", labels=None):
+        """Write the additive bias with its errors to text and LaTeX tables.
+
+        The text table ``<out_base>.txt`` lists c with shape-noise (sn) and
+        jackknife (jk) errors; the LaTeX table ``<out_base>.tex`` lists c in
+        units of 1e-4 with jackknife errors (shape-noise errors without
+        jackknife). LaTeX rows are labelled by ``labels``, else the config's
+        ``label``, else the version name.
+
+        Parameters
+        ----------
+        npatch : int, optional
+            Number of jackknife patches; see
+            :meth:`calculate_additive_bias_errors`.
+        out_base : str, optional
+            Output file name base in the output directory.
+        labels : dict, optional
+            LaTeX row label per version.
+        """
+        if not hasattr(self, "_c_err"):
+            self.calculate_additive_bias_errors(npatch=npatch)
+
+        def fmt(value, error, latex=False):
+            return f"${ufloat(value, error):.1uL}$" if latex else f"{ufloat(value, error):.1u}"
+
+        out_path = self._output_path(f"{out_base}.txt")
+        with open(out_path, "w") as f:
+            print(
+                f"{'# Version':30s} {'c_1 (sn)':25s} {'c_2 (sn)':25s}"
+                + f" {'c_1 (jk)':25s} {'c_2 (jk)':25s}",
+                file=f,
+            )
+            for ver in self.versions:
+                err = self._c_err[ver]
+                row = [fmt(self.c1[ver], err["sn"][0]), fmt(self.c2[ver], err["sn"][1])]
+                if err["jk"] is not None:
+                    row += [fmt(self.c1[ver], err["jk"][0]), fmt(self.c2[ver], err["jk"][1])]
+                print(f"{ver:30s} " + " ".join(f"{r:25s}" for r in row), file=f)
+        self.print_done(f"Additive bias table written to {out_path}")
+
+        lines = [
+            r"\begin{tabular}{lcc}",
+            r"\toprule",
+            r"version & $c_1 / 10^{-4}$ & $c_2 / 10^{-4}$ \\",
+            r"\midrule",
+        ]
+        for ver in self.versions:
+            err = self._c_err[ver]["jk"] or self._c_err[ver]["sn"]
+            label = (labels or {}).get(ver) or self.cc[ver].get("label", ver)
+            label = label.replace("_", r"\_")
+            lines.append(
+                f"{label} & {fmt(self.c1[ver] * 1e4, err[0] * 1e4, latex=True)}"
+                + f" & {fmt(self.c2[ver] * 1e4, err[1] * 1e4, latex=True)} \\\\"
+            )
+        lines += [r"\bottomrule", r"\end{tabular}"]
+        out_path = self._output_path(f"{out_base}.tex")
+        with open(out_path, "w") as f:
+            print("\n".join(lines), file=f)
+        self.print_done(f"Additive bias table written to {out_path}")
+
     @property
     def c1(self):
         if not hasattr(self, "_c1"):

src/sp_validation/cosmo_val/core.py

@@ -657,10 +657,17 @@ def _read_shear_cols(self, ver, *keys):
         ``dat_shear`` directly.
 
         Returns one array per key (a bare array, not a 1-tuple, when a single
-        key is requested).
+        key is requested). A column name ``one`` gives unit values, e.g.
+        ``w_col: one`` for unweighted statistics.
         """
+        dat = self.results[ver].dat_shear
         cols = tuple(
-            self.results[ver].dat_shear[self.cc[ver]["shear"][key]] for key in keys
+            (
+                np.ones(len(dat))
+                if self.cc[ver]["shear"][key] == "one"
+                else dat[self.cc[ver]["shear"][key]]
+            )
+            for key in keys
         )
         return cols[0] if len(cols) == 1 else cols
 

src/sp_validation/statistics.py

@@ -53,6 +53,38 @@ def jackknife_patch_centers(cat, npatch, seed=0, init="tree"):
     return centers
 
 
+def jackknife_weighted_mean(x, w, patch):
+    """Delete-one-patch jackknife of a weighted mean.
+
+    Parameters
+    ----------
+    x, w : array of float
+        Values and weights.
+    patch : array of int
+        Patch index of each value, ``>= 0``.
+
+    Returns
+    -------
+    float
+        Mean of the delete-one estimates over the non-empty patches.
+    float
+        Jackknife error ``sqrt((K - 1) / K sum_k (x_k - mean)^2)``, where
+        ``x_k`` is the weighted mean without patch ``k`` and ``K`` the number of
+        non-empty patches. Unlike resampling single objects, this includes
+        spatially correlated contributions to the error.
+    """
+    x = np.asarray(x, dtype=float)
+    w = np.asarray(w, dtype=float)
+    sw = np.bincount(patch, weights=w)
+    swx = np.bincount(patch, weights=w * x)
+    has_data = sw > 0
+    x_k = (swx.sum() - swx[has_data]) / (sw.sum() - sw[has_data])
+    n_k = len(x_k)
+    mean = np.mean(x_k)
+    err = np.sqrt((n_k - 1) / n_k * np.sum((x_k - mean) ** 2))
+    return float(mean), float(err)
+
+
 def jackknif_weighted_average2(
     data,
     weights,

src/sp_validation/survey.py

@@ -13,6 +13,8 @@
 import healpy as hp
 import numpy as np
 
+from .statistics import jackknife_weighted_mean
+
 
 def get_area(dd, area_tile, verbose=False):
     """Get area.
@@ -295,6 +297,42 @@ def additive_bias(e1, e2, w, R):
     return c1, c2
 
 
+def additive_bias_errors(e1, e2, w, R, patch=None):
+    """Errors of the additive-bias estimates of :func:`additive_bias`.
+
+    Parameters
+    ----------
+    e1, e2 : array of float
+        Ellipticity components (before response correction).
+    w : array of float
+        Per-galaxy weights.
+    R : float
+        Multiplicative shear response.
+    patch : array of int, optional
+        Jackknife patch index of each galaxy; if ``None``, no jackknife error.
+
+    Returns
+    -------
+    dict
+        ``"sn"``: shape-noise errors ``(dc1, dc2)``,
+        ``sqrt(sum w^2 (e/R - c)^2) / sum w``, which is ``sigma_e / sqrt(n)``
+        for unit weights; ``"jk"``: delete-one-patch jackknife errors
+        ``(dc1, dc2)``, or ``None`` without ``patch``.
+    """
+    w = np.asarray(w, dtype=float)
+    errors = {"sn": [], "jk": [] if patch is not None else None}
+    for e in (e1, e2):
+        e = np.asarray(e, dtype=float) / R
+        c = np.average(e, weights=w)
+        errors["sn"].append(float(np.sqrt(np.sum(w**2 * (e - c) ** 2)) / np.sum(w)))
+        if patch is not None:
+            errors["jk"].append(jackknife_weighted_mean(e, w, patch)[1])
+    errors["sn"] = tuple(errors["sn"])
+    if patch is not None:
+        errors["jk"] = tuple(errors["jk"])
+    return errors
+
+
 def effective_survey_stats(e1, e2, w, area_deg2):
     """Effective survey statistics for a shear catalogue.
 

src/sp_validation/tests/test_additive_bias_table.py

@@ -0,0 +1,48 @@
+"""The additive-bias table of CosmologyValidation, on a stub without catalogues."""
+
+import os
+
+from sp_validation.cosmo_val.catalog_characterization import (
+    CatalogCharacterizationMixin,
+)
+
+
+class _Stub(CatalogCharacterizationMixin):
+    def __init__(self, output_dir):
+        self.versions = ["A_v1", "B"]
+        self.cc = {"A_v1": {}, "B": {"label": "Bee"}, "paths": {"output": output_dir}}
+        self._c1 = {"A_v1": -1.6e-4, "B": 1e-4}
+        self._c2 = {"A_v1": 2.3e-4, "B": 2e-4}
+        self._c_err = {
+            "A_v1": {"sn": (3e-5, 3e-5), "jk": (8e-5, 6e-5)},
+            "B": {"sn": (3e-5, 3e-5), "jk": None},
+        }
+
+    def _output_path(self, *parts):
+        return os.path.join(self.cc["paths"]["output"], *parts)
+
+    def print_done(self, msg):
+        pass
+
+
+def test_print_additive_bias_tables(tmp_path):
+    stub = _Stub(str(tmp_path))
+
+    stub.print_additive_bias(labels={"A_v1": "DES weights"})
+
+    txt = (tmp_path / "c_non_tomographic.txt").read_text().splitlines()
+    assert txt[1].split() == [
+        "A_v1",
+        "-0.00016+/-0.00003",
+        "0.00023+/-0.00003",
+        "-0.00016+/-0.00008",
+        "0.00023+/-0.00006",
+    ]
+    # Without jackknife only the shape-noise columns
+    assert len(txt[2].split()) == 3
+
+    tex = (tmp_path / "c_non_tomographic.tex").read_text()
+    # Jackknife errors where available, else shape noise; labels from the
+    # argument, then the config
+    assert r"DES weights & $-1.6 \pm 0.8$ & $2.3 \pm 0.6$ \\" in tex
+    assert r"Bee & $1.0 \pm 0.3$ & $2.0 \pm 0.3$ \\" in tex

src/sp_validation/tests/test_statistics.py

@@ -21,6 +21,7 @@
     cov_from_one_covariance,
     effective_number_of_tests,
     jackknif_weighted_average2,
+    jackknife_weighted_mean,
 )
 
 
@@ -335,3 +336,42 @@ def test_global_pte_is_calibrated_under_the_null():
     p = np.array([cal.global_pte(row)[0] for row in data])
     npt.assert_allclose(np.mean(p <= 0.05), 0.05, atol=0.012)
     assert 1.0 < cal.k_eff < 5.0
+
+
+def test_jackknife_weighted_mean_white_noise_matches_shot_noise():
+    """Uncorrelated values: the patch jackknife gives the shot-noise error."""
+    rng = np.random.default_rng(1)
+    n, npatch, sigma = 200_000, 50, 0.3
+    x = rng.normal(0.0, sigma, n)
+    w = rng.uniform(0.5, 1.5, n)
+    patch = rng.integers(0, npatch, n)
+
+    mean, err = jackknife_weighted_mean(x, w, patch)
+
+    shot = np.sqrt(np.sum(w**2 * (x - np.average(x, weights=w)) ** 2)) / w.sum()
+    npt.assert_allclose(mean, np.average(x, weights=w), atol=1e-5)
+    npt.assert_allclose(err, shot, rtol=0.25)
+
+
+def test_jackknife_weighted_mean_sees_patch_offsets():
+    """A constant offset per patch raises the error by its scatter / sqrt(K)."""
+    rng = np.random.default_rng(2)
+    n, npatch = 200_000, 50
+    patch = rng.integers(0, npatch, n)
+    w = np.ones(n)
+    offset = rng.normal(0.0, 0.01, npatch)
+
+    _, err = jackknife_weighted_mean(offset[patch], w, patch)
+
+    npt.assert_allclose(err, offset.std(ddof=1) / np.sqrt(npatch), rtol=0.05)
+
+
+def test_jackknife_weighted_mean_skips_empty_patches():
+    x = np.array([1.0, 2.0, 3.0, 4.0])
+    w = np.ones(4)
+    with_gap = np.array([0, 0, 2, 2])
+    without_gap = np.array([0, 0, 1, 1])
+    npt.assert_allclose(
+        jackknife_weighted_mean(x, w, with_gap),
+        jackknife_weighted_mean(x, w, without_gap),
+    )

src/sp_validation/tests/test_survey.py

@@ -53,3 +53,19 @@ def test_get_area(self):
             sorted(tile_IDs) == sorted(self._tile_IDs),
             msg=f"{tile_IDs}!={self._tile_IDs}",
         )
+
+
+def test_additive_bias_errors_shape_noise_and_jackknife():
+    rng = np.random.default_rng(3)
+    n, npatch, R = 100_000, 20, 0.7
+    e1 = rng.normal(0.0, 0.25, n)
+    e2 = rng.normal(0.0, 0.25, n)
+    w = np.ones(n)
+    patch = rng.integers(0, npatch, n)
+
+    errors = survey.additive_bias_errors(e1, e2, w, R, patch=patch)
+
+    # Unit weights: sigma_e / sqrt(n), not var / sqrt(n)
+    npt.assert_allclose(errors["sn"], (np.std(e1 / R), np.std(e2 / R)) / np.sqrt(n))
+    npt.assert_allclose(errors["jk"], errors["sn"], rtol=0.4)
+    assert survey.additive_bias_errors(e1, e2, w, R)["jk"] is None