Pure E/B: fixed-quadrature operator with exact covariance

papers/bmodes/config/config.yaml

@@ -58,14 +58,9 @@ fiducial:
 
 var_method: jackknife
 
-# covariance settings for semi-analytical E/B mode covariance
 covariance:
   # Use masked covariance by default when true
   default_masked: true
-  # Enable semi-analytical covariance propagation
-  use_semianalytic: true
-  # Number of Monte Carlo samples for covariance propagation (used by pure E/B mode analysis)
-  n_samples: 2000
   # Cosmology: imported from cs_util.cosmo.PLANCK18 (astropy Planck18)
   # Mask Cl paths: defined in covariance.smk (MASK_CLS_FILES)
 
@@ -135,11 +130,6 @@ plotting:
     fiducial_line_color: "black"
     fiducial_line_width: 1.0
 
-# Pure E/B mode analysis
-pure_eb:
-  # Number of parallel chunks for MC covariance estimation
-  n_chunks: 20
-
 # pixel mask processing for survey geometry and CosmoCov integration
 pixel_mask:
   # Target nside values for downgrading (balance accuracy vs computational cost)

papers/bmodes/rules/figures.smk

@@ -81,14 +81,6 @@ def _reporting_cov_path(version):
     return covariance_path(version, gaussian="ng")
 
 
-def _xi_reporting_path(version):
-    """Path to reporting-scale 2PCF file."""
-    return (
-        f"{COSMO_VAL_OUTPUT}/{version}_xi_minsep={FIDUCIAL['min_sep']}"
-        f"_maxsep={FIDUCIAL['max_sep']}_nbins={FIDUCIAL['nbins']}_npatch={FIDUCIAL['npatch']}.txt"
-    )
-
-
 def _xi_integration_path(version):
     """Path to fine-binned 2PCF integration file. Unpatched: values only, no covariance."""
     return (
@@ -183,50 +175,20 @@ rule cosebis_data_vector:
 # Pure E/B
 # ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
 
-# Number of parallel chunks for MC covariance estimation
-N_PURE_EB_CHUNKS = config["pure_eb"]["n_chunks"]
-
-
-rule precompute_pure_eb_chunk:
-    """Compute a chunk of MC samples for pure E/B covariance (scatter)."""
-    input:
-        cov_integration=lambda w: _cov_integration_path(w.version),
-        xi_reporting=lambda w: _xi_reporting_path(w.version),
-        xi_integration=lambda w: _xi_integration_path(w.version),
-    output:
-        "results/paper_plots/intermediate/chunks/{version}_pure_eb_chunk_{chunk_id}.npz",
-    params:
-        version="{version}",
-        chunk_id="{chunk_id}",
-        n_chunks=N_PURE_EB_CHUNKS,
-        n_samples=config["covariance"]["n_samples"],
-        cosmo_params=PLANCK18,
-        **FIDUCIAL_BINNING,
-    resources:
-        mem_mb=8000,
-    script:
-        "../scripts/precompute_pure_eb_chunk.py"
-
-
-rule precompute_pure_eb:
-    """Gather MC sample chunks and compute final pure E/B covariance."""
+rule pure_eb_modes:
+    """Pure E/B modes and their exact covariance K C_ξ Kᵀ from the fine ξ±."""
     input:
-        chunks=expand(
-            "results/paper_plots/intermediate/chunks/{{version}}_pure_eb_chunk_{chunk_id}.npz",
-            chunk_id=range(N_PURE_EB_CHUNKS),
-        ),
-        xi_reporting=lambda w: _xi_reporting_path(w.version),
         xi_integration=lambda w: _xi_integration_path(w.version),
+        cov_integration=lambda w: _cov_integration_path(w.version),
     output:
-        "results/paper_plots/intermediate/{version}_pure_eb_semianalytic.npz",
+        "results/paper_plots/intermediate/{version}_pure_eb.npz",
     params:
-        version="{version}",
         **FIDUCIAL_BINNING,
     resources:
         mem_mb=8000,
-        runtime=5,
+        runtime=20,
     script:
-        "../scripts/gather_pure_eb_chunks.py"
+        "../scripts/pure_eb_modes.py"
 
 
 rule pure_eb_data_vector:
@@ -239,7 +201,7 @@ rule pure_eb_data_vector:
     """
     input:
         # Per-version inputs: pure_eb_{version} and cov_{version} for all versions
-        **{f"pure_eb_{ver}": f"results/paper_plots/intermediate/{ver}_pure_eb_semianalytic.npz"
+        **{f"pure_eb_{ver}": f"results/paper_plots/intermediate/{ver}_pure_eb.npz"
            for ver in VERSIONS_ALL_FOR_PLOTS},
         **{f"cov_{ver}": _reporting_cov_path(ver) for ver in VERSIONS_ALL_FOR_PLOTS},
     output:
@@ -259,7 +221,7 @@ rule pure_eb_version_comparison:
     input:
         # Pure E/B only for leak-corrected versions
         pure_eb_data=[
-            f"results/paper_plots/intermediate/{ver}_pure_eb_semianalytic.npz"
+            f"results/paper_plots/intermediate/{ver}_pure_eb.npz"
             for ver in VERSIONS_LEAK_CORR
         ],
     params:
@@ -282,7 +244,7 @@ rule pure_eb_covariance:
     - Correlation structure across 6 blocks (E+/E-/B+/B-/amb+/amb-)
     """
     input:
-        pure_eb_data=f"results/paper_plots/intermediate/{FIDUCIAL_VERSION}_pure_eb_semianalytic.npz",
+        pure_eb_data=f"results/paper_plots/intermediate/{FIDUCIAL_VERSION}_pure_eb.npz",
     output:
         evidence=f"{TAPESTRY_DIR}/pure_eb_covariance/evidence.json",
         figure=f"{TAPESTRY_DIR}/pure_eb_covariance/figure.png",
@@ -292,17 +254,13 @@ rule pure_eb_covariance:
 
 
 rule calculate_pure_eb_ptes:
-    """PTE matrices for pure E/B-mode scale-cut robustness.
-
-    The PTEs are Hartlap-debiased by the MC draw count.
-    """
+    """PTE matrices for pure E/B-mode scale-cut robustness."""
     input:
-        pure_eb_data="results/paper_plots/intermediate/{version}_pure_eb_semianalytic.npz",
+        pure_eb_data="results/paper_plots/intermediate/{version}_pure_eb.npz",
     output:
         "results/paper_plots/intermediate/{version}_pure_eb_ptes.npz",
     params:
         version="{version}",
-        n_samples=config["covariance"]["n_samples"],
     resources:
         mem_mb=16000,
         runtime=30,
@@ -455,7 +413,7 @@ rule bb_covariance_nz_independence:
     """
     input:
         # Per-realisation MC-propagated pure E/B covariances
-        **{f"pure_eb_{label}": f"results/paper_plots/intermediate/{ver}_pure_eb_semianalytic.npz"
+        **{f"pure_eb_{label}": f"results/paper_plots/intermediate/{ver}_pure_eb.npz"
            for label, ver in NZ_REALISATIONS.items()},
         # COSEBIS: xi integration file (shared) + per-realisation config-space covariances
         xi_integration=_xi_integration_path(MOCK_VERSION),

papers/bmodes/scripts/bb_covariance_nz_independence.py

@@ -136,7 +136,6 @@ def make_figure(
     cosebis_results,
     reference,
     output_path,
-    n_samples=2000,
 ):
     """Four-panel figure comparing BB vs EE stability across n(z) realisations.
 
@@ -152,40 +151,26 @@ def make_figure(
     color_E = "#E69F00"  # orange for E
     color_B = "#0072B2"  # blue for B
 
-    # Expected 1σ error on ratio of two MC-estimated quantities
-    # σ(ratio) ≈ √(2/N) for ratio ≈ 1
-    ratio_err = np.sqrt(2.0 / n_samples)
-
-    def setup_ratio_panel(ax, xlabel, title, show_mc_band=False):
+    def setup_ratio_panel(ax, xlabel, title):
         ax.axhline(1.0, color="gray", ls="-", lw=0.8, zorder=0)
-        if show_mc_band:
-            ax.axhspan(
-                1 - ratio_err,
-                1 + ratio_err,
-                color="gray",
-                alpha=0.25,
-                label=rf"$\pm\sqrt{{2/N}}$ ($N={n_samples}$)",
-            )
         ax.set_xscale("log")
         ax.set_xlabel(xlabel)
         ax.set_title(title)
 
-    def plot_ratios(ax, x, results, b_key, e_key, b_name, e_name, shift, b_err):
-        """B (with optional MC error bar) and E ratios for every realisation."""
+    def plot_ratios(ax, x, results, b_key, e_key, b_name, e_name, shift):
+        """B and E ratios for every realisation."""
         for i, (label, res) in enumerate(results.items()):
             xi = shift(x, i)
             marker = MARKERS[i % len(MARKERS)]
             pair = f"{label}/{reference}"
-            ax.errorbar(
+            ax.plot(
                 xi,
                 res[b_key]["ratio"],
-                yerr=b_err,
-                fmt=marker,
+                marker,
                 color=color_B,
                 label=f"{b_name} {pair}",
                 markersize=5,
                 alpha=0.8,
-                capsize=0,
             )
             ax.plot(
                 xi,
@@ -212,7 +197,6 @@ def lin_shift(x, i):
             ax,
             r"$\theta$ [arcmin]",
             rf"${name}$: covariance ratio across n(z)",
-            show_mc_band=True,
         )
         plot_ratios(
             ax,
@@ -223,7 +207,6 @@ def lin_shift(x, i):
             "B-mode",
             "E-mode",
             log_shift,
-            ratio_err,
         )
         ax.legend(loc="upper right", fontsize=7, ncol=2)
         ax.set_xlim(1, 300)
@@ -237,17 +220,15 @@ def lin_shift(x, i):
     ax.axhline(1.0, color="gray", ls="-", lw=0.8, zorder=0)
     ax.set_xlabel(r"Mode $n$")
     ax.set_title(r"COSEBIS: covariance ratio across n(z)")
-    plot_ratios(
-        ax, n_arr, cosebis_results, "B", "E", "B-mode", "E-mode", lin_shift, None
-    )
+    plot_ratios(ax, n_arr, cosebis_results, "B", "E", "B-mode", "E-mode", lin_shift)
     ax.set_ylabel("Diagonal ratio")
     ax.legend(loc="upper right", fontsize=7, ncol=2)
     ax.set_ylim(0.85, 1.15)
 
     # --- Panel 4: C_ell^BB vs C_ell^EE ---
     ax = axes[1, 1]
     setup_ratio_panel(ax, r"$\ell$", r"$C_\ell$: covariance ratio across n(z)")
-    plot_ratios(ax, ell_eff, harmonic_results, "BB", "EE", "BB", "EE", log_shift, None)
+    plot_ratios(ax, ell_eff, harmonic_results, "BB", "EE", "BB", "EE", log_shift)
     ax.set_ylabel("Diagonal ratio")
     ax.legend(loc="upper right", fontsize=7, ncol=2)
     ax.set_ylim(0.85, 1.15)
@@ -330,7 +311,6 @@ def ratios_to_reference(data, modes):
     for res in cosebis_results.values():
         res["nmodes"] = nmodes
 
-    n_samples = config["covariance"]["n_samples"]
     make_figure(
         theta,
         ell_eff,
@@ -339,7 +319,6 @@ def ratios_to_reference(data, modes):
         cosebis_results,
         reference,
         figure_path,
-        n_samples=n_samples,
     )
 
     def max_dev(results, mode):
@@ -485,7 +464,7 @@ def _from_cli(argv=None):
     ap.add_argument(
         "--pure-eb-dir",
         required=True,
-        help="Dir with {version}_pure_eb_semianalytic.npz per n(z) realisation",
+        help="Dir with {version}_pure_eb.npz per n(z) realisation",
     )
     ap.add_argument(
         "--cosmo-val-dir",
@@ -517,7 +496,7 @@ def _from_cli(argv=None):
 
     realisations = fid["nz_realisations"]
     pure_eb_paths = {
-        b: os.path.join(a.pure_eb_dir, f"{ver}_pure_eb_semianalytic.npz")
+        b: os.path.join(a.pure_eb_dir, f"{ver}_pure_eb.npz")
         for b, ver in realisations.items()
     }
     harmonic_paths = {

papers/bmodes/scripts/calculate_pure_eb_ptes.py

@@ -1,16 +1,14 @@
 """Calculate PTE matrices for pure E/B-mode scale-cut robustness.
 
-CLI refactor of the former Snakemake ``script:`` rule. Reads the gathered
-pure-E/B ``semianalytic.npz`` (data vectors + MC covariance), evaluates the
+CLI refactor of the former Snakemake ``script:`` rule. Reads the pure-E/B
+``_pure_eb.npz`` (data vectors + analytic covariance), evaluates the
 ξ_+^B / ξ_-^B / joint ξ_tot^B χ² PTE matrices over the scale-cut grid via
-``sp_validation.b_modes.calculate_eb_statistics`` (Hartlap-corrected inverse
-MC covariance, debiased by the draw count), and writes the PTE matrices to
-``{out}/{version}_pure_eb_ptes.npz``.
+``sp_validation.b_modes.calculate_eb_statistics``, and writes the PTE matrices
+to ``{out}/{version}_pure_eb_ptes.npz``.
 
     python calculate_pure_eb_ptes.py \
         --version SP_v1.4.6.3_leak_corr \
-        --pure-eb-data <..._pure_eb_semianalytic.npz> \
-        --n-samples 2000 --out <output_dir>
+        --pure-eb-data <..._pure_eb.npz> --out <output_dir>
 """
 
 import argparse
@@ -24,7 +22,6 @@
 def calculate_ptes(
     version,
     pure_eb_data,
-    n_samples,
     output_dir,
 ):
     dataset = np.load(pure_eb_data)
@@ -33,8 +30,8 @@ def calculate_ptes(
 
     results = {
         "theta": theta,
-        # The MC draws are the realisations behind this covariance.
-        "n_eff": int(n_samples),
+        # The covariance is analytic: no Hartlap factor.
+        "npatch": None,
         "xip_E": dataset["xip_E"],
         "xim_E": dataset["xim_E"],
         "xip_B": dataset["xip_B"],
@@ -50,6 +47,9 @@ def calculate_ptes(
     pte_matrices = results["pte_matrices"]
     output_data = {
         "theta": theta,
+        # The bins the matrices are indexed on, for scale-cut windows.
+        "left_edges": dataset["left_edges"],
+        "right_edges": dataset["right_edges"],
         "pte_xip_B": pte_matrices["xip_B"],
         "pte_xim_B": pte_matrices["xim_B"],
         "pte_combined": pte_matrices["combined"],
@@ -65,14 +65,12 @@ def calculate_ptes(
 def _from_cli(argv=None):
     ap = argparse.ArgumentParser(description=__doc__.split("\n")[0])
     ap.add_argument("--version", required=True)
-    ap.add_argument("--pure-eb-data", required=True, help="Gathered semianalytic .npz")
-    ap.add_argument("--n-samples", type=int, default=2000)
+    ap.add_argument("--pure-eb-data", required=True, help="Pure E/B .npz")
     ap.add_argument("--out", required=True, help="Output directory (lc {output})")
     a = ap.parse_args(argv)
     calculate_ptes(
         version=a.version,
         pure_eb_data=a.pure_eb_data,
-        n_samples=a.n_samples,
         output_dir=a.out,
     )