PSF leakage: joint errors-in-variables ρ/τ fit by default

papers/cosmo_val/config/config.yaml

@@ -51,7 +51,19 @@ cosmo_val:
   noise_bias_method: analytic
   cell_seed: 8192
   path_onecovariance: "/home/guerrini/OneCovariance/"
-  rho_tau_method: lsq
+  # PSF-leakage fit: 'eiv' is the joint errors-in-variables fit of the rho/tau
+  # jackknife measurement below; 'lsq'/'emcee' fit the published tau statistics
+  # with the cov_estimate_method covariance.
+  rho_tau_method: eiv
+  psf_leakage:
+    # k-means patches shared by stars and galaxies in the joint measurement
+    npatch_fine: 800
+    # jackknife patch counts the fit reports (Ward unions of the fine patches)
+    patch_counts: [200, 400]
+    # the patch count whose fit gives the samples used downstream
+    fiducial_npatch: 400
+    # beta(T_PSF) = beta0 + beta1 (T_PSF / <T_PSF> - 1)
+    beta_psf_size: false
 
   # Fiducial scale cut (arcmin) for B-mode PTE extraction. Paper II: [12, 83].
   fiducial_scale_cut: [12, 83]

pyproject.toml

@@ -81,11 +81,8 @@ dependencies = [
     "sacc>=2.4,<3",
     "scipy",
     "seaborn",
-    # Intentionally tracks the mutable `develop` branch (NOT SHA-pinned):
-    # shear_psf_leakage is actively co-developed CosmoStat code we want to stay
-    # current with, so freezing it would only create bump-churn. (Contrast the
-    # getdist feature-branch below, which is an external fork we pin for repro.)
-    "shear_psf_leakage @ git+https://github.com/CosmoStat/shear_psf_leakage.git@develop",
+    # Centred galaxy-position alpha template (shear_psf_leakage PR #59).
+    "shear_psf_leakage @ git+https://github.com/CosmoStat/shear_psf_leakage.git@d0e10d83ecb1978ed56599f0b46326e3d000680e",
     "skyproj",
     # UNIONS-WL fork of DESC Smokescreen. Tracks the fork's main; uv.lock
     # freezes the exact commit, so updates land via

src/sp_validation/correlation.py

@@ -1,11 +1,30 @@
 """Layout-independent full-sample means with patched covariance measurements."""
 
+import hashlib
 import json
 from pathlib import Path
 
+import numpy as np
 import treecorr
 
 
+def selection_digest(*masks):
+    """Digest of row selections; ``None`` stands for every row.
+
+    Recorded beside a cached measurement so a product made from different rows,
+    e.g. another star flag or tomographic bin, is never reused.
+    """
+    digest = hashlib.sha1()
+    for mask in masks:
+        if mask is None:
+            digest.update(b"all;")
+        else:
+            mask = np.asarray(mask, dtype=bool)
+            digest.update(f"{mask.size};".encode())
+            digest.update(np.packbits(mask).tobytes())
+    return digest.hexdigest()
+
+
 def _measurement_configs(config):
     """Return the explicit means settings and default-tolerance patched settings."""
     means = dict(config)
@@ -14,15 +33,19 @@ def _measurement_configs(config):
     return {"means": means, "patched": patched}
 
 
-def measurement_matches(path, config):
-    """Whether a cache records the current means and patched configurations."""
+def measurement_matches(path, config, selection=None, definition=None):
+    """Whether a cache records the current configurations and row selection."""
     metadata = Path(str(path) + ".json")
     if not Path(path).exists() or not metadata.exists():
         return False
     recorded = json.loads(metadata.read_text())
+    if recorded.get("definition") != definition:
+        return False
     recorded_configs = recorded.get("configs")
     if recorded.get("means") != "unpatched" or not isinstance(recorded_configs, dict):
         return False
+    if recorded.get("selection") != selection:
+        return False
     expected_configs = _measurement_configs(config)
     if config.get("min_top") is not None:
         # With the root depth pinned, threads affect reduction round-off,
@@ -33,12 +56,14 @@ def measurement_matches(path, config):
     return recorded_configs == expected_configs
 
 
-def write_measurement_metadata(path, config):
-    """Record both TreeCorr configurations beside a cached product."""
-    Path(str(path) + ".json").write_text(
-        json.dumps({"means": "unpatched", "configs": _measurement_configs(config)})
-        + "\n"
-    )
+def write_measurement_metadata(path, config, selection=None, definition=None):
+    """Record both TreeCorr configurations and the row selection beside a product."""
+    record = {"means": "unpatched", "configs": _measurement_configs(config)}
+    if definition is not None:
+        record["definition"] = definition
+    if selection is not None:
+        record["selection"] = selection
+    Path(str(path) + ".json").write_text(json.dumps(record) + "\n")
 
 
 def _unpatched_catalog(cat):

src/sp_validation/cosmo_val/core.py

@@ -22,7 +22,7 @@
 from .catalog_characterization import CatalogCharacterizationMixin
 from .cosebis import CosebisMixin
 from .pseudo_cl import PseudoClMixin
-from .psf_systematics import PSFSystematicsMixin
+from .psf_systematics import PSF_LEAKAGE_DEFAULTS, PSFSystematicsMixin
 from .pure_eb import PureEBMixin
 from .real_space import RealSpaceMixin
 
@@ -128,8 +128,18 @@ class CosmologyValidation(
         and analysis settings for each version.
     output_dir : str, optional
         Output directory. If None, the catalog config's paths.output.
-    rho_tau_method : {'lsq', 'mcmc'}, default 'lsq'
-        Fitting method for PSF leakage systematics parameters.
+    rho_tau_method : {'eiv', 'lsq', 'emcee', 'none'}, default 'eiv'
+        Fitting method for the PSF-leakage parameters (alpha, beta, eta).
+        'eiv' is the joint errors-in-variables fit of the rho/tau jackknife
+        measurement configured by ``psf_leakage``; 'lsq' and 'emcee' fit the
+        published tau statistics with the ``cov_estimate_method`` covariance.
+    psf_leakage : dict, optional
+        Settings of the PSF-leakage measurement and 'eiv' fit, over
+        ``PSF_LEAKAGE_DEFAULTS``: ``npatch_fine`` (k-means patches shared by
+        stars and galaxies), ``patch_counts`` (jackknife patch counts the fit
+        reports, unions of the fine patches), ``fiducial_npatch`` (the count
+        whose fit gives the samples used downstream), ``beta_psf_size`` (fit
+        beta(T_PSF) = beta0 + beta1 (T_PSF / <T_PSF> - 1)).
     cov_estimate_method : {'th', 'jk', 'sim'}, default 'th'
         Covariance estimation: 'th' for semi-analytic theory, 'jk' for jackknife,
         'sim' for a covariance measured on simulations.
@@ -284,7 +294,8 @@ def __init__(
         versions,
         catalog_config="./cat_config.yaml",
         output_dir=None,
-        rho_tau_method="lsq",
+        rho_tau_method="eiv",
+        psf_leakage=None,
         cov_estimate_method="th",
         compute_cov_rho=True,
         n_cov=100,
@@ -318,7 +329,20 @@ def __init__(
         b_target=0.01,
     ):
         self.b_target = b_target
-        self.rho_tau_method = rho_tau_method
+        self.set_psf_parameter_sampling_method(rho_tau_method)
+        unknown = set(psf_leakage or {}) - set(PSF_LEAKAGE_DEFAULTS)
+        if unknown:
+            raise ValueError(f"Unknown psf_leakage settings: {sorted(unknown)}")
+        self.psf_leakage = {**PSF_LEAKAGE_DEFAULTS, **(psf_leakage or {})}
+        npatch_fine = self.psf_leakage["npatch_fine"]
+        patch_counts = self.psf_leakage["patch_counts"]
+        if not all(1 < n < npatch_fine for n in patch_counts):
+            raise ValueError(
+                f"psf_leakage patch_counts {patch_counts} must each lie strictly "
+                f"between 1 and npatch_fine={npatch_fine}"
+            )
+        if self.psf_leakage["fiducial_npatch"] not in patch_counts:
+            raise ValueError("psf_leakage fiducial_npatch must be one of patch_counts")
         self.cov_estimate_method = cov_estimate_method
         self.compute_cov_rho = compute_cov_rho
         self.n_cov = n_cov

src/sp_validation/cosmo_val/diagnostic_plots.py

@@ -93,20 +93,39 @@ def plot_objectwise_leakage(coefficients, styles, output):
     plt.close(fig)
 
 
+def leakage_contour_axes(nparam):
+    """getdist names, labels and reference markers of the leakage parameters.
+
+    Three parameters are (alpha, beta, eta); four, with the PSF-size slope,
+    (alpha, beta0, eta, beta1).
+    """
+    if nparam == 3:
+        labels = [r"\alpha", r"\beta", r"\eta"]
+    elif nparam == 4:
+        labels = [r"\alpha", r"\beta_0", r"\eta", r"\beta_1"]
+    else:
+        raise ValueError(f"expected 3 or 4 leakage parameters, got {nparam}")
+    names = [f"x{i}" for i in range(nparam)]
+    markers = dict(zip(names, (0, 1, 1, 0)))
+    return names, labels, markers
+
+
 def plot_rho_tau_fits(products, styles, contours, samples, summary):
     """Render posterior contours and propagated PSF errors from cached arrays."""
     from shear_psf_leakage import plots as psf_plots
 
     versions = list(products)
+    sample_list = [products[v]["flat_samples"] for v in versions]
+    names, labels, markers = leakage_contour_axes(sample_list[0].shape[1])
     psf_plots.plot_contours(
-        [products[v]["flat_samples"] for v in versions],
-        names=["x0", "x1", "x2"],
-        labels=[r"\alpha", r"\beta", r"\eta"],
+        sample_list,
+        names=names,
+        labels=labels,
         savefig=str(contours),
         legend_labels=versions,
         legend_loc="upper right",
         contour_colors=[styles[v]["colour"] for v in versions],
-        markers={"x0": 0, "x1": 1, "x2": 1},
+        markers=markers,
         show=False,
         close=True,
     )