diff --git a/papers/bmodes/config/config.yaml b/papers/bmodes/config/config.yaml index 1b31c647..1198d984 100644 --- a/papers/bmodes/config/config.yaml +++ b/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) diff --git a/papers/bmodes/rules/figures.smk b/papers/bmodes/rules/figures.smk index 590735e3..31b65ea0 100644 --- a/papers/bmodes/rules/figures.smk +++ b/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), diff --git a/papers/bmodes/scripts/bb_covariance_nz_independence.py b/papers/bmodes/scripts/bb_covariance_nz_independence.py index 69c0db18..7f8cae02 100644 --- a/papers/bmodes/scripts/bb_covariance_nz_independence.py +++ b/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,9 +220,7 @@ 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) @@ -247,7 +228,7 @@ def lin_shift(x, i): # --- 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 = { diff --git a/papers/bmodes/scripts/calculate_pure_eb_ptes.py b/papers/bmodes/scripts/calculate_pure_eb_ptes.py index feb0be78..ca453dc4 100644 --- a/papers/bmodes/scripts/calculate_pure_eb_ptes.py +++ b/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 + --pure-eb-data <..._pure_eb.npz> --out """ 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, ) diff --git a/papers/bmodes/scripts/config_space_pte_matrices.py b/papers/bmodes/scripts/config_space_pte_matrices.py index 3bfa2811..f4365ec1 100644 --- a/papers/bmodes/scripts/config_space_pte_matrices.py +++ b/papers/bmodes/scripts/config_space_pte_matrices.py @@ -39,15 +39,19 @@ make_pte_norm, ) +from sp_validation.b_modes import bins_from_scale_cut + plt.style.use(PAPER_MPLSTYLE) def resolve_fiducial_bin_window(edges, theta_min, theta_max): - """Return the first and last reporting bins inside a scale-cut window.""" - left, right = edges[:-1], edges[1:] - inside = (left >= theta_min * (1.0 - 1e-2)) & (right <= theta_max * (1.0 + 1e-2)) - bins = np.flatnonzero(inside) - return int(bins[0]), int(bins[-1]) + """Return the first and last reporting bins of a pure-E/B scale cut. + + The cut snaps to the nearest reporting edges, as everywhere pure-E/B PTEs + are read (``sp_validation.b_modes.bins_from_scale_cut``). + """ + start, stop = bins_from_scale_cut(edges[:-1], edges[1:], (theta_min, theta_max)) + return start, stop - 1 def _path_matches_version(path, version): @@ -189,22 +193,25 @@ def load_pure_eb_pte_matrices(pte_files, version, override_path=None): Angular scale grid. pte_combined : ndarray or None PTE matrix for combined ξ_tot^B, or None if not available. + edges : ndarray or None + The reporting edges the matrices are indexed on, or None for a file + that does not record them. """ - if override_path is not None: - data = np.load(override_path) - pte_combined = data["pte_combined"] if "pte_combined" in data else None - return data["pte_xip_B"], data["pte_xim_B"], data["theta"], pte_combined - - for pte_file in pte_files: + if override_path is None: # Filter to this version (exact match, no substring false positives) - if not _path_matches_version(pte_file, version): - continue - - data = np.load(pte_file) - pte_combined = data["pte_combined"] if "pte_combined" in data else None - return data["pte_xip_B"], data["pte_xim_B"], data["theta"], pte_combined - - raise ValueError(f"No PTE file found for version {version}") + matching = [p for p in pte_files if _path_matches_version(p, version)] + if not matching: + raise ValueError(f"No PTE file found for version {version}") + override_path = matching[0] + + data = np.load(override_path) + pte_combined = data["pte_combined"] if "pte_combined" in data else None + edges = ( + np.append(data["left_edges"], data["right_edges"][-1]) + if "left_edges" in data + else None + ) + return data["pte_xip_B"], data["pte_xim_B"], data["theta"], pte_combined, edges def _load_version_pte_data( @@ -215,12 +222,19 @@ def _load_version_pte_data( Returns ------- dict with keys: pte_xip_B, pte_xim_B, pte_combined (or None), - pte_cosebis, pte_cosebis_20, theta_pure_eb, theta_cosebis. + pte_cosebis, pte_cosebis_20, theta_pure_eb, edges_pure_eb, + theta_cosebis. """ pure_eb_override, cosebis_override = _resolve_overrides(version, fiducial_overrides) - pte_xip_B, pte_xim_B, theta_pure_eb, pte_combined = load_pure_eb_pte_matrices( - pure_eb_pte_files, version, override_path=pure_eb_override + pte_xip_B, pte_xim_B, theta_pure_eb, pte_combined, edges_pure_eb = ( + load_pure_eb_pte_matrices( + pure_eb_pte_files, version, override_path=pure_eb_override + ) ) + if edges_pure_eb is None: + # A PTE file without saved edges was binned on the nominal grid. + fid = config["fiducial"] + edges_pure_eb = np.geomspace(fid["min_sep"], fid["max_sep"], fid["nbins"] + 1) pte_cosebis, theta_cosebis = load_cosebis_pte_matrix( cosebis_pte_files, version, @@ -242,6 +256,7 @@ def _load_version_pte_data( "pte_cosebis": pte_cosebis, "pte_cosebis_20": pte_cosebis_20, "theta_pure_eb": theta_pure_eb, + "edges_pure_eb": edges_pure_eb, "theta_cosebis": theta_cosebis, } @@ -522,13 +537,9 @@ def create_3panel_composite( cosebis_fid_start = np.argmin(np.abs(theta_cosebis[:-1] - cosebis_fid[0])) cosebis_fid_stop = np.argmin(np.abs(theta_cosebis[1:] - cosebis_fid[1])) + 1 - reporting_edges = np.geomspace( - config["fiducial"]["min_sep"], - config["fiducial"]["max_sep"], - config["fiducial"]["nbins"] + 1, - ) - xip_start, xip_stop = resolve_fiducial_bin_window(reporting_edges, *xip_fid) - xim_start, xim_stop = resolve_fiducial_bin_window(reporting_edges, *xim_fid) + edges_pure_eb = matrices["edges_pure_eb"] + xip_start, xip_stop = resolve_fiducial_bin_window(edges_pure_eb, *xip_fid) + xim_start, xim_stop = resolve_fiducial_bin_window(edges_pure_eb, *xim_fid) # Create subplot axes ax_xip = fig.add_subplot(gs[0, 0]) @@ -680,13 +691,9 @@ def create_9panel_composite( cosebis_fid_start = np.argmin(np.abs(theta_cosebis[:-1] - cosebis_fid[0])) cosebis_fid_stop = np.argmin(np.abs(theta_cosebis[1:] - cosebis_fid[1])) + 1 - reporting_edges = np.geomspace( - config["fiducial"]["min_sep"], - config["fiducial"]["max_sep"], - config["fiducial"]["nbins"] + 1, - ) - xip_start, xip_stop = resolve_fiducial_bin_window(reporting_edges, *xip_fid) - xim_start, xim_stop = resolve_fiducial_bin_window(reporting_edges, *xim_fid) + edges_pure_eb = matrices["edges_pure_eb"] + xip_start, xip_stop = resolve_fiducial_bin_window(edges_pure_eb, *xip_fid) + xim_start, xim_stop = resolve_fiducial_bin_window(edges_pure_eb, *xim_fid) # Create subplot axes for this row ax_xip = fig.add_subplot(gs[row_idx, 0]) @@ -910,16 +917,12 @@ def main( fiducial_overrides, ) theta_co = matrices["theta_cosebis"] - reporting_edges = np.geomspace( - config["fiducial"]["min_sep"], - config["fiducial"]["max_sep"], - config["fiducial"]["nbins"] + 1, - ) + edges_pure_eb = matrices["edges_pure_eb"] xip_start, xip_stop = resolve_fiducial_bin_window( - reporting_edges, *xip_fid + edges_pure_eb, *xip_fid ) xim_start, xim_stop = resolve_fiducial_bin_window( - reporting_edges, *xim_fid + edges_pure_eb, *xim_fid ) cos_start = np.argmin(np.abs(theta_co[:-1] - cosebis_fid[0])) cos_stop = np.argmin(np.abs(theta_co[1:] - cosebis_fid[1])) + 1 diff --git a/papers/bmodes/scripts/gather_pure_eb_chunks.py b/papers/bmodes/scripts/gather_pure_eb_chunks.py deleted file mode 100644 index b6d5d702..00000000 --- a/papers/bmodes/scripts/gather_pure_eb_chunks.py +++ /dev/null @@ -1,136 +0,0 @@ -"""Gather MC sample chunks and compute the final pure E/B covariance. - -CLI refactor of the former Snakemake ``script:`` gather rule. Reads the actual -ξ± data vectors (reporting + integration grids), computes the pure E/B/ambiguous -decomposition (Schneider 2022), stacks the per-chunk MC sample blocks, and -forms the empirical 6-block covariance. Writes the per-version -``_pure_eb_semianalytic.npz`` consumed by every downstream -pure-mode plot / PTE. - - python gather_pure_eb_chunks.py \ - --version SP_v1.4.6.3_leak_corr \ - --xi-reporting \ - --xi-integration \ - --chunks-dir \ - --min-sep 1.0 --max-sep 250.0 --nbins 20 \ - --min-sep-int 0.5 --max-sep-int 300.0 --nbins-int 1000 \ - --out -""" - -import argparse -import glob -import os - -import numpy as np - - -def _load_xi(path, nbins): - """Load ξ± from a TreeCorr text dump.""" - data = np.loadtxt(path, comments="#", max_rows=nbins) - return {"meanr": data[:, 1], "xip": data[:, 3], "xim": data[:, 4]} - - -def gather( - version, - xi_reporting, - xi_integration, - chunk_files, - min_sep, - max_sep, - nbins, - nbins_int, - output_dir, -): - from cosmo_numba.B_modes.schneider2022 import get_pure_EB_modes - - print(f"Gathering pure E/B for {version}") - - gg = _load_xi(xi_reporting, nbins) - gg_int = _load_xi(xi_integration, nbins_int) - - eb_results = get_pure_EB_modes( - theta=gg["meanr"], - xip=gg["xip"], - xim=gg["xim"], - theta_int=gg_int["meanr"], - xip_int=gg_int["xip"], - xim_int=gg_int["xim"], - tmin=min_sep, - tmax=max_sep, - ) - xip_E, xim_E, xip_B, xim_B, xip_amb, xim_amb = eb_results - - chunk_files = sorted(chunk_files) - all_samples = [] - for chunk_file in chunk_files: - data = np.load(chunk_file) - all_samples.append(data["eb_samples"]) - print(f"Loaded {len(data['eb_samples'])} samples from {chunk_file}") - - eb_samples = np.vstack(all_samples) - print(f"Total samples: {len(eb_samples)}") - - cov_pure_eb = np.cov(eb_samples.T) - - package = { - "theta": gg["meanr"], - "theta_int": gg_int["meanr"], - "xip_total": gg["xip"], - "xim_total": gg["xim"], - "xip_E": xip_E, - "xim_E": xim_E, - "xip_B": xip_B, - "xim_B": xim_B, - "xip_amb": xip_amb, - "xim_amb": xim_amb, - "cov_pure_eb": cov_pure_eb, - } - - os.makedirs(output_dir, exist_ok=True) - out_path = os.path.join(output_dir, f"{version}_pure_eb_semianalytic.npz") - np.savez(out_path, **package) - print(f"Saved to {out_path}") - return out_path - - -def _resolve_chunks(args): - if args.chunks: - files = args.chunks - else: - files = glob.glob(os.path.join(args.chunks_dir, "pure_eb_chunk_*.npz")) - if not files: - raise SystemExit(f"No chunk .npz found (chunks-dir={args.chunks_dir})") - return files - - -def _from_cli(argv=None): - ap = argparse.ArgumentParser(description=__doc__.split("\n")[0]) - ap.add_argument("--version", required=True) - ap.add_argument("--xi-reporting", required=True) - ap.add_argument("--xi-integration", required=True) - ap.add_argument("--chunks-dir", help="Directory holding pure_eb_chunk_*.npz") - ap.add_argument("--chunks", nargs="+", help="Explicit chunk .npz paths") - ap.add_argument("--min-sep", type=float, default=1.0) - ap.add_argument("--max-sep", type=float, default=250.0) - ap.add_argument("--nbins", type=int, default=20) - ap.add_argument("--min-sep-int", type=float, default=0.5) - ap.add_argument("--max-sep-int", type=float, default=300.0) - ap.add_argument("--nbins-int", type=int, default=1000) - ap.add_argument("--npatch", type=int, default=1) - ap.add_argument("--out", required=True, help="Output directory (lc {output})") - a = ap.parse_args(argv) - gather( - version=a.version, - xi_reporting=a.xi_reporting, - xi_integration=a.xi_integration, - chunk_files=_resolve_chunks(a), - min_sep=a.min_sep, - max_sep=a.max_sep, - nbins=a.nbins, - nbins_int=a.nbins_int, - output_dir=a.out, - ) - - -if __name__ == "__main__": - _from_cli() diff --git a/papers/bmodes/scripts/precompute_pure_eb_chunk.py b/papers/bmodes/scripts/precompute_pure_eb_chunk.py deleted file mode 100644 index 2db6bf92..00000000 --- a/papers/bmodes/scripts/precompute_pure_eb_chunk.py +++ /dev/null @@ -1,209 +0,0 @@ -"""Compute one chunk of MC samples for the pure E/B covariance. - -CLI refactor of the former Snakemake ``script:`` rule. The compute is -unchanged: draw ``n_samples // n_chunks`` Gaussian realisations of ξ±(θ) from -the 1000-bin integration-grid CosmoCov covariance (deterministic seed -``42 + chunk_id``), rebin to the reporting grid, and push each draw through the -Schneider-2022 pure-mode integral transforms (``cosmo_numba``). The per-chunk -E/B/amb sample block is written to ``{out}/pure_eb_chunk_{chunk_id}.npz`` for -the gather stage. Each chunk is independent (fresh RNG per chunk_id), so the -20 chunks reproduce the paper's 2000-sample covariance bit-for-bit whether run -in parallel or looped in one process. - - python precompute_pure_eb_chunk.py \ - --chunk-id 0 --n-chunks 20 --n-samples 2000 \ - --version SP_v1.4.6.3_leak_corr \ - --cat-config /path/cosmo_val/cat_config.yaml \ - --xi-reporting \ - --xi-integration \ - --cov-integration \ - --min-sep 1.0 --max-sep 250.0 --nbins 20 \ - --min-sep-int 0.5 --max-sep-int 300.0 --nbins-int 1000 \ - --npatch 1 --out -""" - -import argparse -import os - -import numpy as np -import tqdm -from scipy import sparse - - -def _build_cosmology(cosmo_params): - """Build a CCL cosmology from a PLANCK18-style params dict.""" - import pyccl as ccl - - return ccl.Cosmology( - Omega_c=cosmo_params["Omega_m"] - cosmo_params["Omega_b"], - Omega_b=cosmo_params["Omega_b"], - h=cosmo_params["h"], - sigma8=cosmo_params["sigma_8"], - n_s=cosmo_params["n_s"], - ) - - -def _load_xi(path, min_sep, max_sep, nbins): - """Load ξ± from a TreeCorr text dump and recompute the log bin edges.""" - data = np.loadtxt(path, comments="#", max_rows=nbins) - meanr = data[:, 1] - xip = data[:, 3] - xim = data[:, 4] - bin_edges = np.logspace(np.log10(min_sep), np.log10(max_sep), nbins + 1) - return { - "meanr": meanr, - "xip": xip, - "xim": xim, - "left_edges": bin_edges[:-1], - "right_edges": bin_edges[1:], - } - - -def compute_chunk( - chunk_id, - n_chunks, - n_samples_total, - version, - cat_config, - xi_reporting, - xi_integration, - cov_integration, - min_sep, - max_sep, - nbins, - min_sep_int, - max_sep_int, - nbins_int, - output_dir, - cosmo_params=None, -): - from cosmo_numba.B_modes.schneider2022 import get_pure_EB_modes - from cs_util.cosmo import PLANCK18, get_theo_xi - - from sp_validation.cosmo_val import CosmologyValidation - - if cosmo_params is None: - cosmo_params = dict(PLANCK18) - - samples_per_chunk = n_samples_total // n_chunks - start_idx = chunk_id * samples_per_chunk - end_idx = ( - start_idx + samples_per_chunk if chunk_id < n_chunks - 1 else n_samples_total - ) - n_samples_chunk = end_idx - start_idx - print( - f"Chunk {chunk_id}/{n_chunks}: samples {start_idx}-{end_idx} " - f"({n_samples_chunk} samples)" - ) - - gg = _load_xi(xi_reporting, min_sep, max_sep, nbins) - gg_int = _load_xi(xi_integration, min_sep_int, max_sep_int, nbins_int) - - cv = CosmologyValidation( - versions=[version], - catalog_config=cat_config, - output_dir=output_dir, - ) - z, nz = cv.get_redshift(version) - z_dist = np.column_stack([z, nz]) - - cosmo_cov = _build_cosmology(cosmo_params) - - cov_int = np.loadtxt(cov_integration) - - theta_int = gg_int["meanr"] - reporting_bin_edges = np.concatenate([gg["left_edges"], [gg["right_edges"][-1]]]) - bin_indices = np.digitize(theta_int, reporting_bin_edges) - 1 - valid_mask = (bin_indices >= 0) & (bin_indices < len(gg["meanr"])) - row_indices, col_indices = bin_indices[valid_mask], np.where(valid_mask)[0] - - binning_matrix = sparse.csr_matrix( - (np.ones(len(row_indices)), (row_indices, col_indices)), - shape=(len(gg["meanr"]), nbins_int), - ) - row_sums = np.array(binning_matrix.sum(axis=1)).flatten() - binning_matrix = sparse.diags(1 / row_sums) @ binning_matrix - - # One n(z) gives one tracer pair: get_theo_xi's single (xi+, xi-) entry. - (xi_pm,) = get_theo_xi( - theta=theta_int, - z=z_dist[:, 0], - nz=z_dist[:, 1], - backend="ccl", - cosmo=cosmo_cov, - ).values() - mean_int = np.concatenate(xi_pm) - - rng = np.random.default_rng(seed=42 + chunk_id) - - samples_int = rng.multivariate_normal(mean_int, cov_int, size=n_samples_chunk) - samples_int_xip = samples_int[:, :nbins_int] - samples_int_xim = samples_int[:, nbins_int:] - samples_rep_xip = (binning_matrix @ samples_int_xip.T).T - samples_rep_xim = (binning_matrix @ samples_int_xim.T).T - - transformed_samples = [ - np.concatenate( - get_pure_EB_modes( - theta=gg["meanr"], - theta_int=gg_int["meanr"], - xip=samples_rep_xip[i], - xim=samples_rep_xim[i], - xip_int=samples_int_xip[i], - xim_int=samples_int_xim[i], - tmin=min_sep, - tmax=max_sep, - ) - ) - for i in tqdm.tqdm(range(n_samples_chunk), desc=f"Chunk {chunk_id}") - ] - - eb_samples = np.array(transformed_samples) - - os.makedirs(output_dir, exist_ok=True) - out_path = os.path.join(output_dir, f"pure_eb_chunk_{chunk_id}.npz") - np.savez(out_path, eb_samples=eb_samples, chunk_id=chunk_id) - print(f"Saved {n_samples_chunk} samples to {out_path}") - return out_path - - -def _from_cli(argv=None): - ap = argparse.ArgumentParser(description=__doc__.split("\n")[0]) - ap.add_argument("--chunk-id", type=int, required=True) - ap.add_argument("--n-chunks", type=int, default=20) - ap.add_argument("--n-samples", type=int, default=2000) - ap.add_argument("--version", required=True) - ap.add_argument("--cat-config", required=True) - ap.add_argument("--xi-reporting", required=True) - ap.add_argument("--xi-integration", required=True) - ap.add_argument("--cov-integration", required=True) - ap.add_argument("--min-sep", type=float, default=1.0) - ap.add_argument("--max-sep", type=float, default=250.0) - ap.add_argument("--nbins", type=int, default=20) - ap.add_argument("--min-sep-int", type=float, default=0.5) - ap.add_argument("--max-sep-int", type=float, default=300.0) - ap.add_argument("--nbins-int", type=int, default=1000) - ap.add_argument("--npatch", type=int, default=1) - ap.add_argument("--out", required=True, help="Output directory (lc {output})") - a = ap.parse_args(argv) - compute_chunk( - chunk_id=a.chunk_id, - n_chunks=a.n_chunks, - n_samples_total=a.n_samples, - version=a.version, - cat_config=a.cat_config, - xi_reporting=a.xi_reporting, - xi_integration=a.xi_integration, - cov_integration=a.cov_integration, - min_sep=a.min_sep, - max_sep=a.max_sep, - nbins=a.nbins, - min_sep_int=a.min_sep_int, - max_sep_int=a.max_sep_int, - nbins_int=a.nbins_int, - output_dir=a.out, - ) - - -if __name__ == "__main__": - _from_cli() diff --git a/papers/bmodes/scripts/pure_eb_covariance.py b/papers/bmodes/scripts/pure_eb_covariance.py index cb137c62..08898968 100644 --- a/papers/bmodes/scripts/pure_eb_covariance.py +++ b/papers/bmodes/scripts/pure_eb_covariance.py @@ -226,8 +226,7 @@ def _from_cli(argv=None): ap.add_argument( "--pure-eb-data", required=True, - help="Fiducial _pure_eb_semianalytic.npz " - "(provides the 6-block cov_pure_eb)", + help="Fiducial _pure_eb.npz (provides the 6-block cov_pure_eb)", ) ap.add_argument("--out", required=True, help="Output directory (lc {output})") a = ap.parse_args(argv) diff --git a/papers/bmodes/scripts/pure_eb_data_vector.py b/papers/bmodes/scripts/pure_eb_data_vector.py index d7845361..d79060cc 100644 --- a/papers/bmodes/scripts/pure_eb_data_vector.py +++ b/papers/bmodes/scripts/pure_eb_data_vector.py @@ -7,7 +7,7 @@ CLI: python pure_eb_data_vector.py \ --config config.yaml \ - --pure-eb-data _pure_eb_semianalytic.npz \ + --pure-eb-data _pure_eb.npz \ --reporting-cov /covariance_processed.txt \ --out """ @@ -34,7 +34,7 @@ def _extract_sigma(covariance, block_index, block_size): return np.sqrt(np.clip(np.diag(covariance[block_slice, block_slice]), 0, None)) -def _compute_joint_pte(xip_B, xim_B, cov_xip_B, cov_xim_B, cov_cross, n_samples=None): +def _compute_joint_pte(xip_B, xim_B, cov_xip_B, cov_xim_B, cov_cross): """Compute joint PTE for combined B-mode data vector [xip_B, xim_B].""" data_joint = np.concatenate([xip_B, xim_B]) n_xip, n_xim = len(xip_B), len(xim_B) @@ -45,12 +45,12 @@ def _compute_joint_pte(xip_B, xim_B, cov_xip_B, cov_xim_B, cov_cross, n_samples= cov_joint[:n_xip, n_xip:] = cov_cross cov_joint[n_xip:, :n_xip] = cov_cross.T - chi2, pte, dof = compute_chi2_pte(data_joint, cov_joint, n_samples=n_samples) + chi2, pte, dof = compute_chi2_pte(data_joint, cov_joint) return pte, chi2, dof def _load_pure_eb_data(pure_eb_path, cov_path): - """Load pure E/B decomposition, the 6-block MC covariance, and the + """Load pure E/B decomposition, the 6-block covariance, and the reporting-grid CosmoCov ξ± covariance used for the total-curve error bars.""" dataset = np.load(pure_eb_path) theta = dataset["theta"] @@ -241,9 +241,6 @@ def main(config, pure_eb_path, cov_path, out_dir): cov_pure_eb = data["cov_pure_eb"] xip_B, xim_B = data["xip_B"], data["xim_B"] - # Hartlap correction: MC-propagated covariance uses n_samples from config - n_samples = int(config["covariance"]["n_samples"]) - # Extract B-mode covariance blocks cov_xip_B_full = cov_pure_eb[2 * nbins : 3 * nbins, 2 * nbins : 3 * nbins] cov_xim_B_full = cov_pure_eb[3 * nbins : 4 * nbins, 3 * nbins : 4 * nbins] @@ -264,10 +261,10 @@ def main(config, pure_eb_path, cov_path, out_dir): # Compute PTEs at fiducial scale cuts chi2_xip_fid, pte_xip_fid, dof_xip_fid = compute_chi2_pte( - xip_B[mask_xip], cov_xip_B_cut, n_samples=n_samples + xip_B[mask_xip], cov_xip_B_cut ) chi2_xim_fid, pte_xim_fid, dof_xim_fid = compute_chi2_pte( - xim_B[mask_xim], cov_xim_B_cut, n_samples=n_samples + xim_B[mask_xim], cov_xim_B_cut ) pte_joint_fid, chi2_joint_fid, dof_joint_fid = _compute_joint_pte( xip_B[mask_xip], @@ -275,19 +272,17 @@ def main(config, pure_eb_path, cov_path, out_dir): cov_xip_B_cut, cov_xim_B_cut, cov_cross_cut, - n_samples=n_samples, ) # Compute PTEs at full range - _, pte_xip_full, _ = compute_chi2_pte(xip_B, cov_xip_B_full, n_samples=n_samples) - _, pte_xim_full, _ = compute_chi2_pte(xim_B, cov_xim_B_full, n_samples=n_samples) + _, pte_xip_full, _ = compute_chi2_pte(xip_B, cov_xip_B_full) + _, pte_xim_full, _ = compute_chi2_pte(xim_B, cov_xim_B_full) pte_joint_full, chi2_joint_full, dof_joint_full = _compute_joint_pte( xip_B, xim_B, cov_xip_B_full, cov_xim_B_full, cov_cross_full, - n_samples=n_samples, ) print( @@ -338,8 +333,7 @@ def _from_cli(argv=None): ap.add_argument( "--pure-eb-data", required=True, - help="Fiducial _pure_eb_semianalytic.npz " - "(decomposed ξ± + 6-block MC covariance)", + help="Fiducial _pure_eb.npz (decomposed ξ± + 6-block covariance)", ) ap.add_argument( "--reporting-cov", diff --git a/papers/bmodes/scripts/pure_eb_modes.py b/papers/bmodes/scripts/pure_eb_modes.py new file mode 100644 index 00000000..643fe12c --- /dev/null +++ b/papers/bmodes/scripts/pure_eb_modes.py @@ -0,0 +1,116 @@ +"""Pure E/B modes and their exact covariance for one version. + +Reads the fine-grid ξ± (TreeCorr text dump, with its pair weights) and the +Gaussian ξ± covariance on the same grid, and applies +``sp_validation.b_modes.calculate_pure_eb_correlation``: the fixed-operator +pure-E/B estimator averaged into the reporting bins, with covariance +``K C_ξ Kᵀ``. Writes ``_pure_eb.npz``, the input of every downstream +pure-mode plot / PTE. + + python pure_eb_modes.py \ + --xi-integration \ + --cov-integration \ + --min-sep 1.0 --max-sep 250.0 --nbins 20 \ + --min-sep-int 0.5 --max-sep-int 300.0 --nbins-int 1000 \ + --out _pure_eb.npz +""" + +import argparse +import os + +import numpy as np + +from sp_validation.b_modes import calculate_pure_eb_correlation +from sp_validation.sacc_io import PURE_KEYS + + +def _load_xi(path, nbins): + """``(meanr, xip, xim, weight)`` from a TreeCorr text dump. + + TreeCorr's ASCII header is + r_nom meanr meanlogr xip xim xip_im xim_im sigma_xip sigma_xim weight npairs. + """ + data = np.loadtxt(path, comments="#", max_rows=nbins) + return data[:, 1], data[:, 3], data[:, 4], data[:, 9] + + +def pure_eb_modes( + xi_integration, + cov_integration, + min_sep, + max_sep, + nbins, + min_sep_int, + max_sep_int, + nbins_int, + out_path, +): + results = calculate_pure_eb_correlation( + *_load_xi(xi_integration, nbins_int), + np.geomspace(min_sep_int, max_sep_int, nbins_int + 1), + np.loadtxt(cov_integration), + np.geomspace(min_sep, max_sep, nbins + 1), + ) + package = { + "theta": results["theta"], + "left_edges": results["left_edges"], + "right_edges": results["right_edges"], + "theta_int": results["theta_int"], + "xip_total": results["xip"], + "xim_total": results["xim"], + **{key: results[key] for key in PURE_KEYS}, + "cov_pure_eb": results["cov"], + } + os.makedirs(os.path.dirname(out_path) or ".", exist_ok=True) + np.savez(out_path, **package) + print(f"Saved pure E/B to {out_path}") + return out_path + + +def _from_snakemake(smk): + p = smk.params + pure_eb_modes( + smk.input["xi_integration"], + smk.input["cov_integration"], + p["min_sep"], + p["max_sep"], + p["nbins"], + p["min_sep_int"], + p["max_sep_int"], + p["nbins_int"], + smk.output[0], + ) + + +def _from_cli(argv=None): + ap = argparse.ArgumentParser(description=__doc__.split("\n")[0]) + ap.add_argument("--xi-integration", required=True) + ap.add_argument("--cov-integration", required=True) + ap.add_argument("--min-sep", type=float, default=1.0) + ap.add_argument("--max-sep", type=float, default=250.0) + ap.add_argument("--nbins", type=int, default=20) + ap.add_argument("--min-sep-int", type=float, default=0.5) + ap.add_argument("--max-sep-int", type=float, default=300.0) + ap.add_argument("--nbins-int", type=int, default=1000) + ap.add_argument("--out", required=True, help="Output .npz path") + a = ap.parse_args(argv) + pure_eb_modes( + a.xi_integration, + a.cov_integration, + a.min_sep, + a.max_sep, + a.nbins, + a.min_sep_int, + a.max_sep_int, + a.nbins_int, + a.out, + ) + + +if __name__ == "__main__": + try: + snakemake # noqa: F821 — injected by Snakemake's script: directive + except NameError: + _from_cli() + else: + _from_snakemake(snakemake) # noqa: F821 diff --git a/papers/bmodes/scripts/pure_eb_version_comparison.py b/papers/bmodes/scripts/pure_eb_version_comparison.py index 3e99e68b..6f20c054 100644 --- a/papers/bmodes/scripts/pure_eb_version_comparison.py +++ b/papers/bmodes/scripts/pure_eb_version_comparison.py @@ -243,7 +243,7 @@ def _create_version_comparison_figure( def _pure_eb_npz(results_dir, ver): - return f"{results_dir}/{ver}_pure_eb_semianalytic.npz" + return f"{results_dir}/{ver}_pure_eb.npz" def main( @@ -472,7 +472,7 @@ def _from_cli(argv=None): ap.add_argument( "--results-dir", required=True, - help="Directory holding per-version _pure_eb_semianalytic.npz files", + help="Directory holding per-version _pure_eb.npz files", ) ap.add_argument("--out", required=True, help="Output directory (lc {output})") ap.add_argument( diff --git a/papers/cosmo_val/config/config.yaml b/papers/cosmo_val/config/config.yaml index 8501bfbc..5b33999f 100644 --- a/papers/cosmo_val/config/config.yaml +++ b/papers/cosmo_val/config/config.yaml @@ -57,7 +57,7 @@ cosmo_val: # C_l^BB column belongs in the standard B-mode summary. include_pseudo_cl: true - # Theory cosmology for pseudo-Cl / semi-analytic covariance (astropy Planck18 + # Theory cosmology for the pseudo-Cl covariance (astropy Planck18 # + CAMB nonlinear settings, matching the original run_cosmo_val.py driver). cosmo_params: Omega_m: 0.30966 @@ -110,8 +110,6 @@ var_method: jackknife covariance: default_masked: true - use_semianalytic: true - n_samples: 2000 cl: n_ell_bins: 32 diff --git a/pyproject.toml b/pyproject.toml index 35c9e6ce..2ffcd753 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -38,13 +38,13 @@ dependencies = [ # get_theo_c_ell / get_theo_xi / PLANCK18 in cs_util.cosmo, which land via # CosmoStat/cs_util#76 — so this goes green once #76 merges into develop. "cs_util @ git+https://github.com/CosmoStat/cs_util.git@develop", - # Fast numba B-mode kernels (Schneider et al. 2022): the Schneider E/B split - # and COSEBIS live here, imported in b_modes.py. Tracks aguinot/cosmo-numba - # main (not published on PyPI). main carries the numpy-2 FFT fix via its - # rocket-fft dependency (which teaches numba's nopython mode to handle - # np.fft), and declares numba/numpy/rocket-fft from its requirements.txt so - # those constraints reach the resolver. - "cosmo-numba @ git+https://github.com/aguinot/cosmo-numba.git@main", + # Fast numba B-mode kernels (Schneider et al. 2022): the fixed-quadrature + # pure-E/B operator (schneider2022) and COSEBIS live here, imported + # in b_modes.py. Pinned to a commit on cailmdaley/cosmo-numba, which carries + # the operator on top of aguinot/cosmo-numba main (not published on PyPI): + # main's numpy-2 FFT fix via rocket-fft, and its numba/numpy/rocket-fft + # requirements, which reach the resolver. + "cosmo-numba @ git+https://github.com/cailmdaley/cosmo-numba.git@a64cb2ed13e595a16ad056a128ef24594b811504", "emcee", # numba is the load-bearing pin of this whole environment: its numpy ceiling # (numba 0.66 -> numpy<2.5) is what keeps the resolver from drifting numpy diff --git a/src/sp_validation/b_modes.py b/src/sp_validation/b_modes.py index fdc2ddf8..083dafd1 100644 --- a/src/sp_validation/b_modes.py +++ b/src/sp_validation/b_modes.py @@ -1,8 +1,9 @@ """ B-mode analysis functions for weak lensing validation. -This module contains pure E/B mode decomposition, COSEBIs analysis, -and semi-analytical covariance calculations extracted from CosmologyValidation. +Pure E/B modes (Schneider et al. 2022) as a fixed linear operator on the fine +ξ± grid with their exact covariance, COSEBIs, and the χ²/PTE and plotting +helpers shared by both. """ import warnings @@ -11,19 +12,12 @@ import numpy as np import seaborn as sns import tqdm -import treecorr -from cs_util.cosmo import get_theo_xi from mpl_toolkits.axes_grid1 import make_axes_locatable -from scipy import sparse, stats +from scipy import stats _EB_KEYS = ("xip_E", "xim_E", "xip_B", "xim_B", "xip_amb", "xim_amb") -def _eb_vector(modes): - """Pure-E/B modes concatenated in ``_EB_KEYS`` order, the covariance layout.""" - return np.concatenate([modes[k] for k in _EB_KEYS]) - - def find_conservative_scale_cut_key(results, requested_scale_cut): """ Find scale cut key that conservatively fits within requested range. @@ -97,6 +91,23 @@ def bins_from_edges(left_edges, right_edges, min_scale=None, max_scale=None): return start_bin, stop_bin +def bins_from_scale_cut(left_edges, right_edges, scale_cut): + """Reporting bins inside a pure-E/B scale cut, as ``(start_bin, stop_bin)``. + + Each end of ``scale_cut`` snaps to the nearest reporting edge in log θ — + the rule that snaps reporting edges onto the fine grid — so a cut placed + on a nominal edge selects the same bins whichever fine grid the edges + were snapped to. ``stop_bin`` is exclusive. + """ + edges = np.append(left_edges, right_edges[-1]) + start_bin, stop_bin = ( + int(np.abs(np.log(edges) - np.log(cut)).argmin()) for cut in scale_cut + ) + if stop_bin <= start_bin: + raise RuntimeError(f"scale cut {scale_cut} selects no bins") + return start_bin, stop_bin + + def log_bin_edges(min_sep, max_sep, nbins): """TreeCorr ``Log`` bin edges — the grid a binning defines. @@ -125,215 +136,253 @@ def correlation_from_covariance(covariance): return covariance / np.outer(stdev, stdev) -def calculate_pure_eb_correlation( - gg, - gg_int, - var_method="jackknife", - cov_path_int=None, - cosmo_cov=None, - n_samples=1000, - z_dist=None, -): +def hartlap_factor(npatch, dof): + """Hartlap (2007) debiasing of an inverse covariance. + + ``(N - p - 2) / (N - 1)`` for a jackknife covariance from ``N = npatch`` + patches inverted over ``p = dof`` data points; exactly 1 for an analytic + covariance, which ``npatch=None`` denotes. """ - Calculate pure E/B modes from correlation function objects. + return 1.0 if npatch is None else (npatch - dof - 2) / (npatch - 1) - Parameters - ---------- - gg : treecorr.GGCorrelation - Correlation function for reporting binning (coarser binning for final results) - gg_int : treecorr.GGCorrelation - Correlation function for integration binning (fine binning for numerical - integration) - var_method : str, optional - Variance method ("jackknife" or "bootstrap") - cov_path_int : str, optional - Path to integration covariance matrix for semi-analytical calculation - cosmo_cov : pyccl.Cosmology, optional - Cosmology for theoretical predictions in semi-analytical covariance - n_samples : int, optional - Number of Monte Carlo samples for semi-analytical covariance - z_dist : 2D array, optional - Redshift distribution; - z_dist[:, 0] = z, z_dist[:, 1] = n(z) - Returns - ------- - dict - Dictionary containing pure E/B mode results and covariance - """ - # Calculate min_sep and max_sep from gg object - min_sep, max_sep = gg.left_edges[0], gg.right_edges[-1] - - def pure_EB(corrs): - gg, gg_int = corrs - return pure_eb_from_xi( - theta_report=gg.meanr, - xip_report=gg.xip, - xim_report=gg.xim, - theta_int=gg_int.meanr, - xip_int=gg_int.xip, - xim_int=gg_int.xim, - tmin=min_sep, - tmax=max_sep, - ) +def _reporting_binning(weight_int, edges_int, reporting_edges): + """Reporting bins as unions of fine bins, and their pair-weighted average. - # The results dict is self-describing: the grids it was measured on travel - # with the modes, so every consumer downstream works from values alone. - results = { - "theta": gg.meanr, - "left_edges": gg.left_edges, - "right_edges": gg.right_edges, - "xip": gg.xip, - "xim": gg.xim, - "var_xip": gg.varxip, - "var_xim": gg.varxim, - "theta_int": gg_int.meanr, - "xip_int": gg_int.xip, - "xim_int": gg_int.xim, - "n_eff": n_samples if cov_path_int is not None else gg.npatch1, - } - results.update(pure_EB([gg, gg_int])) + Each requested reporting edge snaps to the nearest fine edge (in log θ), so + every reporting bin is a whole number of fine bins. Row ``i`` of the + ``(n_report, n_fine)`` matrix weights the fine bins of reporting bin ``i`` + by their TreeCorr pair weight ``Σ w_i w_j`` and sums to one. TreeCorr's ξ± + and ``meanr`` are averages over pairs under that same weight, so each row + reproduces what TreeCorr would measure on the snapped bin. - if cov_path_int is not None: - if z_dist is None or cosmo_cov is None: - raise ValueError( - "semi-analytical covariance needs both z_dist and cosmo_cov" - ) - cov, eb_samples = pure_eb_covariance_mc( - theta=gg.meanr, - left_edges=gg.left_edges, - right_edges=gg.right_edges, - theta_int=gg_int.meanr, - cov_int=np.loadtxt(cov_path_int), - z=z_dist[:, 0], - nz=z_dist[:, 1], - cosmo=cosmo_cov, - n_samples=n_samples, + Returns ``(binning, edges)``, ``edges`` being the snapped reporting edges. + """ + weight_int = np.asarray(weight_int, dtype=float) + edges_int = np.asarray(edges_int, dtype=float) + reporting_edges = np.asarray(reporting_edges, dtype=float) + if edges_int.shape != (weight_int.size + 1,): + raise ValueError("edges_int must hold n_fine + 1 fine-grid edges") + if reporting_edges.min() < edges_int[0] or reporting_edges.max() > edges_int[-1]: + raise ValueError( + f"reporting edges [{reporting_edges.min()}, {reporting_edges.max()}] " + f"reach outside the fine grid [{edges_int[0]}, {edges_int[-1]}]" ) - results.update({"cov": cov, "eb_samples": eb_samples}) - else: - # Use existing treecorr covariance estimation - results["cov"] = treecorr.estimate_multi_cov( - [gg, gg_int], - var_method, - func=lambda x: _eb_vector(pure_EB(x)), - cross_patch_weight="match" if var_method == "jackknife" else None, + snap = np.abs(np.log(reporting_edges)[:, None] - np.log(edges_int)).argmin(axis=1) + if np.any(np.diff(snap) <= 0): + raise ValueError( + "reporting edges snap onto the same fine edge; the fine grid is too " + f"coarse for them: {reporting_edges.tolist()} -> " + f"{edges_int[snap].tolist()}" ) - - # Validate covariance matrix - try: - np.linalg.cholesky(results["cov"]) - except np.linalg.LinAlgError: - warnings.warn( - "E/B mode covariance matrix is not positive definite. " - "Chi-squared statistics may be unreliable.", - UserWarning, + binning = np.zeros((snap.size - 1, weight_int.size)) + for i, (lo, hi) in enumerate(zip(snap[:-1], snap[1:])): + binning[i, lo:hi] = weight_int[lo:hi] + weight = binning.sum(axis=1) + if np.any(weight == 0): + empty = np.flatnonzero(weight == 0).tolist() + raise ValueError(f"reporting bins {empty} hold no integration-grid pairs") + return binning / weight[:, None], edges_int[snap] + + +def _fixed_quadrature_operator(edges_int): + """The Schneider (2022) transform on the fine grid, ``_EB_KEYS`` order. + + The single call into cosmo_numba. Its interpolator places the samples on + a regular grid in log θ, so the transform runs on the log-uniform nodes of + the fine bins — their geometric centres, TreeCorr's ``rnom`` — and is + evaluated at every node, with ``[tmin, tmax]`` at their extent. That + makes it a function of the fine grid alone (the zero-padded ξ− window + extends from the evaluation grid). Returns the ``(6, n_fine, 2 * n_fine)`` + stack of the six matrices acting on ``[xi_+; xi_-]``; rows whose + integration support is too small near the grid edges are NaN. + """ + from cosmo_numba.B_modes.schneider2022 import get_pure_EB_operator + + log_edges = np.log(np.asarray(edges_int, dtype=float)) + step = np.diff(log_edges) + if np.ptp(step) > 1e-10 * np.mean(step): + raise ValueError( + "the pure-E/B transform needs log-uniform fine-grid edges (TreeCorr " + f"Log binning); their log spacing varies by {np.ptp(step):.3e}" ) + nodes = np.exp(0.5 * (log_edges[:-1] + log_edges[1:])) + return np.stack( + get_pure_EB_operator( + nodes, + nodes, + tmin=nodes[0] * (1 - 1e-9), + tmax=nodes[-1] * (1 + 1e-9), + local_from_int=True, + ) + ) - return results +def pure_eb_operator(weight_int, edges_int, reporting_edges): + """The pure-E/B estimator as one matrix on the fine ξ± grid. -def pure_eb_from_xi( - theta_report, xip_report, xim_report, theta_int, xip_int, xim_int, tmin, tmax -): - """Pure-E/B correlation functions from ξ± arrays through the pipeline kernel. + The Schneider et al. (2022) transform is evaluated with fixed-quadrature + weights at the log-uniform fine-grid nodes (TreeCorr ``rnom``), + integrating over the whole fine grid, and + the six pure modes are then averaged into the reporting bins with TreeCorr + pair weights. The reporting edges snap to the nearest fine edges, so each + reporting bin is a union of fine bins. Both steps are linear and + independent of the ξ± values, so the estimator is + ``K = (I_6 ⊗ P) · M`` and + + [xip_E; xim_E; xip_B; xim_B; xip_amb; xim_amb] = K @ [xip_int; xim_int] - The one place this module calls cosmo_numba's Schneider (2022) transform. + with ``K`` of shape ``(6 * n_report, 2 * n_fine)``. - ``tmin``/``tmax`` are the reporting correlation's TreeCorr *bin edges* - (``gg.left_edges[0]`` / ``gg.right_edges[-1]``). The reporting grid must be a - strict sub-range of the integration grid: a reporting point on the - integration boundary has no interior support and comes back NaN. + Parameters + ---------- + weight_int : array_like + TreeCorr pair weight ``Σ w_i w_j`` per fine bin (``gg.weight``), the + averaging weights. + edges_int : array_like + The ``n_fine + 1`` log-uniform fine-grid bin edges. The transform's + ``[tmin, tmax]`` is the extent of their centres, so they must reach + beyond the reporting range on both sides. + reporting_edges : array_like + Requested reporting-bin edges, ``n_report + 1`` of them. Returns ------- - dict - Keyed by ``_EB_KEYS`` (xip_E, xim_E, xip_B, xim_B, xip_amb, xim_amb). + operator : numpy.ndarray + ``K``, shape ``(6 * n_report, 2 * n_fine)``. + binning : numpy.ndarray + ``P``, the ``(n_report, n_fine)`` pair-weighted average. + edges : numpy.ndarray + The reporting edges actually used, each a fine-grid edge. """ - from cosmo_numba.B_modes.schneider2022 import get_pure_EB_modes - - modes = get_pure_EB_modes( - theta=np.asarray(theta_report), - xip=np.asarray(xip_report), - xim=np.asarray(xim_report), - theta_int=np.asarray(theta_int), - xip_int=np.asarray(xip_int), - xim_int=np.asarray(xim_int), - tmin=tmin, - tmax=tmax, - parallel=True, - ) - return dict(zip(_EB_KEYS, (np.asarray(m) for m in modes))) + binning, edges = _reporting_binning(weight_int, edges_int, reporting_edges) + nodes = np.flatnonzero(binning.any(axis=0)) + # Only the rows P averages enter; NaN edge rows outside them are dropped, + # since 0 · NaN would still poison the product. + transform = _fixed_quadrature_operator(edges_int)[:, nodes] + undetermined = { + key: int(np.count_nonzero(~np.isfinite(rows).all(axis=1))) + for key, rows in zip(_EB_KEYS, transform) + if not np.isfinite(rows).all() + } + if undetermined: + raise ValueError( + "pure-E/B operator rows are under-determined (reporting bins too " + f"close to the integration-grid edge): {undetermined}" + ) + operator = np.vstack([binning[:, nodes] @ rows for rows in transform]) + return operator, binning, edges -def pure_eb_covariance_mc( - *, - theta, - left_edges, - right_edges, +def calculate_pure_eb_correlation( theta_int, - cov_int, - z, - nz, - cosmo, - n_samples=1000, + xip_int, + xim_int, + weight_int, + edges_int, + cov_xi, + reporting_edges, + *, + npatch=None, ): - """Pure-E/B covariance by Monte Carlo through the same kernel as the modes. + """Pure E/B modes and their exact covariance from fine-grid ξ±. + + The modes are :func:`pure_eb_operator` applied to ``[xip_int; xim_int]``, + and the covariance is ``K C_xi Kᵀ`` — exact for whatever ξ± covariance is + supplied, analytic or jackknife. ``npatch`` records which: the jackknife + patch count behind ``cov_xi``, or ``None`` for an analytic covariance. It + travels in the results and sets the Hartlap factor of every χ² built on + them (:func:`hartlap_factor`). + + The reporting bins are unions of fine bins (the requested edges snap to + the nearest fine edges), and their ``theta``, ``xip``/``xim`` and + variances are the same pair-weighted average of the fine grid, so + ``xi_± = E ± B + amb`` holds bin by bin and ``theta``/``xip``/``xim`` are + what TreeCorr would measure on the snapped bins. - ξ± draws come from ``cov_int``, a ξ± covariance on the integration grid, - around the theory mean for ``(z, nz)`` under ``cosmo``; each draw is binned - down to the reporting grid and pushed through :func:`pure_eb_from_xi`. The - covariance of the transformed draws is the result, so it depends on the - covariance model and the grids, never on the measured data vector. + Parameters + ---------- + theta_int, xip_int, xim_int, weight_int : array_like + Fine-grid ``meanr``, ξ±, and TreeCorr pair weight ``Σ w_i w_j``. + ``meanr`` enters only the reported ``theta``; the transform runs on + the log-uniform nodes of ``edges_int``. + edges_int : array_like + The ``n_fine + 1`` log-uniform fine-grid bin edges. + cov_xi : array_like + ``(2 n_fine, 2 n_fine)`` covariance of ``[xip_int; xim_int]``. + reporting_edges : array_like + Requested reporting-bin edges, ``n_report + 1`` of them. + npatch : int, optional + Jackknife patch count behind ``cov_xi``; ``None`` if it is analytic. - Returns ``(cov, eb_samples)`` — the covariance in ``_EB_KEYS`` order and - the draws behind it. + Returns + ------- + dict + The six ``_EB_KEYS`` mode arrays, ``cov`` (in ``_EB_KEYS`` block + order), ``npatch``, the reporting grid (``theta``, the snapped + ``left_edges``/``right_edges``, ``xip``, ``xim``, ``var_xip``, + ``var_xim``) and the fine-grid inputs (``theta_int``, ``xip_int``, + ``xim_int``, ``weight_int``, ``edges_int``). """ - theta, theta_int = np.asarray(theta), np.asarray(theta_int) - nbins_int = len(theta_int) - - # Each reporting bin averages the integration bins that fall inside it. - reporting_bin_edges = np.concatenate([left_edges, [right_edges[-1]]]) - bin_indices = np.digitize(theta_int, reporting_bin_edges) - 1 - valid_mask = (bin_indices >= 0) & (bin_indices < len(theta)) - row_indices, col_indices = (bin_indices[valid_mask], np.where(valid_mask)[0]) - binning_matrix = sparse.csr_matrix( - (np.ones(len(row_indices)), (row_indices, col_indices)), - shape=(len(theta), nbins_int), + if npatch is not None and npatch < 2: + raise ValueError(f"a jackknife covariance needs npatch > 1, not {npatch}") + theta_int, xip_int, xim_int = ( + np.asarray(a, dtype=float) for a in (theta_int, xip_int, xim_int) ) - row_sums = np.array(binning_matrix.sum(axis=1)).flatten() - binning_matrix = sparse.diags(1 / row_sums) @ binning_matrix - - # One n(z) gives one tracer pair: get_theo_xi's single (xi+, xi-) entry. - (xi_pm,) = get_theo_xi( - theta=theta_int, z=z, nz=nz, backend="ccl", cosmo=cosmo - ).values() - mean_int = np.concatenate(xi_pm) - samples_int = np.random.multivariate_normal(mean_int, cov_int, size=n_samples) - samples_int_xip, samples_int_xim = ( - samples_int[:, :nbins_int], - samples_int[:, nbins_int:], - ) - samples_rep_xip = (binning_matrix @ samples_int_xip.T).T - samples_rep_xim = (binning_matrix @ samples_int_xim.T).T - - def eb_draw(i): - modes = pure_eb_from_xi( - theta_report=theta, - xip_report=samples_rep_xip[i], - xim_report=samples_rep_xim[i], - theta_int=theta_int, - xip_int=samples_int_xip[i], - xim_int=samples_int_xim[i], - tmin=left_edges[0], - tmax=right_edges[-1], + cov_xi = np.asarray(cov_xi, dtype=float) + if np.shape(weight_int) != theta_int.shape: + raise ValueError("weight_int must have one entry per integration bin") + operator, binning, edges = pure_eb_operator(weight_int, edges_int, reporting_edges) + if cov_xi.shape != (operator.shape[1],) * 2: + raise ValueError( + f"cov_xi has shape {cov_xi.shape}; the fine grid needs " + f"{(operator.shape[1],) * 2}" ) - return _eb_vector(modes) - eb_samples = np.array( - [eb_draw(i) for i in tqdm.tqdm(range(n_samples), desc="MC samples")] + n_report = len(edges) - 1 + modes = operator @ np.concatenate([xip_int, xim_int]) + n_fine = theta_int.size + var_xip, var_xim = ( + np.einsum("ij,jk,ik->i", binning, block, binning) + for block in (cov_xi[:n_fine, :n_fine], cov_xi[n_fine:, n_fine:]) ) - return np.cov(eb_samples.T), eb_samples + results = { + "theta": binning @ theta_int, + "left_edges": edges[:-1], + "right_edges": edges[1:], + "xip": binning @ xip_int, + "xim": binning @ xim_int, + "var_xip": var_xip, + "var_xim": var_xim, + "theta_int": theta_int, + "xip_int": xip_int, + "xim_int": xim_int, + "weight_int": np.asarray(weight_int, dtype=float), + "edges_int": np.asarray(edges_int, dtype=float), + "cov": operator @ cov_xi @ operator.T, + "npatch": npatch, + } + for i, key in enumerate(_EB_KEYS): + results[key] = modes[i * n_report : (i + 1) * n_report] + + # The B-mode block is what every χ² inverts; the ambiguous blocks are + # ill-conditioned by construction. + b_block = slice(2 * n_report, 4 * n_report) + try: + np.linalg.cholesky(results["cov"][b_block, b_block]) + except np.linalg.LinAlgError: + warnings.warn( + "B-mode covariance is not positive definite. " + "Chi-squared statistics may be unreliable.", + UserWarning, + ) + + return results + + +def covariance_label(npatch): + """How a pure-E/B covariance was made, from its ``npatch`` record.""" + return "analytic" if npatch is None else f"jackknife ({npatch} patches)" def calculate_cosebis(gg, nmodes=10, scale_cuts=None, cov_path=None): @@ -477,8 +526,9 @@ def calculate_eb_statistics(results): ---------- results : dict Pure E/B results: the six mode arrays, the ``cov`` block, the reporting - ``theta``, and ``n_eff`` — the realisation count behind the covariance - (jackknife patches or MC draws), which sets the Hartlap debiasing + ``theta``, and ``npatch`` — the jackknife patch count behind the + covariance, or ``None`` for an analytic one — which sets the Hartlap + factor (:func:`hartlap_factor`) Returns ------- @@ -486,7 +536,7 @@ def calculate_eb_statistics(results): Updated results dictionary with PTE matrices and statistics """ nbins = len(results["theta"]) - n_eff = results["n_eff"] + npatch = results["npatch"] # Extract covariance blocks and standard deviations cov = results["cov"] @@ -508,7 +558,7 @@ def calculate_eb_statistics(results): for start_bin, stop_bin in combinations: nbins_eff = stop_bin - start_bin - hartlap_factor = (n_eff - nbins_eff - 2) / (n_eff - 1) + hartlap = hartlap_factor(npatch, nbins_eff) # Individual B-mode chi-squared calculations data_slices = [results[f"{xi}_B"][start_bin:stop_bin] for xi in ["xip", "xim"]] @@ -517,7 +567,7 @@ def calculate_eb_statistics(results): for xi in ["xip", "xim"] ] chi2_values = [ - hartlap_factor * (data @ np.linalg.solve(cov, data)) + hartlap * (data @ np.linalg.solve(cov, data)) for data, cov in zip(data_slices, cov_slices) ] @@ -536,7 +586,7 @@ def calculate_eb_statistics(results): cov_combined = np.block( [[cov_xip_block, cov_cross_block], [cov_cross_block.T, cov_xim_block]] ) - chi2_combined = hartlap_factor * ( + chi2_combined = hartlap_factor(npatch, 2 * nbins_eff) * ( data_combined @ np.linalg.solve(cov_combined, data_combined) ) pte_combined[start_bin, stop_bin - 1] = stats.chi2.sf( @@ -557,9 +607,10 @@ def plot_integration_vs_reporting(results, output_path, version): Parameters ---------- results : dict - Pure E/B results carrying both grids (``theta``/``xip``/``xim`` and the - ``theta_int``/``xip_int``/``xim_int`` counterparts), plus the reporting - ``var_xip``/``var_xim`` the error bars use + Pure E/B results carrying the fine grid (``theta_int``/``xip_int``/ + ``xim_int``) and its pair-weighted average into the reporting bins + (``theta``/``xip``/``xim``, with the ``var_xip``/``var_xim`` the error + bars use) output_path : str Output file path for the plot version : str @@ -614,7 +665,7 @@ def plot_integration_vs_reporting(results, output_path, version): def _get_pte_from_scale_cut(pte_matrix, edges, scale_cut): """ - Extract PTE value from matrix based on scale cut range using conservative logic. + Extract PTE value from matrix at a scale cut (:func:`bins_from_scale_cut`). Parameters ---------- @@ -637,13 +688,7 @@ def _get_pte_from_scale_cut(pte_matrix, edges, scale_cut): # Return full-range PTE (first row, last column) return pte_matrix[0, nbins - 1] - min_scale, max_scale = scale_cut - - start_bin, stop_bin = bins_from_edges(left_edges, right_edges, min_scale, max_scale) - - # Ensure valid range, otherwise fallback to full range - if stop_bin <= start_bin or start_bin >= nbins or stop_bin <= 0: - raise RuntimeError("Invalid scale cut range") + start_bin, stop_bin = bins_from_scale_cut(left_edges, right_edges, scale_cut) return pte_matrix[start_bin, stop_bin - 1] @@ -678,12 +723,16 @@ def plot_pure_eb_correlations( # Calculate combined PTE using off-diagonal covariance blocks # Get scale cuts for both xi+ and xi- if fiducial_xip_scale_cut is not None: - xip_start_bin, xip_stop_bin = bins_from_edges(*edges, *fiducial_xip_scale_cut) + xip_start_bin, xip_stop_bin = bins_from_scale_cut( + *edges, fiducial_xip_scale_cut + ) else: xip_start_bin, xip_stop_bin = 0, nbins if fiducial_xim_scale_cut is not None: - xim_start_bin, xim_stop_bin = bins_from_edges(*edges, *fiducial_xim_scale_cut) + xim_start_bin, xim_stop_bin = bins_from_scale_cut( + *edges, fiducial_xim_scale_cut + ) else: xim_start_bin, xim_stop_bin = 0, nbins @@ -713,15 +762,7 @@ def plot_pure_eb_correlations( # Calculate combined chi-squared with Hartlap factor nbins_eff = len(xip_B_data) + len(xim_B_data) - - # Determine effective number of samples for Hartlap correction - if "eb_samples" in results: # Semi-analytical case - n_eff = results["eb_samples"].shape[0] - else: # Jackknife case - n_eff = results["n_eff"] - - hartlap_factor = (n_eff - nbins_eff - 2) / (n_eff - 1) - chi2_combined = hartlap_factor * ( + chi2_combined = hartlap_factor(results["npatch"], nbins_eff) * ( data_combined.T @ np.linalg.solve(cov_combined, data_combined) ) combined_pte = stats.chi2.sf(chi2_combined, nbins_eff) @@ -818,11 +859,10 @@ def plot_pure_eb_correlations( scale_cuts = [(fiducial_xip_scale_cut, 0), (fiducial_xim_scale_cut, 1)] for scale_cut, ax_idx in scale_cuts: if scale_cut is not None: - min_scale, max_scale = scale_cut xlim = original_xlims[ax_idx] - # Use conservative scale_cut_to_bins helper for consistency - start_bin, stop_bin = bins_from_edges(*edges, min_scale, max_scale) + # The same bins the PTEs use + start_bin, stop_bin = bins_from_scale_cut(*edges, scale_cut) # Show excluded regions based on bin edges used in PTE calculation # Lower exclusion: bins 0 to start_bin-1 are excluded @@ -1105,8 +1145,7 @@ def plot_pte_2d_heatmaps( fiducial_scale_cuts = [fiducial_xip_scale_cut, fiducial_xim_scale_cut] for ax_idx, fiducial_scale_cut in enumerate(fiducial_scale_cuts): if fiducial_scale_cut is not None: - min_scale, max_scale = fiducial_scale_cut - start_bin, stop_bin = bins_from_edges(*edges, min_scale, max_scale) + start_bin, stop_bin = bins_from_scale_cut(*edges, fiducial_scale_cut) if stop_bin > start_bin and start_bin < nbins and stop_bin > 0: rect_x = start_bin @@ -1157,7 +1196,7 @@ def plot_pte_2d_heatmaps( plt.savefig(output_path, dpi=300, bbox_inches="tight") -def plot_eb_covariance_matrix(cov_matrix, var_method, output_path, version): +def plot_eb_covariance_matrix(cov_matrix, label, output_path, version): """ Plot E/B mode covariance matrix as correlation matrix. @@ -1165,8 +1204,8 @@ def plot_eb_covariance_matrix(cov_matrix, var_method, output_path, version): ---------- cov_matrix : numpy.ndarray Covariance matrix from E/B mode analysis - var_method : str - Variance method used for the analysis + label : str + How the covariance was made (:func:`covariance_label`) output_path : str Output file path for the plot version : str @@ -1198,7 +1237,7 @@ def plot_eb_covariance_matrix(cov_matrix, var_method, output_path, version): divider = make_axes_locatable(ax) cax = divider.append_axes("right", size="5%", pad=0.1) plt.colorbar(im, cax=cax) - ax.set_title(f"{version}: {var_method} correlation matrix") + ax.set_title(f"{version}: {label} correlation matrix") plt.savefig(output_path, dpi=300, bbox_inches="tight") @@ -1272,7 +1311,9 @@ def save_pure_eb_results(results, output_path): Output .npz file path """ # Data vectors and covariance - save_dict = {"theta": results["theta"], "cov": results["cov"]} + save_dict = { + key: results[key] for key in ("theta", "left_edges", "right_edges", "cov") + } for key in _EB_KEYS: save_dict[key] = results[key] @@ -1280,13 +1321,9 @@ def save_pure_eb_results(results, output_path): for key, matrix in results.get("pte_matrices", {}).items(): save_dict[f"pte_matrices_{key}"] = matrix - # Metadata - save_dict["n_eff"] = np.array(results["n_eff"]) - if "eb_samples" in results: - save_dict["var_method"] = np.array("semi-analytic") - save_dict["n_samples"] = np.array(results["eb_samples"].shape[0]) - else: - save_dict["var_method"] = np.array("jackknife") + # The jackknife patch count, stored only for a jackknife covariance. + if results["npatch"] is not None: + save_dict["npatch"] = np.array(results["npatch"]) np.savez(output_path, **save_dict) print(f"Saved pure E/B results to {output_path}") diff --git a/src/sp_validation/cosmo_val/core.py b/src/sp_validation/cosmo_val/core.py index e4957afe..422c3a9c 100644 --- a/src/sp_validation/cosmo_val/core.py +++ b/src/sp_validation/cosmo_val/core.py @@ -11,6 +11,7 @@ from ..b_modes import ( _get_pte_from_scale_cut, + covariance_label, find_conservative_scale_cut_key, ) from ..statistics import chi2_and_pte @@ -608,11 +609,7 @@ def summarize_bmodes(self, fiducial_scale_cut=(12, 83), versions=None): ) except (KeyError, RuntimeError): pass - cov_methods.add( - "semi-analytic" - if "eb_samples" in res - else f"jackknife ({res['n_eff']} patches)" - ) + cov_methods.add(covariance_label(res["npatch"])) # COSEBIs PTE from stored results if ver in self._cosebis_results: diff --git a/src/sp_validation/cosmo_val/pure_eb.py b/src/sp_validation/cosmo_val/pure_eb.py index 90157456..81bd22f0 100644 --- a/src/sp_validation/cosmo_val/pure_eb.py +++ b/src/sp_validation/cosmo_val/pure_eb.py @@ -10,6 +10,7 @@ from ..b_modes import ( calculate_eb_statistics, calculate_pure_eb_correlation, + covariance_label, plot_eb_covariance_matrix, plot_integration_vs_reporting, plot_pte_2d_heatmaps, @@ -28,114 +29,68 @@ def calculate_pure_eb( min_sep_int=0.08, max_sep_int=300, nbins_int=1000, - npatch=256, - var_method="jackknife", + npatch=None, cov_path_int=None, - cosmo_cov=None, - n_samples=1000, ): """ Calculate the pure E/B modes for the given catalog version. - The class instance's treecorr_config will be used for the "reporting" binning - by default, but any kwargs passed to this function will overwrite the defaults. + + ξ± is measured on the fine integration grid only; the reporting + binning (the instance's treecorr_config unless overridden) enters as + bin edges, snapped onto the fine edges, into which + :func:`~sp_validation.b_modes.pure_eb_operator` averages the modes. Parameters ---------- version : str The catalog version to compute the pure E/B modes for. - min_sep : float, optional - Minimum separation for the reporting binning. Defaults to the value in - self.treecorr_config if not provided. - max_sep : float, optional - Maximum separation for the reporting binning. Defaults to the value in - self.treecorr_config if not provided. - nbins : int, optional - Number of bins for the reporting binning. Defaults to the value in - self.treecorr_config if not provided. - min_sep_int : float, optional - Minimum separation for the integration binning. Defaults to 0.08. - max_sep_int : float, optional - Maximum separation for the integration binning. Defaults to 300. - nbins_int : int, optional - Number of bins for the integration binning. Defaults to 1000. + 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 + extend beyond the reporting range on both sides. npatch : int, optional - Number of patches for the jackknife or bootstrap resampling. Defaults to - the value in self.npatch if not provided. - var_method : str, optional - Variance estimation method. Defaults to "jackknife". + Jackknife patch count. Defaults to self.npatch. cov_path_int : str, optional - Path to the covariance matrix for the reporting binning. Replaces the - treecorr covariance matrix if provided, meaning that var_method has no - effect on the results although it is still passed to - CosmologyValidation.calculate_2pcf. - cosmo_cov : pyccl.Cosmology, optional - Cosmology object to use for theoretical xi+/xi- predictions in the - semi-analytical covariance calculation. Defaults to self.cosmo if not - provided. - n_samples : int, optional - Number of Monte Carlo samples for semi-analytical covariance propagation. - Defaults to 1000. + Analytic ξ± covariance on the integration grid. Without it the + covariance is the jackknife of the integration-grid ξ±. Returns ------- dict - A dictionary containing the following keys: - - - "xip_E": Pure E-mode correlation function for xi+. - - "xim_E": Pure E-mode correlation function for xi-. - - "xip_B": Pure B-mode correlation function for xi+. - - "xim_B": Pure B-mode correlation function for xi-. - - "xip_amb": Ambiguity mode for xi+. - - "xim_amb": Ambiguity mode for xi-. - - "cov": Covariance matrix for the pure E/B modes. - - "theta", "left_edges", "right_edges": Reporting-grid bin centres - and edges. - - "xip", "xim", "var_xip", "var_xim": Reporting-grid xi+/xi- and - their variances. - - "theta_int", "xip_int", "xim_int": Integration-grid xi+/xi-. - - "n_eff": Realisation count behind "cov" (jackknife patches or - MC draws), which sets the Hartlap debiasing. - - "eb_samples": (only when using semi-analytical covariance) Semi-analytic - EB samples used for covariance calculation. Shape: (n_samples, 6*nbins) - - Notes - ----- - - A shared patch file is used for the reporting and integration binning, - and is created if it does not exist. + The results of :func:`~sp_validation.b_modes.calculate_pure_eb_correlation`: + the six pure-mode arrays, their covariance ``cov`` and its + ``npatch`` record (``None`` for an analytic covariance), the + reporting grid and the integration-grid ξ±. """ self.print_start(f"Computing {version} pure E/B") - # Set up parameters with defaults - npatch = npatch or self.npatch - - # Create TreeCorr configurations - treecorr_config = self._binning(min_sep, max_sep, nbins) - treecorr_config_int = self._binning(min_sep_int, max_sep_int, nbins_int) - - # Calculate correlation functions - gg = self.calculate_2pcf(version, npatch=npatch, **treecorr_config) - gg_int = self.calculate_2pcf(version, npatch=npatch, **treecorr_config_int) - - # Get redshift distribution if using analytic covariance - z_dist = ( - np.column_stack(self.get_redshift(version)) - if cov_path_int is not None - else None + reporting = self._binning(min_sep, max_sep, nbins) + gg_int = self.calculate_2pcf( + version, + npatch=npatch, + **self._binning(min_sep_int, max_sep_int, nbins_int), ) - # Delegate to b_modes module - results = calculate_pure_eb_correlation( - gg=gg, - gg_int=gg_int, - var_method=var_method, - cov_path_int=cov_path_int, - cosmo_cov=cosmo_cov, - n_samples=n_samples, - z_dist=z_dist, + if cov_path_int is not None: + cov_xi, npatch = np.loadtxt(cov_path_int), None + else: + cov_xi = gg_int.estimate_cov("jackknife", cross_patch_weight="match") + npatch = gg_int.npatch1 + + return calculate_pure_eb_correlation( + gg_int.meanr, + gg_int.xip, + gg_int.xim, + gg_int.weight, + np.append(gg_int.left_edges, gg_int.right_edges[-1]), + cov_xi, + np.geomspace( + reporting["min_sep"], reporting["max_sep"], reporting["nbins"] + 1 + ), + npatch=npatch, ) - return results - def plot_pure_eb( self, versions=None, @@ -149,12 +104,8 @@ def plot_pure_eb( max_sep_int=300, nbins_int=1000, npatch=None, - var_method="jackknife", cov_path_int=None, - cosmo_cov=None, - n_samples=1000, results=None, - **kwargs, ): """ Generate comprehensive pure E/B mode analysis plots. @@ -182,29 +133,21 @@ def plot_pure_eb( (default: 0.08-300 arcmin, 1000 bins) npatch : int, optional Number of patches for jackknife covariance. Uses self.npatch if None. - var_method : str - Variance method ("jackknife" or "semi-analytic"). - Automatically set to "semi-analytic" when cov_path_int is provided. cov_path_int : str, optional - Path to integration covariance matrix for semi-analytical calculation - cosmo_cov : pyccl.Cosmology, optional - Cosmology for theoretical predictions in semi-analytical covariance - n_samples : int - Number of Monte Carlo samples for semi-analytical covariance (default: 1000) + Analytic ξ± covariance on the integration grid; the jackknife is + used without it. results : dict or list, optional Precalculated results to avoid recomputation. Can be a single results dict for one version, or a list of results dicts for multiple versions. If None (default), results will be calculated using calculate_pure_eb. - **kwargs : dict - Additional arguments passed to calculate_eb_statistics Notes ----- This function orchestrates the full E/B mode analysis workflow: - Uses instance configuration as defaults for unspecified parameters - - Automatically switches to analytical variance when theoretical - covariance provided + - Uses the analytic covariance when cov_path_int is given, the + jackknife otherwise - Generates standardized output file naming based on all analysis parameters - Delegates individual plot generation to specialized functions in @@ -215,9 +158,7 @@ def plot_pure_eb( output_dir = output_dir or self.cc["paths"]["output"] npatch = npatch or self.npatch - # Override var_method to analytic when cov_path_int is provided - if cov_path_int is not None: - var_method = "semi-analytic" + var_method = "jackknife" if cov_path_int is None else "analytic" # Use treecorr_config defaults for reporting scale binning min_sep = min_sep or self.treecorr_config["min_sep"] @@ -266,14 +207,11 @@ def plot_pure_eb( 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, ) # Calculate E/B statistics for all bin combinations - version_results = calculate_eb_statistics(version_results, **kwargs) + version_results = calculate_eb_statistics(version_results) # Integration vs Reporting comparison plot plot_integration_vs_reporting( @@ -303,7 +241,7 @@ def plot_pure_eb( # Covariance matrix plot plot_eb_covariance_matrix( version_results["cov"], - var_method, + covariance_label(version_results["npatch"]), out_stub + "_covariance.png", version, ) diff --git a/src/sp_validation/sacc_io.py b/src/sp_validation/sacc_io.py index 99a0efe2..ff11db0b 100644 --- a/src/sp_validation/sacc_io.py +++ b/src/sp_validation/sacc_io.py @@ -115,8 +115,8 @@ } # PURE_TYPES key order is the insertion order of the six pure-EB blocks — # matches b_modes._EB_KEYS, whose order is the [xip_E; xim_E; xip_B; xim_B; -# xip_amb; xim_amb] layout of the treecorr/MC pure-EB covariance -# (b_modes.calculate_eb_statistics, ~L392). +# xip_amb; xim_amb] block layout of the pure-EB covariance +# (b_modes.calculate_pure_eb_correlation). PURE_KEYS = tuple(PURE_TYPES) RHO_PLUS = "psf_rho{k}_xi_plus" @@ -580,6 +580,11 @@ def get_xi(s, bins, *, grid): return _get_pm(s, XI_PLUS, XI_MINUS, _pair(bins), grid=grid) +def get_xi_weight(s, bins, *, grid): + """Return the TreeCorr pair weights stored with :func:`add_xi`'s ξ+ points.""" + return _tag(s, XI_PLUS, _pair(bins), "weight", grid=grid) + + def get_pseudo_cl(s, bins): """Return ``(ell_eff, cl_ee, cl_bb, cl_eb, window)`` for one tracer pair. diff --git a/src/sp_validation/tests/conftest.py b/src/sp_validation/tests/conftest.py index 40424120..004abc83 100644 --- a/src/sp_validation/tests/conftest.py +++ b/src/sp_validation/tests/conftest.py @@ -10,10 +10,11 @@ @pytest.fixture def pure_eb_xi(): - """Committed ξ± of the synthetic coherent-shear catalogue. + """Committed fine-grid ξ± of the synthetic coherent-shear catalogue. - Exact-binning reporting [15, 70]′ in 6 bins and integration [1, 300]′ in - 600 bins, keyed by ``b_modes.pure_eb_from_xi``'s parameters. + Exact-binning integration grid [1, 300]′ in 600 bins with its pair weights + and edges, and the requested edges of a [15, 70]′ reporting grid in 6 bins, + keyed by ``b_modes.calculate_pure_eb_correlation``'s parameters. """ with np.load(PURE_EB_XI) as npz: return {k: (v.item() if v.ndim == 0 else v) for k, v in npz.items()} diff --git a/src/sp_validation/tests/data/pure_eb_xi_fixture.npz b/src/sp_validation/tests/data/pure_eb_xi_fixture.npz index fa67b006..193715c4 100644 Binary files a/src/sp_validation/tests/data/pure_eb_xi_fixture.npz and b/src/sp_validation/tests/data/pure_eb_xi_fixture.npz differ diff --git a/src/sp_validation/tests/test_b_modes.py b/src/sp_validation/tests/test_b_modes.py index 20a9aa67..2ed8a341 100644 --- a/src/sp_validation/tests/test_b_modes.py +++ b/src/sp_validation/tests/test_b_modes.py @@ -2,8 +2,8 @@ This module pins the numeric behavior of the pure E/B-mode helpers in ``sp_validation.b_modes`` against fixed, deterministic inputs (seeded RNG, -hand-built arrays and one committed ξ± fixture — no cluster data, no catalogue -files). +hand-built arrays and one committed fine-grid ξ± fixture — no cluster data, no +catalogue files). Every pinned literal was produced by an actual run of the estimator inside the container; a future refactor that changes the numbers must fail. @@ -49,9 +49,8 @@ def _grid_gg(): def _eb_inputs(): """Fixed seeded input for calculate_eb_statistics. - nbins=4, npatch=50 so the Hartlap factor (n_eff - nbins_eff - 2)/(n_eff-1) - is well-defined and strictly positive for every scale-cut combination. - n_eff is the jackknife patch count, as it is for a jackknife covariance. + nbins=4, npatch=50 so the Hartlap factor (npatch - p - 2)/(npatch - 1) is + well-defined and strictly positive for every scale-cut combination. The covariance is built SPD via A @ A.T + I; the B-mode vectors are O(1) so the chi-squared (and hence PTE) lands in a meaningful range rather than being saturated at 1.0. @@ -64,7 +63,7 @@ def _eb_inputs(): xim_B = rng.standard_normal(nbins) return { "theta": np.geomspace(1.0, 100.0, nbins), - "n_eff": npatch, + "npatch": npatch, "cov": cov, "xip_B": xip_B, "xim_B": xim_B, @@ -236,8 +235,9 @@ def test_calculate_eb_statistics_pte_matrices(): """Pin representative PTE-matrix entries from the full 2D E/B analysis. Inputs are fixed (seed 12345, nbins=4, npatch=50, SPD cov = A@A.T + I, - O(1) B-mode vectors). The Hartlap correction uses n_eff = 50, the patch - count behind a jackknife covariance. For each of xip_B, xim_B and combined we pin the + O(1) B-mode vectors). The Hartlap correction uses npatch = 50 over the + length of the inverted vector (2x for combined). For each of xip_B, xim_B + and combined we pin the full-range entry [0, nbins-1] (start=0, stop=nbins) and an interior entry [0, 2] (start=0, stop=3). These chi2->sf PTE values are deterministic functions of the seeded input. @@ -253,12 +253,12 @@ def test_calculate_eb_statistics_pte_matrices(): # Full-range entries [0, nbins-1]. npt.assert_allclose(pm["xip_B"][0, nbins - 1], 0.9985059590347458, rtol=1e-9) npt.assert_allclose(pm["xim_B"][0, nbins - 1], 0.9979174764123961, rtol=1e-9) - npt.assert_allclose(pm["combined"][0, nbins - 1], 0.9999930883595443, rtol=1e-9) + npt.assert_allclose(pm["combined"][0, nbins - 1], 0.9999952393605003, rtol=1e-9) # Interior entries [0, 2] (start_bin=0, stop_bin=3). npt.assert_allclose(pm["xip_B"][0, 2], 0.9991074524059739, rtol=1e-9) npt.assert_allclose(pm["xim_B"][0, 2], 0.9896253892931961, rtol=1e-9) - npt.assert_allclose(pm["combined"][0, 2], 0.9999497648089674, rtol=1e-9) + npt.assert_allclose(pm["combined"][0, 2], 0.9999590178816327, rtol=1e-9) # Structural pins: off the valid upper triangle the matrices are NaN. for key in ("xip_B", "xim_B", "combined"): @@ -268,6 +268,28 @@ def test_calculate_eb_statistics_pte_matrices(): assert np.all(np.isfinite(np.diag(m))) # single-bin cuts are valid +def test_calculate_eb_statistics_analytic_covariance_skips_hartlap(): + """An analytic covariance (npatch None) gives the plain χ² PTE. + + The full-range χ² is data·C⁻¹·data with no factor, so its PTE is pinned + directly against scipy; the jackknife PTE on the same input differs. + """ + from scipy import stats + + results, nbins = _eb_inputs() + results["npatch"] = None + pm = b_modes.calculate_eb_statistics(results)["pte_matrices"] + + cov_B = results["cov_xip_B"] + chi2 = results["xip_B"] @ np.linalg.solve(cov_B, results["xip_B"]) + npt.assert_allclose(pm["xip_B"][0, nbins - 1], stats.chi2.sf(chi2, nbins)) + npt.assert_allclose(pm["combined"][0, nbins - 1], 0.9999894806723678, rtol=1e-9) + + jackknife, _ = _eb_inputs() + pm_jk = b_modes.calculate_eb_statistics(jackknife)["pte_matrices"] + assert pm_jk["xip_B"][0, nbins - 1] != pm["xip_B"][0, nbins - 1] + + def test_calculate_eb_statistics_has_teeth(): """Teeth for #4: a 10x-louder B-mode signal must drop the full-range PTE. @@ -295,82 +317,284 @@ def test_calculate_eb_statistics_has_teeth(): # --------------------------------------------------------------------------- -# 5. pure_eb_from_xi on committed ξ± (the transform pin) +# 5. The pure-E/B operator on committed fine-grid ξ± # --------------------------------------------------------------------------- -# pure_eb_from_xi(**fixture); regenerated only when the transform is meant to move. +# calculate_pure_eb_correlation(**fixture) modes; regenerated only when the +# estimator is meant to move. _PURE_EB_PINS = { "xip_E": [ - -2.9831529669542025e-06, - -1.5008524620265777e-05, - 3.221623968725757e-07, - 1.1797672310858565e-05, - 5.715510692557323e-06, - 8.825804523824443e-07, + 0.00012671119742922553, + 0.00011647263203301798, + 0.00011277680462627065, + 0.00011334741455438306, + 9.407280835053578e-05, + 8.500501372536384e-05, ], "xim_E": [ - -4.737558091773235e-05, - -0.00010853189443993388, - -9.094825175032069e-05, - -5.826599101284694e-05, - -4.646405415748759e-05, - -1.9978028925333273e-05, + 7.684864690558884e-06, + -2.9296916553309855e-07, + 1.1150412226648596e-06, + 5.81495542664638e-06, + -5.8802091546090845e-06, + 5.245220040191294e-06, ], "xip_B": [ - 1.7069121242262332e-05, - 3.059889782373755e-05, - -4.8805399253844115e-06, - -6.999262696335271e-06, - -1.2672006989728095e-05, - -1.214149138979614e-06, + -5.872804581564311e-05, + -5.3898600777518044e-05, + -7.01611878338831e-05, + -6.447148201049941e-05, + -6.422794082425538e-05, + -5.7688865883825626e-05, ], "xim_B": [ - -0.00011478091634539627, - -5.445112002141066e-05, - -3.100806652947907e-05, - -1.0940424256759085e-05, - -5.755185146643215e-06, - -1.628217762504557e-06, + -1.2625808436389906e-05, + 8.906960918149856e-07, + 3.6497760325039334e-06, + 9.0173210532702e-06, + 5.367045395488404e-06, + 4.3985417696428825e-06, ], "xip_amb": [ - 0.00014017621792612224, - 0.0001378482153667787, - 0.0001339573019551001, - 0.00012745126271361765, - 0.00011662385105911986, - 9.851844704443032e-05, + 8.975404658793264e-05, + 8.870470395597095e-05, + 8.695112377370022e-05, + 8.401987012227805e-05, + 7.914207248993659e-05, + 7.098386083077631e-05, ], "xim_amb": [ - -4.389203999135455e-05, - 5.279277664928643e-05, - 5.800339397836051e-05, - 4.4242350610114584e-05, - 2.9902912946755567e-05, - 1.912262132836568e-05, + 1.5068975296876175e-06, + 9.035986620575128e-07, + 5.414736226368331e-07, + 3.2422141650144327e-07, + 1.9445163845662097e-07, + 1.1646044607564388e-07, ], } -def test_pure_eb_from_xi_reproduces_pins_on_committed_xi(pure_eb_xi): - """The pure-E/B transform of the committed ξ± reproduces its pins. +def _spd(n, seed): + """A seeded SPD matrix standing in for a ξ± covariance.""" + A = np.random.default_rng(seed).standard_normal((n, n)) + return A @ A.T / n + np.eye(n) + - With ξ± frozen, these pins move only when the transform does. rtol=1e-6 is - far above the 1e-12 reduction-order noise across thread counts. +def test_pure_eb_reproduces_pins_on_committed_xi(pure_eb_xi): + """The pure-E/B estimator on the committed ξ± reproduces its pins. + + With ξ± frozen, these pins move only when the estimator does. The operator + is deterministic linear algebra, so rtol=1e-8 leaves room only for BLAS + summation order. """ - modes = b_modes.pure_eb_from_xi(**pure_eb_xi) + n_fine = len(pure_eb_xi["theta_int"]) + results = b_modes.calculate_pure_eb_correlation( + **pure_eb_xi, cov_xi=np.eye(2 * n_fine) + ) for key in b_modes._EB_KEYS: - npt.assert_allclose(modes[key], _PURE_EB_PINS[key], rtol=1e-6, err_msg=key) + npt.assert_allclose(results[key], _PURE_EB_PINS[key], rtol=1e-8, err_msg=key) - # Teeth: widening the integration interval by 1% leaves the pins. - moved = b_modes.pure_eb_from_xi(**{**pure_eb_xi, "tmax": 1.01 * pure_eb_xi["tmax"]}) + # Teeth: uniform rather than pair weights leave the pins. + moved = b_modes.calculate_pure_eb_correlation( + **{**pure_eb_xi, "weight_int": np.ones(n_fine)}, cov_xi=np.eye(2 * n_fine) + ) assert not np.allclose(moved["xip_E"], _PURE_EB_PINS["xip_E"], rtol=1e-6, atol=0) +def test_pure_eb_modes_sum_to_the_averaged_xi(pure_eb_xi): + """ξ± = E ± B + amb holds in every reporting bin. + + It holds at each fine node by construction of the decomposition, and the + modes and the reported ξ± are the same pair-weighted average of those nodes. + """ + n_fine = len(pure_eb_xi["theta_int"]) + r = b_modes.calculate_pure_eb_correlation(**pure_eb_xi, cov_xi=np.eye(2 * n_fine)) + scale = np.abs(r["xip"]).max() + npt.assert_allclose( + r["xip"], r["xip_E"] + r["xip_B"] + r["xip_amb"], rtol=0, atol=1e-12 * scale + ) + npt.assert_allclose( + r["xim"], r["xim_E"] - r["xim_B"] + r["xim_amb"], rtol=0, atol=1e-12 * scale + ) + + +def test_pure_eb_covariance_is_the_operator_sandwich(pure_eb_xi): + """``cov`` is K C Kᵀ for the supplied ξ± covariance, and records npatch. + + The reported ξ± variances are the same pair-weighted average pushed + through the ξ+ and ξ− blocks of C, and the reported edges are the snapped + ones the operator used. + """ + n_fine = len(pure_eb_xi["theta_int"]) + cov_xi = _spd(2 * n_fine, seed=7) + r = b_modes.calculate_pure_eb_correlation(**pure_eb_xi, cov_xi=cov_xi, npatch=40) + K, P, edges = b_modes.pure_eb_operator( + pure_eb_xi["weight_int"], + pure_eb_xi["edges_int"], + pure_eb_xi["reporting_edges"], + ) + npt.assert_allclose(r["cov"], K @ cov_xi @ K.T, rtol=1e-12) + npt.assert_allclose(r["cov"], r["cov"].T, rtol=1e-12) + npt.assert_allclose(r["var_xip"], np.diag(P @ cov_xi[:n_fine, :n_fine] @ P.T)) + npt.assert_allclose(r["var_xim"], np.diag(P @ cov_xi[n_fine:, n_fine:] @ P.T)) + npt.assert_array_equal(r["left_edges"], edges[:-1]) + npt.assert_array_equal(r["right_edges"], edges[1:]) + assert np.all(np.isin(edges, pure_eb_xi["edges_int"])) + assert r["npatch"] == 40 + + with pytest.raises(ValueError, match="npatch > 1"): + b_modes.calculate_pure_eb_correlation(**pure_eb_xi, cov_xi=cov_xi, npatch=1) + + +def test_pure_eb_reporting_bins_are_unions_of_fine_bins(): + """Requested edges snap to the nearest fine edge; rows pool whole fine bins. + + Rows sum to one and weight their fine bins by the pair weight; a fine bin + with no weight carries none; equal weights give the plain mean. Edges that + snap together, edges outside the fine grid and empty bins raise. + """ + edges_int = np.geomspace(1.0, 100.0, 41) + weight = np.arange(1.0, 41.0) + weight[20] = 0.0 + requested = np.array([2.1, 6.0, 20.0, 49.0]) + P, edges = b_modes._reporting_binning(weight, edges_int, requested) + + snap = [np.argmin(np.abs(np.log(edges_int / e))) for e in requested] + npt.assert_array_equal(edges, edges_int[snap]) + npt.assert_allclose(P.sum(axis=1), 1.0) + for row, (lo, hi) in enumerate(zip(snap[:-1], snap[1:])): + inside = np.zeros(40, dtype=bool) + inside[lo:hi] = True + npt.assert_allclose(P[row, inside], weight[inside] / weight[inside].sum()) + assert np.all(P[row, ~inside] == 0) + assert np.all(P[:, 20] == 0) + + flat, _ = b_modes._reporting_binning(np.ones(40), edges_int, requested) + npt.assert_allclose(flat[0, snap[0] : snap[1]], 1.0 / (snap[1] - snap[0])) + + with pytest.raises(ValueError, match="same fine edge"): + b_modes._reporting_binning(weight, edges_int, [2.0, 2.05, 20.0]) + with pytest.raises(ValueError, match="outside the fine grid"): + b_modes._reporting_binning(weight, edges_int, [2.0, 20.0, 200.0]) + with pytest.raises(ValueError, match="hold no integration-grid pairs"): + b_modes._reporting_binning(np.zeros(40), edges_int, requested) + + +def _weighted_catalogue(treecorr): + """A weighted flat shear catalogue, separations in arcmin.""" + rng = np.random.default_rng(2024) + n_gal = 3000 + x, y = rng.uniform(0.0, 300.0, (2, n_gal)) + g1, g2 = 0.02 + 0.05 * rng.standard_normal((2, n_gal)) + return treecorr.Catalog(x=x, y=y, g1=g1, g2=g2, w=rng.uniform(0.2, 1.0, n_gal)) + + +@pytest.mark.parametrize( + "requested", + [np.geomspace(15.0, 70.0, 7), np.geomspace(14.0, 73.0, 7)], + ids=["nested", "snapped"], +) +def test_pure_eb_binning_reproduces_the_reporting_measurement(requested): + """P on the fine grid gives TreeCorr's ξ± and meanr on the snapped bins. + + TreeCorr's ξ± and meanr in a bin are averages over its pairs weighted by + ``w_i w_j``, and every reporting bin is a union of fine bins, so pooling + them with the fine ``weight`` is the same sum. A weighted catalogue with + exact binning makes that hold to round-off, whether or not the requested + edges fall on fine edges; weighting by ``npairs`` instead does not. + """ + treecorr = pytest.importorskip("treecorr") + cat = _weighted_catalogue(treecorr) + exact = {"bin_slop": 0, "angle_slop": 0} + + # Eight fine bins per [15, 70]′ / 6 reporting bin, plus four either side. + step = np.log(70.0 / 15.0) / 48 + fine = treecorr.GGCorrelation( + min_sep=15.0 * np.exp(-4 * step), + max_sep=70.0 * np.exp(4 * step), + nbins=56, + **exact, + ) + fine.process(cat) + edges_int = np.append(fine.left_edges, fine.right_edges[-1]) + + P, edges = b_modes._reporting_binning(fine.weight, edges_int, requested) + assert np.all(np.isin(edges, edges_int)) + measured = [] + for lo, hi in zip(edges[:-1], edges[1:]): + gg = treecorr.GGCorrelation(min_sep=lo, max_sep=hi, nbins=1, **exact) + gg.process(cat) + measured.append(gg) + for key in ("xip", "xim", "meanr"): + npt.assert_allclose( + P @ getattr(fine, key), + [getattr(gg, key)[0] for gg in measured], + rtol=1e-10, + err_msg=key, + ) + + by_npairs, _ = b_modes._reporting_binning(fine.npairs, edges_int, requested) + assert not np.allclose( + by_npairs @ fine.xip, [gg.xip[0] for gg in measured], rtol=1e-6, atol=0 + ) + + +def test_pure_eb_transform_needs_a_log_uniform_grid(): + """The transform runs on log-uniform nodes and refuses irregular edges. + + cosmo_numba's interpolator places samples on a regular grid in log θ, so + edges that are not log-uniform would silently mis-place them. + """ + edges = np.geomspace(1.0, 300.0, 61) + edges[30] *= 1.001 + with pytest.raises(ValueError, match="log-uniform"): + b_modes._fixed_quadrature_operator(edges) + + +def test_pure_eb_operator_refuses_a_reporting_floor_at_the_grid_edge(pure_eb_xi): + """Reporting bins that reach the fine grid's floor raise, never return NaN. + + The ξ− integrals over [tmin, t] have fewer than interp_order + 1 nodes for + the first few fine nodes, so those operator rows are under-determined. + """ + edges_int = pure_eb_xi["edges_int"] + with pytest.raises(ValueError, match="under-determined"): + b_modes.pure_eb_operator( + pure_eb_xi["weight_int"], + edges_int, + np.geomspace(edges_int[0], 70.0, 7), + ) + + # --------------------------------------------------------------------------- # 6. Grid edges and the COSEBIs covariance seam # --------------------------------------------------------------------------- +def test_scale_cuts_select_the_same_bins_on_every_fine_grid(): + """A cut on nominal edges picks the same bins whatever grid they snapped to. + + [1, 250]′ in 20 bins snapped onto the cosmo_val (0.08–300′) and Paper II + (0.5–300′) fine grids moves edge 9 (11.997′) to either side of 12′, so + exact containment would disagree; snapping the cut to the nearest edge + selects bins 9–15 for [12, 83]′ and all bins for [1, 250]′ on both. + """ + requested = np.geomspace(1.0, 250.0, 21) + for lo in (0.08, 0.5): + fine = np.geomspace(lo, 300.0, 1001) + _, edges = b_modes._reporting_binning(np.ones(1000), fine, requested) + left, right = edges[:-1], edges[1:] + assert b_modes.bins_from_scale_cut(left, right, (12.0, 83.0)) == (9, 16) + assert b_modes.bins_from_scale_cut(left, right, (1.0, 250.0)) == (0, 20) + pte = np.arange(400.0).reshape(20, 20) + assert ( + b_modes._get_pte_from_scale_cut(pte, (left, right), (12, 83)) == pte[9, 15] + ) + + with pytest.raises(RuntimeError, match="selects no bins"): + b_modes.bins_from_scale_cut(left, right, (12.0, 12.5)) + + def test_log_bin_edges_matches_the_grid_stub(): """Edges reconstructed from a binning are the ones TreeCorr would report. @@ -477,13 +701,17 @@ def test_pure_eb_npz_carries_what_the_summary_reads(tmp_path): """The .npz keys cv_summarize_bmodes reads are the ones the writer emits. The two live in different rules, so the contract between them — the PTE - matrices under ``pte_matrices_{stat}`` and the realisation count under - ``n_eff`` — is pinned here rather than discovered on a cluster run. + matrices under ``pte_matrices_{stat}``, the reporting edges they are + indexed on and the jackknife patch count under ``npatch`` — is pinned here + rather than discovered on a cluster run. """ results, nbins = _eb_inputs() results.update( {key: np.zeros(nbins) for key in b_modes._EB_KEYS if key not in results} ) + results["left_edges"], results["right_edges"] = b_modes.log_bin_edges( + 1.0, 100.0, nbins + ) results = b_modes.calculate_eb_statistics(results) out = tmp_path / "pure_eb_data.npz" @@ -493,64 +721,17 @@ def test_pure_eb_npz_carries_what_the_summary_reads(tmp_path): for stat in ("xip_B", "xim_B", "combined"): assert f"pte_matrices_{stat}" in saved assert saved[f"pte_matrices_{stat}"].shape == (nbins, nbins) - assert saved["n_eff"] == results["n_eff"] + assert saved["npatch"] == results["npatch"] npt.assert_allclose(saved["theta"], results["theta"]) for key in b_modes._EB_KEYS: assert key in saved # The summary reads the fiducial cut out of those matrices through the same - # helper the plots use, so a valid cut must resolve to a finite PTE. - edges = b_modes.log_bin_edges(1.0, 100.0, nbins) + # helper the plots use, on the saved edges, so a valid cut must resolve to + # a finite PTE. + edges = (saved["left_edges"], saved["right_edges"]) + npt.assert_array_equal(edges[0], results["left_edges"]) pte = b_modes._get_pte_from_scale_cut( saved["pte_matrices_xip_B"], edges, (1.0, 100.0) ) assert np.isfinite(pte) - - -def test_pure_eb_covariance_mc_draws_around_the_theory_mean(monkeypatch): - """The MC draws centre on cs_util's theory ξ±, binned to the reporting grid. - - ``get_theo_xi`` returns ``{pair: (ξ+, ξ−)}``; one n(z) is one pair. With a - zero covariance every draw is the theory mean, so a stub kernel that echoes - its reporting-grid ξ± pins the unpack, the [ξ+; ξ−] order and the binning. - """ - nrep = 4 - left, right = b_modes.log_bin_edges(2.0, 50.0, nrep) - theta = np.sqrt(left * right) - theta_int = np.geomspace(1.0, 100.0, 40) - xip_th, xim_th = theta_int.copy(), 2 * theta_int - - monkeypatch.setattr( - b_modes, "get_theo_xi", lambda **kw: {"W0xW0": (xip_th, xim_th)} - ) - - def _echo(theta, theta_int, xip, xim, xip_int, xim_int, tmin, tmax, parallel): - zeros = np.zeros_like(xip) - return xip, xim, zeros, zeros, zeros, zeros - - module = types.ModuleType("cosmo_numba.B_modes.schneider2022") - module.get_pure_EB_modes = _echo - monkeypatch.setitem( - __import__("sys").modules, "cosmo_numba.B_modes.schneider2022", module - ) - - cov, eb_samples = b_modes.pure_eb_covariance_mc( - theta=theta, - left_edges=left, - right_edges=right, - theta_int=theta_int, - cov_int=np.zeros((2 * len(theta_int), 2 * len(theta_int))), - z=np.linspace(0.01, 2.0, 50), - nz=np.ones(50), - cosmo=None, - n_samples=3, - ) - - inside = [(theta_int >= lo) & (theta_int < hi) for lo, hi in zip(left, right)] - expected_xip = np.array([xip_th[m].mean() for m in inside]) - expected_xim = np.array([xim_th[m].mean() for m in inside]) - assert eb_samples.shape == (3, 6 * nrep) - for draw in eb_samples: - npt.assert_allclose(draw[:nrep], expected_xip) - npt.assert_allclose(draw[nrep : 2 * nrep], expected_xim) - npt.assert_allclose(cov, 0.0, atol=1e-20) diff --git a/src/sp_validation/tests/test_cosmo_val.py b/src/sp_validation/tests/test_cosmo_val.py index 54aed4a5..c08c055e 100644 --- a/src/sp_validation/tests/test_cosmo_val.py +++ b/src/sp_validation/tests/test_cosmo_val.py @@ -581,17 +581,13 @@ def test_calculate_scale_dependent_leakage_runs_on_synthetic_catalog( assert hasattr(res, "C_sys_p") and hasattr(res, "C_sys_m") def test_calculate_pure_eb_runs_on_synthetic_catalog(self, tmp_path, pure_eb_xi): - """calculate_pure_eb carries ξ± through cosmo_numba's pure-E/B split. + """calculate_pure_eb measures the fine ξ± and pushes it through the operator. - The ξ± it measures equal the committed ``pure_eb_xi``, its modes are - ``pure_eb_from_xi`` of those ξ± and edges, and every reporting bin is - finite. ``test_b_modes`` pins the transform itself on the same ξ±, so a - failure names the step that moved: measurement, wiring or transform. - - Finiteness: the Schneider (2022) integrals are near-singular where a - reporting bin meets the integration boundary, so the integration grid - [1, 300]′ brackets the reporting grid [15, 70]′ on both ends and is fine - (600 bins); about 80 integration bins NaN the edge bins. + The ξ± it measures equal the committed ``pure_eb_xi`` (``test_b_modes`` + pins the operator on the same ξ±, so a failure names the step that + moved), and the jackknife covariance of the modes is the jackknife ξ± + covariance through the operator: for a linear estimator that equals + TreeCorr's per-patch jackknife of the modes themselves. ξ±: exact binning (bin_slop = angle_slop = 0) makes ξ± a plain pair sum, independent of the tree and so of the jackknife patches, whose k-means @@ -599,6 +595,8 @@ def test_calculate_pure_eb_runs_on_synthetic_catalog(self, tmp_path, pure_eb_xi) """ pytest.importorskip("treecorr") pytest.importorskip("cosmo_numba") + import treecorr + from sp_validation import b_modes # Coherent shear -> smooth xi+/-, so the pure-E/B integral is well-posed. @@ -627,29 +625,39 @@ def test_calculate_pure_eb_runs_on_synthetic_catalog(self, tmp_path, pure_eb_xi) ) measured = { - "theta_report": results["theta"], - "xip_report": results["xip"], - "xim_report": results["xim"], - "theta_int": results["theta_int"], - "xip_int": results["xip_int"], - "xim_int": results["xim_int"], - "tmin": results["left_edges"][0], - "tmax": results["right_edges"][-1], + key: results[key] + for key in ("theta_int", "xip_int", "xim_int", "weight_int", "edges_int") } + measured["reporting_edges"] = np.geomspace(15.0, 70.0, nbins + 1) # Regenerate the fixture with np.savez(conftest.PURE_EB_XI, **measured). for key, value in measured.items(): np.testing.assert_allclose( value, pure_eb_xi[key], rtol=1e-10, atol=0, err_msg=key ) - modes = b_modes.pure_eb_from_xi(**measured) - for key in b_modes._EB_KEYS: + operator, _, edges = b_modes.pure_eb_operator( + *(measured[k] for k in ("weight_int", "edges_int", "reporting_edges")) + ) + np.testing.assert_array_equal(results["left_edges"], edges[:-1]) + modes = operator @ np.concatenate([measured["xip_int"], measured["xim_int"]]) + for i, key in enumerate(b_modes._EB_KEYS): vec = np.asarray(results[key]) assert vec.shape == (nbins,) assert np.all(np.isfinite(vec)), f"{key} not finite" - np.testing.assert_allclose(vec, modes[key], rtol=1e-10, err_msg=key) + np.testing.assert_allclose( + vec, modes[i * nbins : (i + 1) * nbins], rtol=1e-12, err_msg=key + ) - # Jackknife covariance over the 6 stats (xip/xim x E/B/amb) x nbins. cov = np.asarray(results["cov"]) assert cov.shape == (6 * nbins, 6 * nbins) - assert results["n_eff"] == npatch + assert results["npatch"] == npatch + gg_int = cv.cat_ggs[version] + jackknife_of_modes = treecorr.estimate_multi_cov( + [gg_int], + "jackknife", + func=lambda corrs: operator @ np.concatenate([corrs[0].xip, corrs[0].xim]), + cross_patch_weight="match", + ) + np.testing.assert_allclose( + cov, jackknife_of_modes, rtol=0, atol=1e-10 * np.abs(cov).max() + ) diff --git a/src/sp_validation/tests/test_sacc_io.py b/src/sp_validation/tests/test_sacc_io.py index 8f8564ae..872a8198 100644 --- a/src/sp_validation/tests/test_sacc_io.py +++ b/src/sp_validation/tests/test_sacc_io.py @@ -100,6 +100,7 @@ def test_xi_roundtrip(tmp_path): assert np.array_equal(th, theta) assert np.array_equal(p, xip) assert np.array_equal(m, xim) + assert np.array_equal(sio.get_xi_weight(s2, (0, 0), grid="reporting"), weight) # extra tags survive idx = s2.indices(sio.XI_PLUS, ("source_0", "source_0"), grid="reporting") tags = s2.data[idx[0]].tags diff --git a/uv.lock b/uv.lock index 94341fce..8ab6ed4c 100644 --- a/uv.lock +++ b/uv.lock @@ -588,8 +588,8 @@ wheels = [ [[package]] name = "cosmo-numba" -version = "1.0.1.dev6+g188d272c6" -source = { git = "https://github.com/aguinot/cosmo-numba.git?rev=main#188d272c67d6d699d7cff7a7a53ded8dec759bb0" } +version = "0.1.dev108+ga64cb2ed1" +source = { git = "https://github.com/cailmdaley/cosmo-numba.git?rev=a64cb2ed13e595a16ad056a128ef24594b811504#a64cb2ed13e595a16ad056a128ef24594b811504" } dependencies = [ { name = "mpmath" }, { name = "numba" }, @@ -3821,7 +3821,7 @@ requires-dist = [ { name = "camb", specifier = ">=2.0" }, { name = "clmm" }, { name = "colorama" }, - { name = "cosmo-numba", git = "https://github.com/aguinot/cosmo-numba.git?rev=main" }, + { name = "cosmo-numba", git = "https://github.com/cailmdaley/cosmo-numba.git?rev=a64cb2ed13e595a16ad056a128ef24594b811504" }, { name = "cosmology", marker = "extra == 'glass'", specifier = "==2022.10.9" }, { name = "cosmosis", marker = "extra == 'workflow'", specifier = ">=3.25" }, { name = "cryptography" }, diff --git a/workflow/rules/cosmo_val.smk b/workflow/rules/cosmo_val.smk index 1bff5b43..8d9839c6 100644 --- a/workflow/rules/cosmo_val.smk +++ b/workflow/rules/cosmo_val.smk @@ -7,8 +7,7 @@ # products they write under COSMO_VAL (= cosmo_val/output): # # catalogue ──→ xi (one job per grid: reporting, integration) -# ├─ reporting part ──────→ pure_eb (part, npz, figures) -# ├─ integration part ─┬──→ pure_eb +# ├─ integration part ─┬──→ pure_eb (part, npz, figures) # │ └──→ cosebis (part, npz, figures) # └─ reporting .txt ──→ 2pcf plot, ratio_xi_sys_xi # CosmoCov ξ±, integration grid ──→ pure_eb, cosebis (their covariances) @@ -80,7 +79,7 @@ def _pure_eb_stub(version): f"{version}_eb_minsep={CV['theta_min']}_maxsep={CV['theta_max']}" f"_nbins={CV['nbins']}_minsepint={eb['min_sep']}" f"_maxsepint={eb['max_sep']}_nbinsint={eb['nbins']}" - f"_npatch={CV['npatch']}_varmethod=semi-analytic" + f"_npatch={CV['npatch']}_varmethod=analytic" ) ) @@ -396,13 +395,12 @@ rule cv_plot_pseudo_cl: # --------------------------------------------------------------------------- rule cv_pure_eb: - """Pure E/B-mode decomposition for one version, from its ξ± parts. + """Pure E/B-mode decomposition for one version, from its integration-grid part. - The modes come from the two parts; the covariance is Monte Carlo from the - integration-grid covariance model, so no patched estimator run is involved. + The modes are a fixed linear operator on the part's ξ±; the covariance is + the integration-grid ξ± covariance pushed exactly through it. """ input: - xi_reporting=lambda w: cv_xi_sacc(w.version, "reporting"), xi_integration=lambda w: cv_xi_sacc(w.version, "integration"), cov_integration=lambda w: cv_xi_cov_integration(w.version), output: @@ -414,13 +412,11 @@ rule cv_pure_eb: min_sep=CV["theta_min"], max_sep=CV["theta_max"], nbins=CV["nbins"], - n_samples=CV.get("n_mc_samples", 1000), - cosmo_params=CV["cosmo_params"], + integration=XI_GRIDS["integration"], fiducial_scale_cut=CV["fiducial_scale_cut"], - threads: 24 resources: - mem_mb=40000, - runtime=360, + mem_mb=8000, + runtime=20, script: "../scripts/cv_pure_eb.py" @@ -472,9 +468,6 @@ rule cv_summarize_bmodes: params: versions=CV_VERSIONS, fiducial_scale_cut=CV["fiducial_scale_cut"], - min_sep=CV["theta_min"], - max_sep=CV["theta_max"], - nbins=CV["nbins"], include_pseudo_cl=CV.get("include_pseudo_cl", False), resources: mem_mb=8000, diff --git a/workflow/scripts/cv_pure_eb.py b/workflow/scripts/cv_pure_eb.py index 1edaea93..3f978d3a 100644 --- a/workflow/scripts/cv_pure_eb.py +++ b/workflow/scripts/cv_pure_eb.py @@ -1,28 +1,25 @@ """Rule cv_pure_eb: pure E/B-mode decomposition for one version. -A consumer of the two ξ± parts plus one covariance file — nothing here touches -a catalogue. The modes come from the reporting and integration parts through -the pipeline kernel; the covariance is Monte Carlo through that same kernel, -drawn from the CosmoCov integration-grid ξ± covariance around a theory mean, so -it depends on the covariance model and the grids rather than on the measured -vector. A jackknife of the transformed modes would need per-patch realisations, -which are never persisted. +A consumer of the integration-grid ξ± part plus a ξ± covariance on that grid — +nothing here touches a catalogue. The estimator is one fixed linear operator +on the fine ξ± (b_modes.pure_eb_operator), averaged into the reporting bins with +the part's TreeCorr pair weights into reporting bins snapped onto the fine +edges, so its covariance is the supplied ξ± covariance pushed +exactly through that operator. """ import numpy as np -from cs_util.cosmo import get_cosmo from cv_runner import _unbuffer_streams, verify_outputs from sp_validation import sacc_io from sp_validation.b_modes import ( calculate_eb_statistics, - log_bin_edges, + calculate_pure_eb_correlation, + covariance_label, plot_eb_covariance_matrix, plot_integration_vs_reporting, plot_pte_2d_heatmaps, plot_pure_eb_correlations, - pure_eb_covariance_mc, - pure_eb_from_xi, save_pure_eb_results, ) from sp_validation.cosmo_val.sacc_writers import pure_eb_to_sacc @@ -32,48 +29,22 @@ version = p["version"] fiducial_scale_cut = tuple(p["fiducial_scale_cut"]) -reporting = sacc_io.load(snakemake.input["xi_reporting"]) -integration = sacc_io.load(snakemake.input["xi_integration"]) -theta, xip, xim = sacc_io.get_xi(reporting, (0, 0), grid="reporting") -theta_int, xip_int, xim_int = sacc_io.get_xi(integration, (0, 0), grid="integration") -left_edges, right_edges = log_bin_edges(p["min_sep"], p["max_sep"], p["nbins"]) +part = sacc_io.load(snakemake.input["xi_integration"]) +theta_int, xip_int, xim_int = sacc_io.get_xi(part, (0, 0), grid="integration") +weight_int = sacc_io.get_xi_weight(part, (0, 0), grid="integration") +# A part stores bin centres; the edges come from the grid it was measured on. +grid = p["integration"] +edges_int = np.geomspace(grid["min_sep"], grid["max_sep"], grid["nbins"] + 1) -# The reporting grid must sit strictly inside the integration grid: a reporting -# point on the boundary has no interior support and comes back NaN. -modes = pure_eb_from_xi( - theta, xip, xim, theta_int, xip_int, xim_int, left_edges[0], right_edges[-1] +results = calculate_pure_eb_correlation( + theta_int, + xip_int, + xim_int, + weight_int, + edges_int, + np.loadtxt(snakemake.input["cov_integration"]), + np.geomspace(p["min_sep"], p["max_sep"], p["nbins"] + 1), ) - -z, nz = sacc_io.get_nz(reporting, 0) -cov, eb_samples = pure_eb_covariance_mc( - theta=theta, - left_edges=left_edges, - right_edges=right_edges, - theta_int=theta_int, - cov_int=np.loadtxt(snakemake.input["cov_integration"]), - z=z, - nz=nz, - cosmo=get_cosmo(**p["cosmo_params"]), - n_samples=p["n_samples"], -) - -variances = reporting.covariance.dense.diagonal() -results = { - "theta": theta, - "left_edges": left_edges, - "right_edges": right_edges, - "xip": xip, - "xim": xim, - "var_xip": variances[: len(theta)], - "var_xim": variances[len(theta) :], - "theta_int": theta_int, - "xip_int": xip_int, - "xim_int": xim_int, - "n_eff": p["n_samples"], - "cov": cov, - "eb_samples": eb_samples, - **modes, -} results = calculate_eb_statistics(results) plot_integration_vs_reporting( @@ -94,19 +65,22 @@ fiducial_xim_scale_cut=fiducial_scale_cut, ) plot_eb_covariance_matrix( - cov, "semi-analytic", snakemake.output["figure_covariance"], version + results["cov"], + covariance_label(results["npatch"]), + snakemake.output["figure_covariance"], + version, ) save_pure_eb_results(results, snakemake.output["npz"]) # The part inherits the ξ± part's provenance; `type` is re-stamped on save. -metadata = {k: v for k, v in reporting.metadata.items() if k != "type"} +metadata = {k: v for k, v in part.metadata.items() if k != "type"} s = pure_eb_to_sacc( - {0: (z, nz)}, + {0: sacc_io.get_nz(part, 0)}, metadata, - theta, + results["theta"], {key: results[key] for key in sacc_io.PURE_KEYS}, - covariance=cov, + covariance=results["cov"], ) sacc_io.save(s, snakemake.output["sacc"], type="data") diff --git a/workflow/scripts/cv_summarize_bmodes.py b/workflow/scripts/cv_summarize_bmodes.py index ec44142a..3f9d9e82 100644 --- a/workflow/scripts/cv_summarize_bmodes.py +++ b/workflow/scripts/cv_summarize_bmodes.py @@ -13,15 +13,13 @@ from cv_runner import _unbuffer_streams, verify_outputs from sp_validation import sacc_io -from sp_validation.b_modes import _get_pte_from_scale_cut, log_bin_edges +from sp_validation.b_modes import _get_pte_from_scale_cut, covariance_label from sp_validation.cosmo_val.core import print_bmode_summary from sp_validation.statistics import chi2_and_pte _unbuffer_streams() p = snakemake.params fiducial_scale_cut = tuple(p["fiducial_scale_cut"]) -edges = log_bin_edges(p["min_sep"], p["max_sep"], p["nbins"]) - summary = {} cov_methods = set() @@ -29,6 +27,8 @@ row = {} pure_eb = np.load(snakemake.input["pure_eb"][i]) + # The bins the PTE matrices are indexed on: the snapped reporting edges. + edges = (pure_eb["left_edges"], pure_eb["right_edges"]) for stat in ("xip_B", "xim_B", "combined"): try: row[stat] = _get_pte_from_scale_cut( @@ -36,7 +36,8 @@ ) except (KeyError, RuntimeError): pass - cov_methods.add(f"pure-E/B: semi-analytic ({int(pure_eb['n_eff'])} draws)") + npatch = int(pure_eb["npatch"]) if "npatch" in pure_eb else None + cov_methods.add(f"pure-E/B: {covariance_label(npatch)}") # The COSEBIs .npz is written at the fiducial cut, so its PTE is the one # this table wants. diff --git a/workflow/tests/test_dag.py b/workflow/tests/test_dag.py index 90272dca..b11acf3b 100644 --- a/workflow/tests/test_dag.py +++ b/workflow/tests/test_dag.py @@ -52,7 +52,7 @@ def test_one_integration_grid(toy): part = str(toy.cosmo_val / f"{version}_xi_{tag}.sacc") covariance = str(toy.covariances[version, "g"]) assert by_rule["cv_cosebis"] == {part, covariance}, by_rule["cv_cosebis"] - assert {part, covariance} <= by_rule["cv_pure_eb"], by_rule["cv_pure_eb"] + assert by_rule["cv_pure_eb"] == {part, covariance}, by_rule["cv_pure_eb"] @pytest.mark.parametrize("named", [True, False], ids=["named", "unnamed"])