From 9eb2b187874bd69e6c7379d4b515cb75ae1b1a9f Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Thu, 8 Oct 2026 14:33:01 +0200 Subject: [PATCH 1/2] cosmo_val: shape-noise and spatial jackknife errors on the additive bias - statistics.jackknife_weighted_mean: delete-one-patch jackknife of a weighted mean. - survey.additive_bias_errors: shape-noise error sqrt(sum w^2 (e-c)^2)/sum w and the patch jackknife error of c1, c2. - CosmologyValidation.calculate_additive_bias_errors: jackknife over the seeded k-means patches of the 2PCF, so the error includes spatially correlated contributions (cosmic shear, PSF residuals); print_additive_bias writes c with both errors (txt) and with jackknife errors in units of 1e-4 (tex). - _read_shear_cols: column name 'one' gives unit weights (w_col: one). Co-Authored-By: Claude Opus 5.5 --- .../cosmo_val/catalog_characterization.py | 102 ++++++++++++++++++ src/sp_validation/cosmo_val/core.py | 11 +- src/sp_validation/statistics.py | 32 ++++++ src/sp_validation/survey.py | 38 +++++++ .../tests/test_additive_bias_table.py | 48 +++++++++ src/sp_validation/tests/test_statistics.py | 40 +++++++ src/sp_validation/tests/test_survey.py | 16 +++ 7 files changed, 285 insertions(+), 2 deletions(-) create mode 100644 src/sp_validation/tests/test_additive_bias_table.py diff --git a/src/sp_validation/cosmo_val/catalog_characterization.py b/src/sp_validation/cosmo_val/catalog_characterization.py index e2df6d8e..d871cd9d 100644 --- a/src/sp_validation/cosmo_val/catalog_characterization.py +++ b/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 ``.txt`` lists c with shape-noise (sn) and + jackknife (jk) errors; the LaTeX table ``.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"): diff --git a/src/sp_validation/cosmo_val/core.py b/src/sp_validation/cosmo_val/core.py index 8159a1a0..55417fb1 100644 --- a/src/sp_validation/cosmo_val/core.py +++ b/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 diff --git a/src/sp_validation/statistics.py b/src/sp_validation/statistics.py index 17e05219..eb01263d 100644 --- a/src/sp_validation/statistics.py +++ b/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, diff --git a/src/sp_validation/survey.py b/src/sp_validation/survey.py index e9d61025..2822254a 100644 --- a/src/sp_validation/survey.py +++ b/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. diff --git a/src/sp_validation/tests/test_additive_bias_table.py b/src/sp_validation/tests/test_additive_bias_table.py new file mode 100644 index 00000000..68e598e3 --- /dev/null +++ b/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 diff --git a/src/sp_validation/tests/test_statistics.py b/src/sp_validation/tests/test_statistics.py index d9e08277..979e1fdd 100644 --- a/src/sp_validation/tests/test_statistics.py +++ b/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), + ) diff --git a/src/sp_validation/tests/test_survey.py b/src/sp_validation/tests/test_survey.py index e786c8f8..40ded7e3 100644 --- a/src/sp_validation/tests/test_survey.py +++ b/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 From 70f7694181f730c24756faed535c52488be533fc Mon Sep 17 00:00:00 2001 From: martinkilbinger Date: Thu, 8 Oct 2026 14:38:43 +0200 Subject: [PATCH 2/2] cat_config: star catalogue for the SP_v1.4.6.3_uncal entries CosmologyValidation reads a star section for every version; the three uncalibrated v1.4.6.3 entries (DES, inverse-variance and unit weights) had none and failed with KeyError 'star'. Co-Authored-By: Claude Opus 5.5 --- cosmo_val/cat_config.yaml | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/cosmo_val/cat_config.yaml b/cosmo_val/cat_config.yaml index 41eb3281..45155a43 100644 --- a/cosmo_val/cat_config.yaml +++ b/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