Derive reporting xi means from an exactly nested fine grid

papers/bmodes/config/config.yaml

@@ -34,9 +34,10 @@ fiducial:
     B: SP_v1.4.6.3_B_leak_corr
     C: SP_v1.4.6.3_C_leak_corr
 
-  min_sep_int: 0.5
-  max_sep_int: 300.0
-  nbins_int: 1000
+  # k=34, with outward-rounded 0.08–300 arcmin coverage.
+  min_sep_int: 0.07939053012074941
+  max_sep_int: 301.3330747561176
+  nbins_int: 1015
 
   # COSEBIS parameters
   # Compute 20 modes, evaluate PTEs for both 6 and 20

papers/cosmo_val/config/config.yaml

@@ -81,9 +81,11 @@ cosmo_val:
 
   # The fine ξ± grid both B-mode statistics (COSEBIs, pure-E/B) integrate over.
   integration:
-    min_sep: 0.08
-    max_sep: 300
-    nbins: 1000
+    # k=34; reporting edges are fine edges 312, 346, ..., 992.
+    # Coverage 0.08–300 arcmin is rounded outward to whole fine bins.
+    min_sep: 0.07939053012074941
+    max_sep: 301.3330747561176
+    nbins: 1015
 
   # COSEBIs on the integration grid: mode count and the scale cuts scanned.
   cosebis:
@@ -103,9 +105,9 @@ fiducial:
   min_sep: 1.0
   max_sep: 250.0
   nbins: 20
-  min_sep_int: 0.5
-  max_sep_int: 300.0
-  nbins_int: 1000
+  min_sep_int: 0.07939053012074941
+  max_sep_int: 301.3330747561176
+  nbins_int: 1015
   nmodes: 20
   mock_version: SP_v1.4.6.3
   fiducial_min_scale: 12

src/sp_validation/angular_binning.py

@@ -0,0 +1,45 @@
+"""Logarithmic angular-grid contracts, usable by the host workflow."""
+
+import math
+
+
+def validate_nested_grids(reporting, integration):
+    """Require reporting bins to be exact unions of uniform fine bins.
+
+    Both grids use ``min_sep``, ``max_sep`` and ``nbins`` in the same units.
+    Return the integer oversampling factor. The tolerance allows floating-point
+    round-off in logarithms, not snapping a requested edge onto a nearby one.
+    """
+
+    def grid(config):
+        lo, hi, n = (
+            float(config["min_sep"]),
+            float(config["max_sep"]),
+            int(config["nbins"]),
+        )
+        if config.get("bin_type", "Log") != "Log":
+            raise ValueError("nested ξ grids require logarithmic binning")
+        if not (0 < lo < hi and n > 0 and n == float(config["nbins"])):
+            raise ValueError(
+                "log grids require 0 < min_sep < max_sep and positive integer nbins"
+            )
+        return lo, hi, n, math.log(hi / lo) / n
+
+    rlo, _, rn, rb = grid(reporting)
+    flo, _, fn, fb = grid(integration)
+    factor = rb / fb
+    offset = math.log(rlo / flo) / fb
+    tol = 1e-10
+    prefix = "integration grid must nest into reporting grid: "
+    if factor < 1 or not math.isclose(factor, round(factor), abs_tol=tol, rel_tol=0):
+        raise ValueError(
+            prefix
+            + f"reporting/fine bin_size ratio {factor:.12g} is not a positive integer"
+        )
+    if not math.isclose(offset, round(offset), abs_tol=tol, rel_tol=0):
+        raise ValueError(
+            prefix + f"reporting min_sep lies at non-integer fine edge {offset:.12g}"
+        )
+    if offset < -tol or offset + rn * factor > fn + tol:
+        raise ValueError(prefix + "the fine grid does not cover every reporting edge")
+    return int(round(factor))

src/sp_validation/correlation.py

@@ -3,8 +3,12 @@
 import json
 from pathlib import Path
 
+import numpy as np
 import treecorr
 
+from sp_validation.angular_binning import validate_nested_grids
+from sp_validation.b_modes import _reporting_binning
+
 
 def _measurement_configs(config):
     """Return the explicit means settings and default-tolerance patched settings."""
@@ -64,7 +68,7 @@ def measure_with_patches(measure, catalogs, config, means=None):
     its explicit ``bin_slop``. The patched pass omits ``bin_slop`` so TreeCorr
     uses its standard default. A cached ``means`` can be reused for repeated
     covariance draws of the same catalogue. With no patches, one means-config
-    measurement supplies both products.
+    measurement supplies both products unless supplied means are reused.
     """
     configs = _measurement_configs(config)
     if all(cat.npatch == 1 for cat in catalogs.values()):
@@ -84,13 +88,39 @@ def measure_with_patches(measure, catalogs, config, means=None):
     return means, patched
 
 
-def process_gg(config, cat1, cat2=None):
-    """Measure GG means with ``config`` and patched resampling at TreeCorr's default.
+def rebin_gg_means(fine, config):
+    """Pair-weighted fine ξ± means on exactly nested reporting edges.
 
-    The shared two-pass mechanism omits ``bin_slop`` only from the patched
-    configuration. Covariance estimates, including joint and derived-statistic
-    jackknives, retain TreeCorr's per-patch results. Published pair counts,
-    weights, separations and complex correlations come from the means pass.
+    The same operator used by pure-E/B averages the complex correlations and
+    separations. Pair counts and weights are sums over the corresponding fine
+    bins. This object supplies means only, not a covariance or patch results.
+    """
+    reporting = treecorr.GGCorrelation({**config, "var_method": "shot"})
+    if fine.sep_units != reporting.sep_units:
+        raise ValueError("fine and reporting ξ± must use the same separation units")
+    validate_nested_grids(
+        config, dict(min_sep=fine.min_sep, max_sep=fine.max_sep, nbins=fine.nbins)
+    )
+    fine_edges = np.append(fine.left_edges, fine.right_edges[-1])
+    report_edges = np.append(reporting.left_edges, reporting.right_edges[-1])
+    operator, edges = _reporting_binning(fine.weight, fine_edges, report_edges)
+    for name in ("xip", "xim", "xip_im", "xim_im", "meanr", "meanlogr"):
+        getattr(reporting, name)[:] = operator @ getattr(fine, name)
+    indices = np.searchsorted(fine_edges, edges)
+    for i, (lo, hi) in enumerate(zip(indices[:-1], indices[1:])):
+        reporting.weight[i] = np.sum(fine.weight[lo:hi])
+        reporting.npairs[i] = np.sum(fine.npairs[lo:hi])
+    return reporting
+
+
+def process_gg(config, cat1, cat2=None, means=None):
+    """Measure GG means and patched resampling products at distinct tolerances.
+
+    The configured ``bin_slop`` applies to unpatched means; the shared
+    two-pass mechanism omits it from the patched configuration so TreeCorr's
+    default sets the covariance and resampling tolerance. Supplied ``means``
+    replace the unpatched reporting pass, for example with pair-weighted fine-grid
+    means. The patched covariance is materialized before the means fields change.
     """
     catalogs = {"cat1": cat1}
     if cat2 is not None:
@@ -104,10 +134,9 @@ def measure(cats, measurement_config):
         gg.process(cats["cat1"], cat2=cats.get("cat2"))
         return gg
 
-    means, gg = measure_with_patches(measure, catalogs, config)
+    means, gg = measure_with_patches(measure, catalogs, config, means=means)
     if means is not gg:
-        # TreeCorr estimates covariance lazily from fields that the means pass
-        # replaces below, so freeze the patched-pass estimates first.
+        # Freeze the patched covariance and variances before replacing means.
         _ = gg.varxip, gg.varxim, gg.cov
         for name in (
             "xip",

src/sp_validation/cosmo_val/core.py

@@ -11,6 +11,7 @@
 from shear_psf_leakage import leakage, run_object, run_scale
 
 from .. import io
+from ..angular_binning import validate_nested_grids
 from ..b_modes import (
     _get_pte_from_scale_cut,
     covariance_label,
@@ -146,6 +147,10 @@ class CosmologyValidation(
         Maximum angular separation in arcminutes for correlation function binning.
     nbins : int, default 20
         Number of angular bins for TreeCorr real-space correlation functions.
+    integration : dict, optional
+        Unpatched fine ξ± grid (min_sep, max_sep, nbins) supplying the published
+        reporting means. Its bins must nest into the constructor's reporting
+        grid. Without it, reporting means use an unpatched reporting pass.
     b_target : float, default 0.01
         Maximum logarithmic-bin tolerance for full-sample means, with bin_slop
         capped at 1 per grid. Patched covariance passes use TreeCorr's default.
@@ -315,6 +320,7 @@ def __init__(
         compute_tomography=False,
         force_run=False,
         b_target=0.01,
+        integration=None,
     ):
         self.b_target = b_target
         self.rho_tau_method = rho_tau_method
@@ -388,6 +394,18 @@ def __init__(
         }
 
         self.treecorr_config = self._binning()
+        self.integration = None
+        if integration is not None:
+            if integration.get("npatch", 1) != 1:
+                raise ValueError(
+                    "the integration grid supplying ξ± means must use npatch=1"
+                )
+            self.integration = {
+                "min_sep": float(integration["min_sep"]),
+                "max_sep": float(integration["max_sep"]),
+                "nbins": int(integration["nbins"]),
+            }
+            validate_nested_grids(self.treecorr_config, self.integration)
 
         self.catalog_config_path = Path(catalog_config)
         with self.catalog_config_path.open("r") as file:
@@ -624,6 +642,15 @@ def _binning(self, min_sep=None, max_sep=None, nbins=None, **extra):
         config["bin_slop"] = extra.get("bin_slop", min(1, self.b_target / bin_size))
         return config
 
+    def _integration_binning(self, min_sep, max_sep, nbins, *, default):
+        """Configured fine grid, with a statistic's standalone fallback bounds."""
+        grid = self.integration or dict(zip(("min_sep", "max_sep", "nbins"), default))
+        return self._binning(
+            grid["min_sep"] if min_sep is None else min_sep,
+            grid["max_sep"] if max_sep is None else max_sep,
+            grid["nbins"] if nbins is None else nbins,
+        )
+
     def _read_shear_cols(self, ver, *keys):
         """Read shear-catalog columns by their config-key names.
 

src/sp_validation/cosmo_val/cosebis.py

@@ -21,9 +21,9 @@ class CosebisMixin:
     def calculate_cosebis(
         self,
         version,
-        min_sep_int=0.5,
-        max_sep_int=500,
-        nbins_int=1000,
+        min_sep_int=None,
+        max_sep_int=None,
+        nbins_int=None,
         npatch=None,
         nmodes=10,
         cov_path=None,
@@ -48,13 +48,13 @@ def calculate_cosebis(
             The catalog version to compute the COSEBIs for.
         min_sep_int : float, optional
             Minimum separation for integration binning (fine binning for COSEBIs).
-            Defaults to 0.5 arcmin.
+            Defaults to the configured integration grid, or 0.5 arcmin.
         max_sep_int : float, optional
             Maximum separation for integration binning (fine binning for COSEBIs).
-            Defaults to 500 arcmin.
+            Defaults to the configured integration grid, or 500 arcmin.
         nbins_int : int, optional
             Number of bins for integration binning (fine binning for COSEBIs).
-            Defaults to 1000.
+            Defaults to the configured integration grid, or 1000.
         npatch : int, optional
             Number of patches for the jackknife resampling. Defaults to self.npatch.
         nmodes : int, optional
@@ -102,7 +102,12 @@ def calculate_cosebis(
         npatch = npatch or self.npatch
 
         # Always use integration binning for COSEBIs calculation (fine binning)
-        treecorr_config = self._binning(min_sep_int, max_sep_int, nbins_int)
+        treecorr_config = self._integration_binning(
+            min_sep_int, max_sep_int, nbins_int, default=(0.5, 500, 1000)
+        )
+        min_sep_int, max_sep_int, nbins_int = (
+            treecorr_config[k] for k in ("min_sep", "max_sep", "nbins")
+        )
 
         # Calculate single fine-binned correlation function for COSEBIs
         print(
@@ -171,9 +176,9 @@ def plot_cosebis(
         self,
         version=None,
         output_dir=None,
-        min_sep_int=0.5,
-        max_sep_int=500,
-        nbins_int=1000,  # Integration binning
+        min_sep_int=None,
+        max_sep_int=None,
+        nbins_int=None,  # Integration binning
         npatch=None,
         nmodes=10,
         cov_path=None,
@@ -201,7 +206,7 @@ def plot_cosebis(
             Output directory for plots. Defaults to self.cc['paths']['output'].
         min_sep_int, max_sep_int, nbins_int : float, float, int
             Integration binning parameters for correlation function
-            (default: 0.5, 500, 1000)
+            (the configured grid, otherwise 0.5, 500, 1000)
         npatch : int, optional
             Number of patches for jackknife covariance. Defaults to instance value.
         nmodes : int
@@ -234,6 +239,13 @@ def plot_cosebis(
         # Determine variance method based on whether theoretical covariance is used
         var_method = "analytic" if cov_path is not None else "jackknife"
 
+        integration = self._integration_binning(
+            min_sep_int, max_sep_int, nbins_int, default=(0.5, 500, 1000)
+        )
+        min_sep_int, max_sep_int, nbins_int = (
+            integration[k] for k in ("min_sep", "max_sep", "nbins")
+        )
+
         # Create output filename with integration parameters to match Snakemake
         out_stub = (
             f"{output_dir}/{version}_cosebis_minsep={min_sep_int}_"

src/sp_validation/cosmo_val/pure_eb.py

@@ -7,6 +7,7 @@
 
 import numpy as np
 
+from ..angular_binning import validate_nested_grids
 from ..b_modes import (
     calculate_eb_statistics,
     calculate_pure_eb_correlation,
@@ -26,9 +27,9 @@ def calculate_pure_eb(
         min_sep=None,
         max_sep=None,
         nbins=None,
-        min_sep_int=0.08,
-        max_sep_int=300,
-        nbins_int=1000,
+        min_sep_int=None,
+        max_sep_int=None,
+        nbins_int=None,
         npatch=None,
         compute_tomography=False,
         cov_path_int=None,
@@ -48,7 +49,8 @@ def calculate_pure_eb(
         min_sep, max_sep, nbins : float, float, int, optional
             Reporting binning. Default to the values in self.treecorr_config.
         min_sep_int, max_sep_int, nbins_int : float, float, int, optional
-            Integration binning (default: 0.08-300 arcmin, 1000 bins). It must
+            Integration binning (the configured grid, otherwise 0.08–300 arcmin,
+            1000 bins). It must
             extend beyond the reporting range on both sides.
         npatch : int, optional
             Jackknife patch count. Defaults to self.npatch.
@@ -83,11 +85,16 @@ def calculate_pure_eb(
         reporting_edges = np.geomspace(
             reporting["min_sep"], reporting["max_sep"], reporting["nbins"] + 1
         )
+        integration = self._integration_binning(
+            min_sep_int, max_sep_int, nbins_int, default=(0.08, 300, 1000)
+        )
+        if self.integration is not None:
+            validate_nested_grids(reporting, integration)
         ggs_int = self.calculate_2pcf_version(
             version,
             npatch=npatch,
             compute_tomography=compute_tomography,
-            **self._binning(min_sep_int, max_sep_int, nbins_int),
+            **integration,
         )
 
         results = {}
@@ -119,9 +126,9 @@ def plot_pure_eb(
         min_sep=None,
         max_sep=None,
         nbins=None,
-        min_sep_int=0.08,
-        max_sep_int=300,
-        nbins_int=1000,
+        min_sep_int=None,
+        max_sep_int=None,
+        nbins_int=None,
         npatch=None,
         compute_tomography=False,
         cov_path_int=None,
@@ -150,7 +157,7 @@ def plot_pure_eb(
             Binning parameters for reporting scale. Uses treecorr_config if None.
         min_sep_int, max_sep_int, nbins_int : float, float, int
             Binning parameters for integration scale
-            (default: 0.08-300 arcmin, 1000 bins)
+            (the configured grid, otherwise 0.08–300 arcmin, 1000 bins)
         npatch : int, optional
             Number of patches for jackknife covariance. Uses self.npatch if None.
         compute_tomography : bool, optional
@@ -182,6 +189,12 @@ def plot_pure_eb(
         output_dir = output_dir or self.cc["paths"]["output"]
         npatch = npatch or self.npatch
 
+        integration = self._integration_binning(
+            min_sep_int, max_sep_int, nbins_int, default=(0.08, 300, 1000)
+        )
+        min_sep_int, max_sep_int, nbins_int = (
+            integration[k] for k in ("min_sep", "max_sep", "nbins")
+        )
         var_method = "jackknife" if cov_path_int is None else "analytic"
 
         # Use treecorr_config defaults for reporting scale binning