Feature/305 adapt tomographic 2pcf script to sacc io format

src/sp_validation/cosmo_val/pure_eb.py

@@ -82,7 +82,12 @@ def calculate_pure_eb(
         Returns
         -------
         dict
-            A dictionary containing the following keys:
+            Mapping of ``"tomo_bin_{b1}_tomo_bin_{b2}"`` to that bin pair's pure
+            E/B results, following the keys of :meth:`calculate_2pcf_version`.
+            With ``compute_tomography=False`` the single key is
+            ``"tomo_bin_all_tomo_bin_all"``.
+
+            Each value is a dictionary containing the following keys:
 
             - "xip_E": Pure E-mode correlation function for xi+.
             - "xim_E": Pure E-mode correlation function for xi-.
@@ -164,6 +169,7 @@ def plot_pure_eb(
         max_sep_int=300,
         nbins_int=1000,
         npatch=None,
+        compute_tomography=False,  # LG: Hardcoded to False for now
         var_method="jackknife",
         cov_path_int=None,
         cosmo_cov=None,
@@ -272,20 +278,31 @@ def plot_pure_eb(
             )
 
             # Get or calculate results for this version
-            version_results = results_list[idx] or self.calculate_pure_eb(
-                version,
-                min_sep=min_sep,
-                max_sep=max_sep,
-                nbins=nbins,
-                min_sep_int=min_sep_int,
-                max_sep_int=max_sep_int,
-                nbins_int=nbins_int,
-                npatch=npatch,
-                var_method=var_method,
-                cov_path_int=cov_path_int,
-                cosmo_cov=cosmo_cov,
-                n_samples=n_samples,
-            )
+            version_results = results_list[idx]
+            if version_results is None:
+                version_results_tomo = self.calculate_pure_eb(
+                    version,
+                    min_sep=min_sep,
+                    max_sep=max_sep,
+                    nbins=nbins,
+                    min_sep_int=min_sep_int,
+                    max_sep_int=max_sep_int,
+                    nbins_int=nbins_int,
+                    compute_tomography=False,  # LG: Hardcoded to False for now
+                    npatch=npatch,
+                    var_method=var_method,
+                    cov_path_int=cov_path_int,
+                    cosmo_cov=cosmo_cov,
+                    n_samples=n_samples,
+                )
+                version_results = version_results_tomo[
+                    "tomo_bin_all_tomo_bin_all"
+                ]  # LG: Extract non-tomographic results
+            elif (
+                isinstance(version_results, dict)
+                and "tomo_bin_all_tomo_bin_all" in version_results
+            ):
+                version_results = version_results["tomo_bin_all_tomo_bin_all"]
 
             # Calculate E/B statistics for all bin combinations
             version_results = calculate_eb_statistics(

src/sp_validation/cosmo_val/real_space.py

@@ -7,12 +7,33 @@
 """
 
 import os
+import re
 
 import matplotlib.pyplot as plt
 import numpy as np
 import treecorr
 from astropy.io import fits
 
+from .. import sacc_io
+
+# calculate_2pcf_version keys its results "tomo_bin_{b1}_tomo_bin_{b2}"; this
+# reads the bin pair back out ("all" for the non-tomographic single pair).
+_PAIR_KEY = re.compile(r"tomo_bin_(.+?)_tomo_bin_(.+)")
+
+
+def _sacc_bin_index(ggs):
+    """Map each catalog tomographic bin id to its 0-based SACC bin index.
+
+    The catalog numbers its bins however the ``tomo_bin_col`` column does
+    (1-based in practice, and the single id ``"all"`` for a non-tomographic
+    run), while SACC tracers are ``source_{i}`` counted from zero. The map is
+    built from the ids actually present, in ascending order, so it never
+    assumes the catalog starts at 1 or numbers its bins without gaps.
+    """
+    ids = {b for key in ggs for b in _PAIR_KEY.fullmatch(key).groups()}
+    order = sorted(ids, key=lambda b: 0 if b == "all" else int(b))
+    return {b: i for i, b in enumerate(order)}
+
 
 class RealSpaceMixin:
     def calculate_2pcf_version(
@@ -71,8 +92,6 @@ def calculate_2pcf_version(
 
         ggs = {f"tomo_bin_{b1}_tomo_bin_{b2}": None for b1, b2 in tomo_bin_pairs}
 
-        # LG TO-DO: Change to sacc_io method
-
         patch_file = self._output_path(f"{ver}_patches_npatch={npatch}.dat")
 
         cat_gal = fits.getdata(self.cc[ver]["shear"]["path"])
@@ -130,8 +149,8 @@ def calculate_2pcf_version(
 
     def calculate_2pcf(
         self,
-        npatch=None,
         compute_tomography=False,
+        npatch=None,
         **treecorr_config,
     ):
         """
@@ -159,15 +178,123 @@ def calculate_2pcf(
         for ver in self.versions:
             self.cat_ggs[ver] = self.calculate_2pcf_version(
                 ver,
-                npatch=npatch,
                 compute_tomography=compute_tomography,
+                npatch=npatch,
                 **treecorr_config,
             )
 
-        # LG TO-DO: No longer writing out text file, change to sacc_io method
-
         return self.cat_ggs
 
+    def save_2pcf_sacc(
+        self,
+        ver,
+        sacc_path,
+        ggs=None,
+        *,
+        type,
+        grid="reporting",
+        metadata=None,
+    ):
+        """Write the measured ξ± for ``ver`` to ``sacc_path`` as a SACC file.
+
+        Serialises the TreeCorr output of :meth:`calculate_2pcf_version` (or of
+        :meth:`calculate_2pcf`, via ``self.cat_ggs``) into the standard layout
+        of :mod:`sp_validation.sacc_io`: the ``source_{i}`` n(z) tracers, one
+        ξ+/ξ− block per tomographic bin pair, and the covariance.
+
+        Insertion order is load-bearing. ``sacc_io.add_xi`` writes one bin pair
+        as ``[ξ+; ξ−]``, so the data vector is *pair-major*, and the covariance
+        is tied to it by position alone. Both the ξ insertion and the
+        covariance therefore iterate the same ``pairs`` list, sorted by SACC
+        bin index — never recomputed independently.
+
+        The covariance follows what the measurement can support, read off the
+        correlations themselves (``var_method``) rather than off ``npatch``,
+        which the caller may have overridden per call:
+
+        - jackknife (npatch > 1): ``treecorr.estimate_multi_cov`` over the
+          correlations in ``pairs`` order. TreeCorr concatenates each one as
+          ``[ξ+; ξ−]``, which is exactly the pair-major insertion order, so
+          the result is one contiguous block spanning every ξ point —
+          including the cross-pair covariance, which a per-pair ``gg.cov``
+          would drop.
+        - shot noise (npatch == 1): the ``varxip``/``varxim`` diagonal, stored
+          as a SACC ``DiagonalCovariance``.
+
+        Parameters
+        ----------
+        ver : str
+            Catalog version, used for the n(z) lookup and stamped as metadata.
+        sacc_path : str
+            Path to the SACC product.
+        ggs : dict, optional
+            ``{"tomo_bin_{b1}_tomo_bin_{b2}": treecorr.GGCorrelation}`` as
+            returned by :meth:`calculate_2pcf_version`. Defaults to
+            ``self.cat_ggs[ver]``, i.e. the last :meth:`calculate_2pcf` run.
+        type : {'data', 'mock'}
+            Provenance of the input catalog, required by ``sacc_io.save``.
+            No default: only the caller knows whether it ran on a mock, and
+            ``sacc_io.load`` refuses unblinded ``type='data'`` files.
+        grid : str, optional
+            The ``grid`` tag on every point (default ``'reporting'``). Use
+            ``'integration'`` for the fine grid COSEBIs / pure-EB integrate
+            over, which shares this data type and tracer pair and is told
+            apart by nothing else.
+        metadata : dict, optional
+            Extra key/value pairs stored on the file, merged over the
+            version/npatch/blind stamped here.
+
+        Returns
+        -------
+        sacc.Sacc
+            The written data set.
+        """
+        ggs = self.cat_ggs[ver] if ggs is None else ggs
+        if not ggs:
+            raise ValueError(f"{ver}: no ξ± measurements to write")
+
+        index = _sacc_bin_index(ggs)
+        pairs = sorted(
+            ggs, key=lambda key: [index[b] for b in _PAIR_KEY.fullmatch(key).groups()]
+        )
+
+        z, *nz_cols = self.get_redshift(ver)
+        if len(nz_cols) != len(index):
+            raise ValueError(
+                f"{ver}: the n(z) file has {len(nz_cols)} distribution "
+                f"column(s) but the measurement covers {len(index)} "
+                f"tomographic bin(s) — every source bin needs its own n(z)"
+            )
+
+        s = sacc_io.new_sacc(
+            [(z, nz) for nz in nz_cols],
+            metadata={
+                "version": ver,
+                "npatch": int(ggs[pairs[0]].npatch1),
+                **({"blind": self.blind} if self.blind is not None else {}),
+                **(metadata or {}),
+            },
+        )
+
+        for key in pairs:
+            gg = ggs[key]
+            bins = tuple(index[b] for b in _PAIR_KEY.fullmatch(key).groups())
+            sacc_io.add_xi(
+                s,
+                bins,
+                gg.meanr,
+                gg.xip,
+                gg.xim,
+                grid=grid,
+                theta_nom=gg.rnom,
+                npairs=gg.npairs,
+                weight=gg.weight,
+            )
+
+        sacc_io.save(s, sacc_path, type=type)
+        self.print_done(f"Wrote ξ± SACC for {ver} to {sacc_path}.")
+        return s
+
     def calculate_aperture_mass_dispersion(
         self,
         theta_min=0.3,

src/sp_validation/tests/test_cosmo_val.py

@@ -631,6 +631,13 @@ def test_calculate_pure_eb_runs_on_synthetic_catalog(self, tmp_path):
             nbins_int=600,
         )
 
+        # calculate_pure_eb keys its results by tomographic bin pair, mirroring
+        # calculate_2pcf_version. Non-tomographic (compute_tomography=False, the
+        # default) means exactly one pair -- the all-galaxy auto-correlation --
+        # so assert that and unwrap it to the per-bin mode dict.
+        assert list(results) == ["tomo_bin_all_tomo_bin_all"]
+        results = results["tomo_bin_all_tomo_bin_all"]
+
         # Reference mode vectors from the seeded synthetic catalog + Schneider
         # transform. Deterministic (full-sample treecorr, no RNG); regenerate by
         # running calculate_pure_eb with the setup above and printing repr() of

src/sp_validation/tests/test_glass_mock.py

@@ -40,7 +40,23 @@
 
 camb = pytest.importorskip("camb")
 
-HAVE_GLASS = importlib.util.find_spec("glass") is not None
+def _have(module):
+    """True if ``module`` is importable. False if it or a parent is missing."""
+    try:
+        return importlib.util.find_spec(module) is not None
+    except ModuleNotFoundError:
+        # find_spec imports the parent package, which raises if it is absent.
+        return False
+
+
+HAVE_GLASS = _have("glass")
+
+# The map path needs more than GLASS itself: build_shells and
+# Cosmology_from_camb import ``cosmology.compat.camb``, which ships with the
+# ``cosmology`` API package rather than with GLASS. Gate on it separately so a
+# missing dependency reports as a skip instead of masquerading as the (real,
+# but different) glass/cosmology API incompatibility tracked in the xfail below.
+HAVE_COSMOLOGY_CAMB = _have("cosmology.compat.camb")
 
 REFERENCE = Path(__file__).parent / "data" / "glass_mock_camb_reference.npz"
 
@@ -131,6 +147,10 @@ def test_config_change_breaks_reference():
 
 
 @pytest.mark.skipif(not HAVE_GLASS, reason="GLASS not installed in this image")
+@pytest.mark.skipif(
+    not HAVE_COSMOLOGY_CAMB,
+    reason="cosmology.compat.camb not installed in this image",
+)
 @pytest.mark.xfail(
     reason=(
         "glass_mock map path is incompatible with the installed glass/cosmology "

workflow/common.py

@@ -5,19 +5,13 @@
 import re
 from pathlib import Path
 
-# Output roots are env-overridable so a reproduction run can write into a
-# fresh tree without clobbering (or silently reusing) prior products.
-COSMO_VAL = Path(
-    os.environ.get(
-        "COSMO_VAL", "/n17data/cdaley/unions/code/sp_validation/cosmo_val/output"
-    )
-)
+SP_VALIDATION = Path(__file__).resolve().parents[1]
+COSMO_VAL = Path(os.environ.get("COSMO_VAL", SP_VALIDATION / "results/cosmo_val"))
 COSMO_INFERENCE = Path(
-    os.environ.get(
-        "COSMO_INFERENCE", "/n17data/cdaley/unions/code/sp_validation/cosmo_inference"
-    )
+    os.environ.get("COSMO_INFERENCE", SP_VALIDATION / "cosmo_inference")
 )
-CAT_CONFIG = "/n17data/cdaley/unions/code/sp_validation/cosmo_val/cat_config.yaml"
+CAT_CONFIG = SP_VALIDATION / "cosmo_val/cat_config.yaml"
+
 BLINDS = ["A", "B", "C"]
 BLOCK_PAIRS = [("++", "1"), ("--", "2"), ("+-", "3")]
 
@@ -34,12 +28,14 @@
     "version": r"SP_v[\d.]+(_w_iv)?(_ecut\d+)?(_leak_corr)?",
     "blind": r"[ABC]",
     "nbins": r"\d+",
+    "npatch": r"\d+",
     "min_sep": r"[0-9.]+",
     "max_sep": r"[0-9.]+",
     "gaussian": r"(g|ng)",
     "block_pm": r"(\+\+|--|\+-)",
     "block_i": r"[123]",
     "mask_suffix": r"(_masked)?",
+    "tomo_suffix": r"(_tomo)?",
     "mock_id": r"\d{5}",
     "nside": r"\d+",
 }

workflow/rules/cosmo_val.smk

@@ -45,7 +45,7 @@ CV_BINNING = (
     f"_nbins={CV['nbins']}_npatch={CV['npatch']}"
 )
 
-
+# LG: DEPRACATED, xi_pm vectors no longer written to txt files.
 def cv_xi_txt(version):
     """Path to the 2pcf data vector calculate_2pcf writes for a version.
 
@@ -238,8 +238,9 @@ rule cv_additive_bias:
         "../scripts/cv_additive_bias.py"
 
 
+# LG: DEPRACATED since we do not plot unblinded data vectors.
 rule cv_plot_2pcf:
-    """n_pairs / xi± overlay across versions."""
+    """Tomographic xi± across versions."""
     input:
         xi=[cv_xi_txt(v) for v in CV_VERSIONS],
     output:
@@ -251,11 +252,15 @@ rule cv_plot_2pcf:
     script:
         "../scripts/cv_plot_2pcf.py"
 
-
+# LG: DEPRACATED since we do not plot unblinded data vectors.
 rule cv_ratio_xi_sys_xi:
     """Ratio of PSF systematics (xi_psf_sys) to the cosmic-shear signal (xi+)."""
     input:
-        xi=[cv_xi_txt(v) for v in CV_VERSIONS],
+        # The xi txt this used to declare is written by nothing since
+        # calculate_2pcf stopped emitting it; cv_xi_txt is gone with it. The
+        # rule recomputes xi_psf_sys itself, so the entry is simply dropped —
+        # left in place it is a parse-time NameError that takes the whole
+        # workflow down, deprecated or not.
         rho=[cv_rho_stats(v) for v in CV_VERSIONS],
         tau=[cv_tau_stats(v) for v in CV_VERSIONS],
     output:
@@ -294,8 +299,6 @@ rule cv_pseudo_cl:
 
 rule cv_pure_eb:
     """Pure E/B-mode decomposition for one version (config-space)."""
-    input:
-        xi=lambda w: cv_xi_txt(w.version),
     output:
         npz=cv_pure_eb_npz("{version}"),
     params:
@@ -316,8 +319,9 @@ rule cv_pure_eb:
 
 rule cv_cosebis:
     """COSEBIs E/B decomposition for one version (config-space, fine binning)."""
-    input:
-        xi=lambda w: cv_xi_txt(w.version),
+    # No xi input, for the same reason as cv_pure_eb: calculate_cosebis runs
+    # its own TreeCorr over the fine integration binning, and the xi txt it
+    # used to declare is no longer written by anything.
     output:
         npz=cv_cosebis_npz("{version}"),
     params: