From 37ac49aa6804aeb80edc79f881f3ccc5a1fcca8b Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 24 Oct 2026 11:56:37 -0400 Subject: [PATCH 01/18] Split E/B and COSEBIs rendering from sealed transform products --- src/sp_validation/b_modes.py | 53 +++++++++++++++++++++++++---- workflow/rules/cosmo_val.smk | 30 ++++++++++++++-- workflow/scripts/cv_cosebis.py | 22 +----------- workflow/scripts/cv_plot_cosebis.py | 21 ++++++++++++ workflow/scripts/cv_plot_pure_eb.py | 26 ++++++++++++++ workflow/scripts/cv_pure_eb.py | 29 ---------------- 6 files changed, 123 insertions(+), 58 deletions(-) create mode 100644 workflow/scripts/cv_plot_cosebis.py create mode 100644 workflow/scripts/cv_plot_pure_eb.py diff --git a/src/sp_validation/b_modes.py b/src/sp_validation/b_modes.py index 083dafd1..ea9a80ca 100644 --- a/src/sp_validation/b_modes.py +++ b/src/sp_validation/b_modes.py @@ -1310,12 +1310,13 @@ def save_pure_eb_results(results, output_path): output_path : str Output .npz file path """ - # Data vectors and covariance + # Include the sealed fine/reporting vectors and variances used by the + # comparison figures, not just the transformed modes. save_dict = { - key: results[key] for key in ("theta", "left_edges", "right_edges", "cov") + key: value + for key, value in results.items() + if key not in ("pte_matrices", "npatch") and value is not None } - for key in _EB_KEYS: - save_dict[key] = results[key] # PTE matrices for key, matrix in results.get("pte_matrices", {}).items(): @@ -1329,6 +1330,40 @@ def save_pure_eb_results(results, output_path): print(f"Saved pure E/B results to {output_path}") +def load_pure_eb_results(path): + """Read a sealed plot cache without transforms or catalogue access.""" + with np.load(path, allow_pickle=False) as data: + results = {k: data[k] for k in data.files if not k.startswith("pte_matrices_")} + results["pte_matrices"] = { + k.removeprefix("pte_matrices_"): data[k] + for k in data.files + if k.startswith("pte_matrices_") + } + results["npatch"] = int(results["npatch"]) if "npatch" in results else None + return results + + +def load_cosebis_results(path): + """Read all cached scale cuts and their integration-bin edges.""" + with np.load(path, allow_pickle=False) as data: + results = {} + for lo, hi in data["scale_cuts"]: + suffix = f"_{lo}_{hi}" + # Match the writer's scalar formatting (including integer cuts). + suffix = next( + k.removeprefix("En") + for k in data.files + if k.startswith("En_") + and tuple(map(float, k[3:].split("_"))) == (lo, hi) + ) + results[lo, hi] = { + k: data[k + suffix] + for k in ("En", "Bn", "cov", "chi2_E", "chi2_B", "pte_B") + } + edges = data["left_edges"], data["right_edges"] + return results, edges + + def _cosebis_result_to_dict(r, suffix=""): """Build save dict entries from a single COSEBIs result.""" return { @@ -1341,7 +1376,7 @@ def _cosebis_result_to_dict(r, suffix=""): } -def save_cosebis_results(results, output_path, fiducial_scale_cut=None): +def save_cosebis_results(results, output_path, fiducial_scale_cut=None, edges=None): """ Save COSEBIs data vectors and covariance to .npz. @@ -1359,10 +1394,14 @@ def save_cosebis_results(results, output_path, fiducial_scale_cut=None): is_multi = all(isinstance(k, tuple) for k in results) if is_multi and fiducial_scale_cut is not None: - # Select the single best-matching scale cut + # Keep the untagged fiducial entries for summary consumers, and every + # scanned cut for the plot-only heatmap. key = find_conservative_scale_cut_key(results, fiducial_scale_cut) save_dict = _cosebis_result_to_dict(results[key]) save_dict["scale_cut"] = np.array(key) + for sc, r in results.items(): + save_dict.update(_cosebis_result_to_dict(r, f"_{sc[0]}_{sc[1]}")) + save_dict["scale_cuts"] = np.array(list(results.keys())) elif is_multi: # Save all scale cuts with tagged keys save_dict = {} @@ -1375,6 +1414,8 @@ def save_cosebis_results(results, output_path, fiducial_scale_cut=None): if "scale_cut" in results: save_dict["scale_cut"] = np.array(results["scale_cut"]) + if edges is not None: + save_dict["left_edges"], save_dict["right_edges"] = edges np.savez(output_path, **save_dict) print(f"Saved COSEBIs results to {output_path}") diff --git a/workflow/rules/cosmo_val.smk b/workflow/rules/cosmo_val.smk index b416eac1..c6078f14 100644 --- a/workflow/rules/cosmo_val.smk +++ b/workflow/rules/cosmo_val.smk @@ -388,7 +388,6 @@ rule cv_pure_eb: output: npz=cv_pure_eb_npz("{version}"), sacc=cv_pure_eb_sacc("{version}"), - **cv_pure_eb_figures("{version}"), params: version="{version}", blind=lambda w: blind_of(w.version), @@ -416,7 +415,6 @@ rule cv_cosebis: output: npz=cv_cosebis_npz("{version}"), sacc=cv_cosebis_sacc("{version}"), - **cv_cosebis_figures("{version}"), params: version="{version}", blind=lambda w: blind_of(w.version), @@ -434,6 +432,34 @@ rule cv_cosebis: "../scripts/cv_cosebis.py" +rule cv_plot_pure_eb: + """Render pure E/B figures without rerunning the transform.""" + input: + npz=lambda w: cv_pure_eb_npz(w.version), + plot_code=str(REPO_ROOT / "src/sp_validation/b_modes.py"), + output: + **cv_pure_eb_figures("{version}"), + params: + version="{version}", + fiducial_scale_cut=CV["fiducial_scale_cut"], + script: + "../scripts/cv_plot_pure_eb.py" + + +rule cv_plot_cosebis: + """Render COSEBIs figures without rerunning the transform.""" + input: + npz=lambda w: cv_cosebis_npz(w.version), + plot_code=str(REPO_ROOT / "src/sp_validation/b_modes.py"), + output: + **cv_cosebis_figures("{version}"), + params: + version="{version}", + fiducial_scale_cut=CV["fiducial_scale_cut"], + script: + "../scripts/cv_plot_cosebis.py" + + rule cv_summarize_bmodes: """Collect B-mode PTEs (pure E/B, COSEBIs, pseudo-Cl) into one summary.""" input: diff --git a/workflow/scripts/cv_cosebis.py b/workflow/scripts/cv_cosebis.py index 0223de82..97308d4d 100644 --- a/workflow/scripts/cv_cosebis.py +++ b/workflow/scripts/cv_cosebis.py @@ -15,9 +15,6 @@ cosebis_scan_from_xi, find_conservative_scale_cut_key, log_bin_edges, - plot_cosebis_covariance_matrix, - plot_cosebis_modes, - plot_cosebis_scale_cut_heatmap, save_cosebis_results, ) from sp_validation.cosmo_val.sacc_writers import cosebis_to_sacc @@ -44,24 +41,7 @@ fiducial_key = find_conservative_scale_cut_key(results, fiducial_scale_cut) fiducial = results[fiducial_key] -plot_cosebis_modes( - fiducial, - version, - snakemake.output["figure_modes"], - fiducial_scale_cut=fiducial_scale_cut, -) -plot_cosebis_covariance_matrix( - fiducial, version, "analytic", snakemake.output["figure_covariance"] -) -plot_cosebis_scale_cut_heatmap( - results, - edges, - version, - snakemake.output["figure_scalecut_ptes"], - fiducial_scale_cut=fiducial_scale_cut, -) - -save_cosebis_results(results, snakemake.output["npz"], fiducial_scale_cut) +save_cosebis_results(results, snakemake.output["npz"], fiducial_scale_cut, edges) s = cosebis_to_sacc({0: sacc_io.get_nz(part, 0)}, part.metadata, fiducial, fiducial_key) sacc_io.save(s, snakemake.output["sacc"], derived_from=[part]) diff --git a/workflow/scripts/cv_plot_cosebis.py b/workflow/scripts/cv_plot_cosebis.py new file mode 100644 index 00000000..7e3f10da --- /dev/null +++ b/workflow/scripts/cv_plot_cosebis.py @@ -0,0 +1,21 @@ +"""Render every COSEBIs figure from the already-sealed multi-cut NPZ.""" + +from sp_validation.b_modes import ( + find_conservative_scale_cut_key, + load_cosebis_results, + plot_cosebis_covariance_matrix, + plot_cosebis_modes, + plot_cosebis_scale_cut_heatmap, +) + +p = snakemake.params +results, edges = load_cosebis_results(snakemake.input.npz) +cut = tuple(p.fiducial_scale_cut) +fiducial = results[find_conservative_scale_cut_key(results, cut)] +plot_cosebis_modes(fiducial, p.version, snakemake.output.figure_modes, cut) +plot_cosebis_covariance_matrix( + fiducial, p.version, "analytic", snakemake.output.figure_covariance +) +plot_cosebis_scale_cut_heatmap( + results, edges, p.version, snakemake.output.figure_scalecut_ptes, cut +) diff --git a/workflow/scripts/cv_plot_pure_eb.py b/workflow/scripts/cv_plot_pure_eb.py new file mode 100644 index 00000000..49699699 --- /dev/null +++ b/workflow/scripts/cv_plot_pure_eb.py @@ -0,0 +1,26 @@ +"""Render pure E/B figures from the already-sealed NPZ only.""" + +from sp_validation.b_modes import ( + covariance_label, + load_pure_eb_results, + plot_eb_covariance_matrix, + plot_integration_vs_reporting, + plot_pte_2d_heatmaps, + plot_pure_eb_correlations, +) + +p = snakemake.params +results = load_pure_eb_results(snakemake.input.npz) +kwargs = dict( + fiducial_xip_scale_cut=tuple(p.fiducial_scale_cut), + fiducial_xim_scale_cut=tuple(p.fiducial_scale_cut), +) +plot_integration_vs_reporting( + results, snakemake.output.figure_integration_vs_reporting, p.version +) +plot_pure_eb_correlations(results, snakemake.output.figure_xis, p.version, **kwargs) +plot_pte_2d_heatmaps(results, p.version, snakemake.output.figure_ptes, **kwargs) +plot_eb_covariance_matrix( + results["cov"], covariance_label(results["npatch"]), + snakemake.output.figure_covariance, p.version, +) diff --git a/workflow/scripts/cv_pure_eb.py b/workflow/scripts/cv_pure_eb.py index 25ee3f72..cf561e08 100644 --- a/workflow/scripts/cv_pure_eb.py +++ b/workflow/scripts/cv_pure_eb.py @@ -15,11 +15,6 @@ from sp_validation.b_modes import ( calculate_eb_statistics, calculate_pure_eb_correlation, - covariance_label, - plot_eb_covariance_matrix, - plot_integration_vs_reporting, - plot_pte_2d_heatmaps, - plot_pure_eb_correlations, save_pure_eb_results, ) from sp_validation.cosmo_val.sacc_writers import pure_eb_to_sacc @@ -46,30 +41,6 @@ ) results = calculate_eb_statistics(results) -plot_integration_vs_reporting( - results, snakemake.output["figure_integration_vs_reporting"], version -) -plot_pure_eb_correlations( - results, - snakemake.output["figure_xis"], - version, - fiducial_xip_scale_cut=fiducial_scale_cut, - fiducial_xim_scale_cut=fiducial_scale_cut, -) -plot_pte_2d_heatmaps( - results, - version, - snakemake.output["figure_ptes"], - fiducial_xip_scale_cut=fiducial_scale_cut, - fiducial_xim_scale_cut=fiducial_scale_cut, -) -plot_eb_covariance_matrix( - results["cov"], - covariance_label(results["npatch"]), - snakemake.output["figure_covariance"], - version, -) - save_pure_eb_results(results, snakemake.output["npz"]) sacc_io.save( From f13b53b2994acb053a719727adb2c4ea6d67103a Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 24 Oct 2026 12:45:22 -0400 Subject: [PATCH 02/18] Cache catalogue and PSF diagnostics; render the full suite without compute --- .../cosmo_val/catalog_characterization.py | 82 +++--- .../cosmo_val/diagnostic_plots.py | 252 ++++++++++++++++++ .../cosmo_val/diagnostic_products.py | 45 ++++ .../cosmo_val/psf_systematics.py | 18 +- workflow/rules/cosmo_val.smk | 198 ++++++++++---- workflow/rules/covariance.smk | 14 +- workflow/rules/twopoint.smk | 3 + workflow/scripts/covariance_plot.py | 7 + workflow/scripts/cv_footprints.py | 17 +- workflow/scripts/cv_objectwise_leakage.py | 27 +- workflow/scripts/cv_plot_2pcf.py | 15 -- workflow/scripts/cv_plot_diagnostic.py | 66 +++++ workflow/scripts/cv_plot_pure_eb.py | 6 +- workflow/scripts/cv_plot_rho_stats.py | 13 - workflow/scripts/cv_plot_tau_stats.py | 18 -- workflow/scripts/cv_ratio_xi_sys_xi.py | 18 -- workflow/scripts/cv_rho_tau_fits.py | 38 +-- workflow/scripts/cv_runner.py | 35 +-- workflow/scripts/cv_weights.py | 16 +- workflow/scripts/run_rho_tau.py | 1 + 20 files changed, 667 insertions(+), 222 deletions(-) create mode 100644 src/sp_validation/cosmo_val/diagnostic_plots.py create mode 100644 src/sp_validation/cosmo_val/diagnostic_products.py create mode 100644 workflow/scripts/covariance_plot.py delete mode 100644 workflow/scripts/cv_plot_2pcf.py create mode 100644 workflow/scripts/cv_plot_diagnostic.py delete mode 100644 workflow/scripts/cv_plot_rho_stats.py delete mode 100644 workflow/scripts/cv_plot_tau_stats.py delete mode 100644 workflow/scripts/cv_ratio_xi_sys_xi.py diff --git a/src/sp_validation/cosmo_val/catalog_characterization.py b/src/sp_validation/cosmo_val/catalog_characterization.py index e2df6d8e..e6983f87 100644 --- a/src/sp_validation/cosmo_val/catalog_characterization.py +++ b/src/sp_validation/cosmo_val/catalog_characterization.py @@ -280,25 +280,32 @@ def calculate_ellipticity_dispersion(self, tomography=False): f"Ellipticity dispersion for tomo_bin_{tomo_bin_id} = {self._ellipticity_dispersion[ver][f'tomo_bin_{tomo_bin_id}']:.4f}" ) - def plot_footprints(self): - self.print_start("Plotting footprints:") - for ver in self.versions: - self.print_magenta(ver) - - fp = cs_plots.FootprintPlotter() + def footprint_maps(self): + """Build one sparse count map per version, shared by every region.""" + from .diagnostic_products import footprint_counts - for region in fp._regions: - out_path = self._output_path(f"footprint_{ver}_{region}.png") - if os.path.exists(out_path): - self.print_done(f"Skipping footprint plot, {out_path} exists") - else: - with self.results[ver].temporarily_read_data(): - hsp_map = fp.create_hsp_map( + for ver in self.versions: + with self.results[ver].temporarily_read_data(): + yield ( + ver, + footprint_counts( self.results[ver].dat_shear["RA"], self.results[ver].dat_shear["Dec"], - ) - fp.plot_region(hsp_map, fp._regions[region], outpath=out_path) - self.print_done("Footprint plot saved to " + out_path) + ), + ) + + def plot_footprints(self): + """Notebook convenience wrapper; workflow plots use persisted maps.""" + from .diagnostic_plots import plot_footprint + + for ver, count_map in self.footprint_maps(): + plot_footprint( + count_map, + { + region: self._output_path(f"footprint_{ver}_{region}.png") + for region in cs_plots.FootprintPlotter._regions + }, + ) def plot_ellipticity(self, nbins=200): out_path = self._output_path("ell_hist.png") @@ -349,37 +356,22 @@ def plot_ellipticity(self, nbins=200): cs_plots.show() self.print_done("Ellipticity histograms saved to " + out_path) - def plot_weights(self, nbins=200): - out_path = self._output_path("weight_hist.png") - if os.path.exists(out_path): - self.print_done(f"Skipping weight histograms, {out_path} exists") - else: - self.print_start("Computing weight histograms:") - - plt.figure(figsize=(10, 7)) - for ver in self.versions: - self.print_magenta(ver) - with self.results[ver].temporarily_read_data(): - w = self._read_shear_cols(ver, "w_col") + def weight_histograms(self, nbins=200): + """Return weight-weighted counts and edges, with independent bins per version.""" + histograms = {} + for ver in self.versions: + with self.results[ver].temporarily_read_data(): + w = self._read_shear_cols(ver, "w_col") + histograms[ver] = np.histogram(w, bins=nbins, weights=w) + return histograms - plt.hist( - w, - bins=nbins, - density=False, - histtype="step", - weights=w, - label=ver, - color=self.cc[ver]["colour"], - ) + def plot_weights(self, nbins=200): + """Notebook convenience wrapper; workflow plots read cached histograms.""" + from .diagnostic_plots import plot_weights - plt.xlabel("$w$") - plt.ylabel("frequency") - plt.yscale("log") - plt.legend() - # plt.xlim([-0.01, 1.2]) - cs_plots.savefig(out_path, close_fig=False) - cs_plots.show() - self.print_done("Weight histograms saved to " + out_path) + plot_weights( + self.weight_histograms(nbins), self.cc, self._output_path("weight_hist.png") + ) def plot_separation(self, nbins=200): self.print_start("Separation histograms") diff --git a/src/sp_validation/cosmo_val/diagnostic_plots.py b/src/sp_validation/cosmo_val/diagnostic_plots.py new file mode 100644 index 00000000..59ba45ef --- /dev/null +++ b/src/sp_validation/cosmo_val/diagnostic_plots.py @@ -0,0 +1,252 @@ +"""Catalogue-free renderers for cached validation diagnostics. + +Inputs are arrays, summary dictionaries or a sparse count map, never a +CosmologyValidation instance. No renderer skips an existing output. +""" + +import matplotlib.pyplot as plt +import numpy as np +from cs_util.plots import FootprintPlotter + + +def plot_footprint(count_map, outputs): + """Render every requested region from one cached HealSparse count map.""" + plotter = FootprintPlotter() + for region, path in outputs.items(): + plotter.plot_region(count_map, plotter._regions[region], outpath=str(path)) + plt.close() + + +def plot_weights(histograms, styles, output): + """Draw weight-weighted bin sums; binning is fixed in the compute step.""" + fig, ax = plt.subplots(figsize=(10, 7)) + for ver, (counts, edges) in histograms.items(): + ax.stairs(counts, edges, label=ver, color=styles[ver]["colour"]) + ax.set(xlabel="$w$", ylabel="frequency", yscale="log") + ax.legend() + fig.savefig(output, bbox_inches="tight") + plt.close(fig) + + +def plot_objectwise_leakage(coefficients, styles, output): + """Draw cached regression and scale-dependent alpha summaries. + + Each coefficient is a (nominal value, standard deviation) pair. + """ + fig, ax = plt.subplots(figsize=(15, 15)) + summaries = { + "alpha_mean": (r"$\bar\alpha$", "-", "full"), + "alpha_1": (r"$\alpha(\theta_{\rm min})$", "--", "none"), + "alpha_0": (r"$\alpha(0)$", ":", "left"), + } + handles = [] + for ver, bins in coefficients.items(): + style = styles[ver] + for coeff in bins.values(): + x, dx = coeff["aii_mean"] + for key, (_, ls, fillstyle) in summaries.items(): + y, dy = coeff[key] + eb = ax.errorbar( + x, + y, + xerr=dx, + yerr=dy, + fmt=style["marker"], + color=style["colour"], + fillstyle=fillstyle, + ) + eb[-1][0].set_linestyle(ls) + handles.append( + plt.Line2D( + [], + [], + marker=style["marker"], + ls="", + color=style["colour"], + label=style.get("label", ver), + ) + ) + handles += [ + plt.Line2D( + [], + [], + marker="o", + ls=ls, + color="grey", + fillstyle=fs, + label=label, + ) + for label, ls, fs in summaries.values() + ] + ax.plot([-0.02, 0.02], [-0.02, 0.02], "k:", linewidth=0.5) + ax.legend(handles=handles) + ax.set( + xlabel=r"$\langle a_{ii} \rangle$ (object-wise)", + ylabel=r"$\alpha$ (scale-dependent, $\tau_0 / \rho_0$)", + ) + fig.savefig(output, bbox_inches="tight") + plt.close(fig) + + +def plot_rho_tau_fits(products, styles, contours, samples, summary): + """Render posterior contours and propagated PSF errors from cached arrays.""" + from shear_psf_leakage import plots as psf_plots + + versions = list(products) + psf_plots.plot_contours( + [products[v]["flat_samples"] for v in versions], + names=["x0", "x1", "x2"], + labels=[r"\alpha", r"\beta", r"\eta"], + savefig=str(contours), + legend_labels=versions, + legend_loc="upper right", + contour_colors=[styles[v]["colour"] for v in versions], + markers={"x0": 0, "x1": 1, "x2": 1}, + show=False, + close=True, + ) + for path, draw_samples in ((samples, True), (summary, False)): + fig, axes = plt.subplots(1, 2, figsize=(14, 6)) + for ver, data in products.items(): + theta, colour = data["theta"], styles[ver]["colour"] + for ax, component, sign in zip(axes, ("plus", "minus"), ("+", "-")): + if draw_samples: + ax.plot( + theta, + data[f"samples_{component}"][-100:].T, + color=colour, + alpha=0.3, + ) + else: + ax.plot(theta, data[f"mean_{component}"], color=colour, label=ver) + lo, hi = data[f"quantiles_{component}"] + ax.plot(theta, lo, color=colour) + ax.plot(theta, hi, color=colour) + ax.fill_between(theta, lo, hi, color=colour, alpha=0.3) + ax.set( + xscale="log", + yscale="log", + xlabel=r"$\theta$ [arcmin]", + ylabel=rf"$\xi^{{\rm PSF,sys}}_{sign}(\theta)$", + ) + if not draw_samples: + axes[0].legend() + fig.savefig(path, bbox_inches="tight") + plt.close(fig) + + +def plot_ratio_xi_sys_xi( + products, correlations, styles, output, offset=0.1, threshold=0.1 +): + """Draw PSF/signal ratios directly from sealed reporting correlations. + + Ratios aren't persisted separately: the only cosmological input stays the + existing sealed SACC. Propagating cached variances is just presentation. + """ + fig, axes = plt.subplots(1, 2, figsize=(14, 6)) + for idx, (ver, data) in enumerate(products.items()): + gg = correlations[ver] + theta = gg.meanr * (1 + offset * (idx - (len(products) - 1) / 2)) + for ax, xi, var_xi, component, sign in zip( + axes, + (gg.xip, gg.xim), + (gg.varxip, gg.varxim), + ("plus", "minus"), + ("+", "-"), + ): + mean, var = data[f"mean_{component}"], data[f"var_{component}"] + ax.errorbar( + theta, + mean / xi, + yerr=np.sqrt(var / xi**2 + mean**2 * var_xi / xi**4), + fmt=styles[ver]["marker"], + color=styles[ver]["colour"], + markersize=3, + capsize=2, + label=ver, + ) + ax.set( + xscale="log", + xlabel=r"$\theta$ [arcmin]", + ylabel=rf"$\xi^{{\rm PSF,sys}}_{sign}/\xi_{sign}$", + ) + for ax in axes: + ax.axhspan(-threshold, threshold, color="black", alpha=0.1) + ax.legend() + fig.savefig(output, bbox_inches="tight") + plt.close(fig) + + +def plot_2pcf(correlations, styles, outputs): + """Draw sealed reporting xi±, both plain and theta-weighted.""" + for times_theta, output in enumerate(outputs): + fig, axes = plt.subplots(1, 2, figsize=(14, 6)) + for ver, gg in correlations.items(): + factor = gg.meanr if times_theta else 1 + for ax, xi, var, sign in zip( + axes, (gg.xip, gg.xim), (gg.varxip, gg.varxim), ("+", "-") + ): + ax.errorbar( + gg.meanr, + xi * factor, + yerr=np.sqrt(var) * factor, + fmt=styles[ver]["marker"], + color=styles[ver]["colour"], + label=ver, + ) + ax.set( + xscale="log", + xlabel=r"$\theta$ [arcmin]", + ylabel=(r"$\theta$ " if times_theta else "") + + rf"$\xi_{sign}(\theta)$", + ) + axes[0].legend() + fig.savefig(output, bbox_inches="tight") + plt.close(fig) + + +def plot_psf_stats(tables, styles, output, kind, covariances=()): + """Overlay cached rho/tau FITS columns, with optional tau covariance errors.""" + indices = (0, 1, 2, 3, 4, 5) if kind == "rho" else (0, 2, 5) + fig, axes = plt.subplots( + 2 if kind == "rho" else 1, + 3, + figsize=(15, 10 if kind == "rho" else 6), + squeeze=False, + ) + for i, stat in enumerate(indices): + ax = axes.flat[i] + for j, (ver, table) in enumerate(tables.items()): + theta = table["theta"] + key = f"{kind}_{stat}_p" + variance = table[f"var{key}"] + if covariances: + n = len(theta) + variance = np.diag(covariances[j])[i * n : (i + 1) * n] + ax.errorbar( + theta, + table[key], + yerr=np.sqrt(variance), + color=styles[ver]["colour"], + label=ver, + capsize=2, + ) + ax.set( + xscale="log", + xlabel=r"$\theta$ [arcmin]", + ylabel=rf"$\{kind}_{{{stat},+}}(\theta)$", + ) + ax.legend() + fig.savefig(output, bbox_inches="tight") + plt.close(fig) + + +def plot_covariance(matrix, output): + """Draw a correlation matrix from the cached analytic covariance.""" + from sp_validation.statistics import corr_from_cov + + fig, ax = plt.subplots(figsize=(9, 9)) + image = ax.imshow(corr_from_cov(matrix), vmin=-1, vmax=1, cmap="RdBu_r") + fig.colorbar(image, ax=ax, label="correlation") + fig.savefig(output, bbox_inches="tight") + plt.close(fig) diff --git a/src/sp_validation/cosmo_val/diagnostic_products.py b/src/sp_validation/cosmo_val/diagnostic_products.py new file mode 100644 index 00000000..a70fb772 --- /dev/null +++ b/src/sp_validation/cosmo_val/diagnostic_products.py @@ -0,0 +1,45 @@ +"""Small non-cosmological diagnostic products shared by compute and plot steps.""" + +import healpy as hp +import healsparse +import numpy as np + + +def footprint_counts(ra, dec, nside_coverage=32, nside_sparse=2048): + """Bin coordinates once, using a finite HealSparse sentinel for valid pixels.""" + pixels, counts = np.unique( + hp.ang2pix(nside_sparse, ra, dec, lonlat=True, nest=True), return_counts=True + ) + count_map = healsparse.HealSparseMap.make_empty( + nside_coverage, nside_sparse, dtype=np.float32 + ) + count_map[pixels] = counts.astype(np.float32) + return count_map + + +def save_footprint(count_map, path): + """Store occupied nested pixels, counts and both HealSparse resolutions. + + A sparse NPZ avoids writing the plotter's NaN sentinel to a FITS header + (CFITSIO cannot encode it). No per-object coordinates are retained. + """ + pixels = count_map.valid_pixels + np.savez_compressed( + path, + pixels=pixels, + counts=count_map[pixels], + nside_coverage=count_map.nside_coverage, + nside_sparse=count_map.nside_sparse, + ) + + +def load_footprint(path): + """Reconstruct the sparse map from occupied pixels, without a catalogue.""" + with np.load(path, allow_pickle=False) as data: + count_map = healsparse.HealSparseMap.make_empty( + int(data["nside_coverage"]), + int(data["nside_sparse"]), + dtype=np.float32, + ) + count_map[data["pixels"]] = data["counts"] + return count_map diff --git a/src/sp_validation/cosmo_val/psf_systematics.py b/src/sp_validation/cosmo_val/psf_systematics.py index 991567fb..ef0e004b 100644 --- a/src/sp_validation/cosmo_val/psf_systematics.py +++ b/src/sp_validation/cosmo_val/psf_systematics.py @@ -303,7 +303,23 @@ def _objectwise_leakage_bin(self, results_obj, ver, tomo_bin_id, mix, order): else: self.print_cyan(f"Computing object-wise leakage regression, {suffix}") with results_obj.temporarily_read_data(selection=selection): - results_obj.PSF_leakage(mix=mix, order=order, suffix=suffix) + # Use the regression primitive, not PSF_leakage: that driver + # also fits unrelated 1D models and writes exploratory plots. + e, weights = results_obj.get_ellipticity_weights() + x = [ + results_obj._dat[results_obj._params[key]] + for key in ("e1_PSF_col", "e2_PSF_col") + ] + results_obj.par_best_fit = leakage.corr_2d( + x, + e, + weights=weights, + order=order, + mix=mix, + stats_file=results_obj._stats_file, + verbose=results_obj._params["verbose"], + ) + leakage.save_to_file(results_obj.par_best_fit, out_path) par_best_fit = results_obj.par_best_fit a11 = ufloat(par_best_fit["a11"].value, par_best_fit["a11"].stderr) diff --git a/workflow/rules/cosmo_val.smk b/workflow/rules/cosmo_val.smk index c6078f14..eb37511e 100644 --- a/workflow/rules/cosmo_val.smk +++ b/workflow/rules/cosmo_val.smk @@ -22,13 +22,11 @@ # NaMaster pseudo-Cℓ) are what assembly puts in the terminal file, replacing # the estimates a part was born with. # -# Methods that only emit figures, or whose figure paths derive from internal -# handler state (rho/tau plots, rho_tau_fits, objectwise leakage, 2pcf -# overlay), declare a sentinel under COSMO_VAL/snakemake_sentinels so they stay -# DAG-trackable. Lazy cv state the original code never persists (c1/c2, -# xi_psf_sys) is either materialized to a small JSON (additive bias) or -# recomputed in the one rule that needs it — recompute is cheap next to the -# science it depends on. See workflow/scripts/cv_runner.py. +# Every compute rule writes a data product. Every figure is declared by a +# separate rule consuming those products, with no catalogue inputs or lazy cv +# state. cosmo_val_plots gathers all figures; cosmo_val_all also assembles the +# analysis products. Cosmological caches contain only already-sealed vectors. +# Plot modules are explicit inputs so edits to them invalidate figures only. CV = config["cosmo_val"] CV_VERSIONS = config["versions"] @@ -183,6 +181,31 @@ def cv_analysis_sacc(version): # The CosmologyValidation constructor kwargs every cv_runner rule passes. CV_INIT = cv_init_params(config) +CV_STYLES = { + v: {k: catalogue_entry(v).get(k, default) for k, default in + (("colour", "C0"), ("marker", "o"), ("label", v))} + for v in CV_VERSIONS +} +CV_PLOT_CODE = str(REPO_ROOT / "src/sp_validation/cosmo_val/diagnostic_plots.py") +CV_REGIONS = ("NGC", "SGC", "fullsky") +CV_PSF_PRODUCTS = [str(COSMO_VAL / f"psf_fit_{v}_{CV_BINNING}.npz") for v in CV_VERSIONS] +CV_FOOTPRINT_MAPS = [str(COSMO_VAL / f"footprint_{v}.npz") for v in CV_VERSIONS] +CV_FOOTPRINT_FIGURES = [ + str(COSMO_VAL / f"footprint_{v}_{r}.png") for v in CV_VERSIONS for r in CV_REGIONS +] + + +def cv_shear_catalogues(): + """Catalogue inputs belong to compute rules only.""" + return [str(Path(catalogue_entry(v)["subdir"]) / catalogue_entry(v)["shear"]["path"]) + for v in CV_VERSIONS] + + +def cv_tau_covariances(): + method = CV.get("cov_estimate_method", "th") + return [str(COSMO_VAL / "rho_tau_stats" / + f"cov_tau_{cv_basename(v, CV_FIDUCIAL)}_{method}.npy") + for v in CV_VERSIONS] # --------------------------------------------------------------------------- @@ -198,29 +221,35 @@ CV_INIT = cv_init_params(config) rule cv_plot_rho_stats: """Overlay rho statistics across all versions.""" input: - rho=[cv_rho_stats(v, CV_FIDUCIAL) for v in CV_VERSIONS], + stats=[cv_rho_stats(v, CV_FIDUCIAL) for v in CV_VERSIONS], + covariances=[], + plot_code=CV_PLOT_CODE, output: - sentinel=str(CV_SENTINELS / "plot_rho_stats.done"), + figure=str(COSMO_VAL / "rho_tau_stats/rho_stats.png"), params: - cv_init=CV_INIT, + kind="rho", + styles=CV_STYLES, resources: runtime=20, script: - "../scripts/cv_plot_rho_stats.py" + "../scripts/cv_plot_diagnostic.py" rule cv_plot_tau_stats: """Overlay tau statistics across all versions.""" input: - tau=[cv_tau_stats(v, CV_FIDUCIAL) for v in CV_VERSIONS], + stats=[cv_tau_stats(v, CV_FIDUCIAL) for v in CV_VERSIONS], + covariances=cv_tau_covariances(), + plot_code=CV_PLOT_CODE, output: - sentinel=str(CV_SENTINELS / "plot_tau_stats.done"), + figure=str(COSMO_VAL / "rho_tau_stats/tau_stats.png"), params: - cv_init=CV_INIT, + kind="tau", + styles=CV_STYLES, resources: runtime=20, script: - "../scripts/cv_plot_tau_stats.py" + "../scripts/cv_plot_diagnostic.py" rule cv_rho_tau_fits: @@ -228,8 +257,11 @@ rule cv_rho_tau_fits: input: rho=[cv_rho_stats(v, CV_FIDUCIAL) for v in CV_VERSIONS], tau=[cv_tau_stats(v, CV_FIDUCIAL) for v in CV_VERSIONS], + covariances=cv_tau_covariances(), + rho_covariances=[str(COSMO_VAL / "rho_tau_stats" / + f"cov_rho_{cv_basename(v, CV_FIDUCIAL)}_jk.npy") for v in CV_VERSIONS], output: - sentinel=str(CV_SENTINELS / "rho_tau_fits.done"), + products=CV_PSF_PRODUCTS, params: cv_init=CV_INIT, resources: @@ -244,9 +276,11 @@ rule cv_rho_tau_fits: # --------------------------------------------------------------------------- rule cv_footprints: - """Per-version survey footprint maps.""" + """Cache per-version sparse survey count maps.""" + input: + catalogues=cv_shear_catalogues(), output: - sentinel=str(CV_SENTINELS / "footprints.done"), + maps=CV_FOOTPRINT_MAPS, params: cv_init=CV_INIT, resources: @@ -257,12 +291,14 @@ rule cv_footprints: rule cv_objectwise_leakage: - """Object-wise PSF-leakage regression vs scale-dependent alpha (all versions).""" + """Cache object-wise regression and scale-dependent alpha summaries.""" input: + catalogues=cv_shear_catalogues(), + covariances=cv_tau_covariances(), rho=[cv_rho_stats(v, CV_FIDUCIAL) for v in CV_VERSIONS], tau=[cv_tau_stats(v, CV_FIDUCIAL) for v in CV_VERSIONS], output: - sentinel=str(CV_SENTINELS / "objectwise_leakage.done"), + coefficients=str(COSMO_VAL / "leakage_coefficients.json"), params: cv_init=CV_INIT, threads: 12 @@ -274,9 +310,11 @@ rule cv_objectwise_leakage: rule cv_weights: - """Weighted shear-weight histograms across versions.""" + """Cache weighted shear-weight histogram counts and edges.""" + input: + catalogues=cv_shear_catalogues(), output: - weight_hist=str(COSMO_VAL / "weight_hist.png"), + histograms=str(COSMO_VAL / "weight_hist.npz"), params: cv_init=CV_INIT, resources: @@ -295,6 +333,61 @@ rule cv_weights: # standalone diagnostic (the original driver's cv.calculate_additive_bias() # call); the xi rule does not depend on it. +rule cv_plot_footprints: + input: + maps=CV_FOOTPRINT_MAPS, + plot_code=CV_PLOT_CODE, + output: + figures=CV_FOOTPRINT_FIGURES, + params: + kind="footprints", + styles=CV_STYLES, + regions=CV_REGIONS, + script: + "../scripts/cv_plot_diagnostic.py" + + +rule cv_plot_weights: + input: + histograms=rules.cv_weights.output.histograms, + plot_code=CV_PLOT_CODE, + output: + figure=str(COSMO_VAL / "weight_hist.png"), + params: + kind="weights", + styles=CV_STYLES, + script: + "../scripts/cv_plot_diagnostic.py" + + +rule cv_plot_objectwise_leakage: + input: + coefficients=rules.cv_objectwise_leakage.output.coefficients, + plot_code=CV_PLOT_CODE, + output: + figure=str(COSMO_VAL / "leakage_coefficients.png"), + params: + kind="objectwise_leakage", + styles=CV_STYLES, + script: + "../scripts/cv_plot_diagnostic.py" + + +rule cv_plot_rho_tau_fits: + input: + products=CV_PSF_PRODUCTS, + plot_code=CV_PLOT_CODE, + output: + contours=str(COSMO_VAL / "rho_tau_stats/contours_tau_stat.png"), + samples=str(COSMO_VAL / "rho_tau_stats/xi_psf_sys_tomography_False_sample.png"), + summary=str(COSMO_VAL / "rho_tau_stats/xi_psf_sys_tomography_False_mean_and_std.png"), + params: + kind="rho_tau_fits", + styles=CV_STYLES, + script: + "../scripts/cv_plot_diagnostic.py" + + rule cv_additive_bias: """Weighted mean ellipticity c1/c2 per version (standalone diagnostic).""" output: @@ -309,35 +402,40 @@ rule cv_additive_bias: rule cv_plot_2pcf: - """n_pairs / xi± overlay across versions, from the reporting parts.""" + """Xi± overlays across versions, from the sealed reporting parts.""" input: xi=[cv_xi_sacc(v, "reporting") for v in CV_VERSIONS], + plot_code=CV_PLOT_CODE, output: - sentinel=str(CV_SENTINELS / "plot_2pcf.done"), + figures=[str(COSMO_VAL / name) for name in + ("xi_pm_tomography_False.png", "xi_pm_theta_tomography_False.png")], params: - cv_init=CV_INIT, + kind="xi", + styles=CV_STYLES, + blinds=[blind_of(v) for v in CV_VERSIONS], resources: runtime=20, script: - "../scripts/cv_plot_2pcf.py" + "../scripts/cv_plot_diagnostic.py" rule cv_ratio_xi_sys_xi: """Ratio of PSF systematics (xi_psf_sys) to the cosmic-shear signal (xi+).""" input: xi=[cv_xi_sacc(v, "reporting") for v in CV_VERSIONS], - rho=[cv_rho_stats(v, CV_FIDUCIAL) for v in CV_VERSIONS], - tau=[cv_tau_stats(v, CV_FIDUCIAL) for v in CV_VERSIONS], + products=CV_PSF_PRODUCTS, + plot_code=CV_PLOT_CODE, output: - ratio=str(COSMO_VAL / "ratio_xi_sys_xi.png"), + figure=str(COSMO_VAL / "ratio_xi_sys_xi.png"), params: - offset=0.1, - cv_init=CV_INIT, + kind="ratio", + styles=CV_STYLES, + blinds=[blind_of(v) for v in CV_VERSIONS], resources: mem_mb=16000, runtime=120, script: - "../scripts/cv_ratio_xi_sys_xi.py" + "../scripts/cv_plot_diagnostic.py" # --------------------------------------------------------------------------- @@ -357,6 +455,7 @@ rule cv_plot_pseudo_cl: input: pseudo_cl=[cv_pseudo_cl_analysis_sacc(v) for v in CV_VERSIONS], pseudo_cl_cov=[cv_pseudo_cl_cov(v) for v in CV_VERSIONS], + plot_code=str(REPO_ROOT / "src/sp_validation/cosmo_val/pseudo_cl.py"), output: **cv_pseudo_cl_figures(), params: @@ -547,23 +646,28 @@ rule assemble_sacc_all: # Aggregate target: the whole validation suite # --------------------------------------------------------------------------- +rule cosmo_val_plots: + """Render every validation figure from cached data products.""" + input: + rules.cv_plot_rho_stats.output, + rules.cv_plot_tau_stats.output, + rules.cv_plot_rho_tau_fits.output, + rules.cv_plot_footprints.output, + rules.cv_plot_objectwise_leakage.output, + rules.cv_plot_weights.output, + rules.cv_plot_2pcf.output, + rules.cv_ratio_xi_sys_xi.output, + [path for v in CV_VERSIONS for path in cv_pure_eb_figures(v).values()], + [path for v in CV_VERSIONS for path in cv_cosebis_figures(v).values()], + list(cv_pseudo_cl_figures().values()) if CV.get("include_pseudo_cl", False) else [], + [cv_xi_cov(v).removesuffix(".dat") + "_corr_plot.pdf" for v in CV_VERSIONS], + [cv_xi_cov_integration(v).removesuffix(".dat") + "_corr_plot.pdf" for v in CV_VERSIONS], + + rule cosmo_val_all: - """Run the full CosmologyValidation diagnostic suite.""" + """Run the full suite, including every plot-only rule.""" input: - # PSF diagnostics - str(CV_SENTINELS / "footprints.done"), - str(CV_SENTINELS / "plot_rho_stats.done"), - str(CV_SENTINELS / "plot_tau_stats.done"), - str(CV_SENTINELS / "rho_tau_fits.done"), - # Shear diagnostics - str(CV_SENTINELS / "objectwise_leakage.done"), - str(COSMO_VAL / "weight_hist.png"), + rules.cosmo_val_plots.input, str(COSMO_VAL / "additive_bias.json"), - # Two-point - str(CV_SENTINELS / "plot_2pcf.done"), - str(COSMO_VAL / "ratio_xi_sys_xi.png"), - # B-modes str(COSMO_VAL / "bmode_summary.json"), - list(cv_pseudo_cl_figures().values()) if CV.get("include_pseudo_cl", False) else [], - # Terminal analysis file: the assembled {version}.sacc per version [cv_analysis_sacc(v) for v in CV_VERSIONS], diff --git a/workflow/rules/covariance.smk b/workflow/rules/covariance.smk index c3079760..84e9f14a 100644 --- a/workflow/rules/covariance.smk +++ b/workflow/rules/covariance.smk @@ -123,7 +123,6 @@ rule covariance_ini: params: outdir=lambda w: covariance_dir(*_onecov_args(w)), out_filename=lambda w: covariance_base(*_onecov_args(w)) + ".dat", - out_plot_filename=lambda w: covariance_base(*_onecov_args(w)) + "_corr_plot.pdf", ng_value=lambda w: "1" if w.gaussian == "ng" else "0", omega_m=PLANCK18["Omega_m"], omega_v=PLANCK18["Omega_v"], @@ -173,7 +172,7 @@ directory = {params.outdir} file = {params.out_filename} style = matrix save_configs = False -corrmatrix_plot = {params.out_plot_filename} +corrmatrix_plot = False [covTHETAspace settings] theta_min = {wildcards.min_sep} @@ -287,6 +286,17 @@ rule covariance_run: """ +rule covariance_plot: + """Render the correlation matrix from cached OneCovariance output only.""" + input: + matrix=cov_output(".dat"), + plot_code=str(REPO_ROOT / "src/sp_validation/cosmo_val/diagnostic_plots.py"), + output: + figure=cov_output("_corr_plot.pdf"), + script: + "../scripts/covariance_plot.py" + + def fiducial_covariance_outputs(mask_suffix="", tomo_suffix=""): """Return the covariance matrices for the fiducial version: non-Gaussian on the reporting grid, Gaussian on the integration grid.""" diff --git a/workflow/rules/twopoint.smk b/workflow/rules/twopoint.smk index 5d6450e6..13e1a3f6 100644 --- a/workflow/rules/twopoint.smk +++ b/workflow/rules/twopoint.smk @@ -48,11 +48,14 @@ rule rho_tau_stats: output: rho_stats=str(COSMO_VAL / "rho_tau_stats/rho_stats_{version}_tomo_bin_all_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits"), tau_stats=str(COSMO_VAL / "rho_tau_stats/tau_stats_{version}_tomo_bin_all_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.fits"), + covariance=str(COSMO_VAL / ("rho_tau_stats/cov_tau_{version}_tomo_bin_all_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}_" + config.get("cosmo_val", {}).get("cov_estimate_method", "th") + ".npy")), + rho_covariance=str(COSMO_VAL / "rho_tau_stats/cov_rho_{version}_tomo_bin_all_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}_jk.npy"), # Born-as-SACC ρ/τ part, written alongside the FITS. rho_tau=str(COSMO_VAL / "rho_tau_stats/rho_tau_{version}_tomo_bin_all_minsep={min_sep}_maxsep={max_sep}_nbins={nbins}_npatch={npatch}.sacc"), threads: 48 params: ver="{version}", + cov_estimate_method=config.get("cosmo_val", {}).get("cov_estimate_method", "th"), min_sep="{min_sep}", max_sep="{max_sep}", nbins="{nbins}", diff --git a/workflow/scripts/covariance_plot.py b/workflow/scripts/covariance_plot.py new file mode 100644 index 00000000..2af67d05 --- /dev/null +++ b/workflow/scripts/covariance_plot.py @@ -0,0 +1,7 @@ +"""Plot an analytic covariance without invoking its producer.""" + +import numpy as np + +from sp_validation.cosmo_val.diagnostic_plots import plot_covariance + +plot_covariance(np.loadtxt(snakemake.input.matrix), snakemake.output.figure) diff --git a/workflow/scripts/cv_footprints.py b/workflow/scripts/cv_footprints.py index e4a1af6a..90d36015 100644 --- a/workflow/scripts/cv_footprints.py +++ b/workflow/scripts/cv_footprints.py @@ -1,13 +1,12 @@ -"""Rule cv_footprints: plot per-version survey footprints. +"""Build one sparse count map per catalogue; no figures.""" -Leaf rule. Reads the shear catalog (RA/Dec) for each version; writes -footprint_{version}_{region}.png. Sentinel output keeps the DAG trackable -because region names are internal to FootprintPlotter. -""" +from cv_runner import make_cv, verify_outputs -from cv_runner import _unbuffer_streams, make_cv, touch_sentinels +from sp_validation.cosmo_val.diagnostic_products import save_footprint -_unbuffer_streams() cv = make_cv(snakemake) -cv.plot_footprints() -touch_sentinels(snakemake) +for (_, count_map), path in zip( + cv.footprint_maps(), snakemake.output.maps, strict=True +): + save_footprint(count_map, path) +verify_outputs(snakemake) diff --git a/workflow/scripts/cv_objectwise_leakage.py b/workflow/scripts/cv_objectwise_leakage.py index ae0f012a..2c570191 100644 --- a/workflow/scripts/cv_objectwise_leakage.py +++ b/workflow/scripts/cv_objectwise_leakage.py @@ -1,15 +1,20 @@ -"""Rule cv_objectwise_leakage: object-wise vs scale-dependent PSF leakage. +"""Cache regression coefficients and alpha summaries, without rendering.""" -plot_objectwise_leakage triggers calculate_objectwise_leakage, which itself -triggers calculate_scale_dependent_leakage — the whole leakage chain runs here. -Writes leakage_{version}/ products (the object-wise regression .pkl is the -durable artifact) and leakage_coefficients.png. Sentinel-tracked because the -per-version leakage product paths are built inside shear_psf_leakage. -""" +import json -from cv_runner import _unbuffer_streams, make_cv, touch_sentinels +from cv_runner import make_cv, prepare_cached_psf, verify_outputs -_unbuffer_streams() cv = make_cv(snakemake) -cv.plot_objectwise_leakage() -touch_sentinels(snakemake) +prepare_cached_psf(cv) +cv.calculate_objectwise_leakage(tomography=False) +cv.calculate_alpha_leakage_summaries(tomography=False) +coefficients = { + ver: { + bin_id: {k: [v.nominal_value, v.std_dev] for k, v in coeff.items()} + for bin_id, coeff in bins.items() + } + for ver, bins in cv.leakage_coeff.items() +} +with open(snakemake.output.coefficients, "w") as stream: + json.dump(coefficients, stream, indent=2) +verify_outputs(snakemake) diff --git a/workflow/scripts/cv_plot_2pcf.py b/workflow/scripts/cv_plot_2pcf.py deleted file mode 100644 index 5bb3ac32..00000000 --- a/workflow/scripts/cv_plot_2pcf.py +++ /dev/null @@ -1,15 +0,0 @@ -"""Rule cv_plot_2pcf: xi± overlay across versions. - -Draws each version's reporting ξ± part (the declared inputs) for the -non-tomographic ("all", "all") pair; plot_2pcf draws the parts it is handed -rather than measuring. Writes figures under the output dir. Sentinel-tracked: -plot_2pcf emits several figures whose names are internal. -""" - -from cv_runner import _unbuffer_streams, make_cv, take_reporting_parts, touch_sentinels - -_unbuffer_streams() -cv = make_cv(snakemake) -take_reporting_parts(cv, snakemake.input["xi"]) -cv.plot_2pcf(tomography=False, show=False) -touch_sentinels(snakemake) diff --git a/workflow/scripts/cv_plot_diagnostic.py b/workflow/scripts/cv_plot_diagnostic.py new file mode 100644 index 00000000..7c7683b2 --- /dev/null +++ b/workflow/scripts/cv_plot_diagnostic.py @@ -0,0 +1,66 @@ +"""Catalogue-free plotting entry point for cached scalar/PSF/map diagnostics.""" + +import json + +import numpy as np +from astropy.io import fits + +from sp_validation import sacc_io +from sp_validation.cosmo_val import diagnostic_plots as plots +from sp_validation.cosmo_val.diagnostic_products import load_footprint + +p, inputs, outputs = snakemake.params, snakemake.input, snakemake.output +styles = dict(p.styles) +versions = list(styles) + + +def psf_products(): + products = {} + for ver, path in zip(versions, inputs.products, strict=True): + with np.load(path, allow_pickle=False) as data: + products[ver] = dict(data) + return products + + +if p.kind == "weights": + with np.load(inputs.histograms, allow_pickle=False) as data: + histograms = { + ver: (data[f"counts_{i}"], data[f"edges_{i}"]) + for i, ver in enumerate(data["versions"].tolist()) + } + plots.plot_weights(histograms, styles, outputs.figure) +elif p.kind == "footprints": + for i, (ver, path) in enumerate(zip(versions, inputs.maps, strict=True)): + plots.plot_footprint( + load_footprint(path), + { + region: outputs.figures[i * len(p.regions) + j] + for j, region in enumerate(p.regions) + }, + ) +elif p.kind == "objectwise_leakage": + with open(inputs.coefficients) as stream: + coefficients = json.load(stream) + plots.plot_objectwise_leakage(coefficients, styles, outputs.figure) +elif p.kind == "rho_tau_fits": + plots.plot_rho_tau_fits( + psf_products(), styles, outputs.contours, outputs.samples, outputs.summary + ) +elif p.kind in ("ratio", "xi"): + correlations = { + ver: sacc_io.xi_correlation(sacc_io.load(path, blind=blind), grid="reporting") + for ver, path, blind in zip(versions, inputs.xi, p.blinds, strict=True) + } + if p.kind == "ratio": + plots.plot_ratio_xi_sys_xi(psf_products(), correlations, styles, outputs.figure) + else: + plots.plot_2pcf(correlations, styles, outputs.figures) +elif p.kind in ("rho", "tau"): + tables = { + ver: fits.getdata(path) + for ver, path in zip(versions, inputs.stats, strict=True) + } + covariances = [np.load(path) for path in inputs.covariances] + plots.plot_psf_stats(tables, styles, outputs.figure, p.kind, covariances) +else: + raise ValueError(f"Unknown diagnostic: {p.kind}") diff --git a/workflow/scripts/cv_plot_pure_eb.py b/workflow/scripts/cv_plot_pure_eb.py index 49699699..fec9a0e5 100644 --- a/workflow/scripts/cv_plot_pure_eb.py +++ b/workflow/scripts/cv_plot_pure_eb.py @@ -21,6 +21,8 @@ plot_pure_eb_correlations(results, snakemake.output.figure_xis, p.version, **kwargs) plot_pte_2d_heatmaps(results, p.version, snakemake.output.figure_ptes, **kwargs) plot_eb_covariance_matrix( - results["cov"], covariance_label(results["npatch"]), - snakemake.output.figure_covariance, p.version, + results["cov"], + covariance_label(results["npatch"]), + snakemake.output.figure_covariance, + p.version, ) diff --git a/workflow/scripts/cv_plot_rho_stats.py b/workflow/scripts/cv_plot_rho_stats.py deleted file mode 100644 index a55a1377..00000000 --- a/workflow/scripts/cv_plot_rho_stats.py +++ /dev/null @@ -1,13 +0,0 @@ -"""Rule cv_plot_rho_stats: rho-statistics overlay across versions. - -Reads rho_stats_{base}.fits for every version (declared inputs); writes -rho_stats.png into the leakage output dir. Sentinel-tracked because the figure -lands beside the first version's leakage products, not at a fixed path. -""" - -from cv_runner import _unbuffer_streams, make_cv, touch_sentinels - -_unbuffer_streams() -cv = make_cv(snakemake) -cv.plot_rho_stats(savefig="rho_stats.png", show=False) -touch_sentinels(snakemake) diff --git a/workflow/scripts/cv_plot_tau_stats.py b/workflow/scripts/cv_plot_tau_stats.py deleted file mode 100644 index a2c6f365..00000000 --- a/workflow/scripts/cv_plot_tau_stats.py +++ /dev/null @@ -1,18 +0,0 @@ -"""Rule cv_plot_tau_stats: tau-statistics overlay across versions. - -Reads tau_stats_{base}.fits for every version (declared inputs), with error -bars from the version's tau covariance of the configured type; writes -tau_stats.png into the leakage output dir. Sentinel-tracked (see rho_stats). -""" - -from cv_runner import _unbuffer_streams, make_cv, touch_sentinels - -_unbuffer_streams() -cv = make_cv(snakemake) -cv.plot_tau_stats( - tomography=False, - cov_type=cv.cov_estimate_method, - savefig="tau_stats.png", - show=False, -) -touch_sentinels(snakemake) diff --git a/workflow/scripts/cv_ratio_xi_sys_xi.py b/workflow/scripts/cv_ratio_xi_sys_xi.py deleted file mode 100644 index 9e1e4396..00000000 --- a/workflow/scripts/cv_ratio_xi_sys_xi.py +++ /dev/null @@ -1,18 +0,0 @@ -"""Rule cv_ratio_xi_sys_xi: ratio of PSF systematics to cosmic-shear signal. - -Joins two upstream chains for the non-tomographic ("all", "all") pair: each -version's reporting ξ± part and xi_psf_sys (recomputed in memory from the -rho/tau FITS via the lazy cv.xi_psf_sys property — the PSF-error fit is not -persisted). Both are declared as inputs so the DAG shows the join. Writes -ratio_xi_sys_xi.png at a fixed path (declared output). -""" - -from cv_runner import _unbuffer_streams, make_cv, take_reporting_parts, verify_outputs - -_unbuffer_streams() -cv = make_cv(snakemake) -take_reporting_parts(cv, snakemake.input["xi"]) -cv.plot_ratio_xi_sys_xi( - tomography=False, offset=snakemake.params.get("offset", 0.1), show=False -) -verify_outputs(snakemake) diff --git a/workflow/scripts/cv_rho_tau_fits.py b/workflow/scripts/cv_rho_tau_fits.py index 6e02ba10..d2e49378 100644 --- a/workflow/scripts/cv_rho_tau_fits.py +++ b/workflow/scripts/cv_rho_tau_fits.py @@ -1,22 +1,24 @@ -"""Rule cv_rho_tau_fits: PSF-error model fit and xi_psf_sys contours. +"""Fit PSF parameters and persist all non-cosmological plot arrays.""" -Conditional in the original driver on rho_tau_method != "none". Reads -rho_stats/tau_stats FITS (declared inputs), runs the MCMC fit -(calculate_rho_tau_fits via plot_rho_tau_fits), and writes contour / -xi_psf_sys figures into the leakage output dir. Sentinel-tracked: figure paths -derive from internal handler state. The fitted xi_psf_sys it produces is -recomputed (not persisted) by the cv_ratio_xi_sys_xi rule. -""" +import numpy as np +from cv_runner import make_cv, prepare_cached_psf, verify_outputs -from cv_runner import _unbuffer_streams, make_cv, touch_sentinels - -_unbuffer_streams() cv = make_cv(snakemake) -if cv.rho_tau_method != "none": - cv.plot_rho_tau_fits( - tomography=False, - savefig_contours="contours_tau_stat.png", - savefig_xi_psf_sys="xi_psf_sys", - show=False, +prepare_cached_psf(cv) +for ver, path in zip(cv.versions, snakemake.output.products, strict=True): + params = cv.set_params_rho_tau(ver, cv.results[ver]._params, cv.cc[ver]["psf"]) + flat_samples = cv.get_samples(ver, params, "all") + plus, minus = cv.get_xi_psf_sys_samples(ver, params, "all", "all") + data = dict( + theta=cv.psf_fitter.rho_stat_handler.rho_stats["theta"], + flat_samples=flat_samples, ) -touch_sentinels(snakemake) + for component, samples in (("plus", plus), ("minus", minus)): + data[f"samples_{component}"] = samples + data[f"mean_{component}"] = samples.mean(axis=0) + data[f"var_{component}"] = samples.var(axis=0) + data[f"quantiles_{component}"] = np.quantile( + samples, [1 - cv.quantile, cv.quantile], axis=0 + ) + np.savez_compressed(path, **data) +verify_outputs(snakemake) diff --git a/workflow/scripts/cv_runner.py b/workflow/scripts/cv_runner.py index b4ffcea2..14430166 100644 --- a/workflow/scripts/cv_runner.py +++ b/workflow/scripts/cv_runner.py @@ -6,12 +6,10 @@ owns, and (3) verifies the declared outputs landed. This module carries the common machinery so each rule script stays a few readable lines. -The in-memory ``cv`` object is re-instantiated per rule. The DAG links rules -through the *real* data products each method writes (rho/tau FITS, ξ± parts, -pseudo-Cl FITS, E/B npz, leakage pkl); the lazy ``cv`` properties that aren't -persisted (``c1``/``c2``, ``xi_psf_sys``, ``rho_tau_fits``) are recomputed in -whichever rule needs them, which is cheap next to the science compute they -depend on. See common.cv_init_params for the constructor contract. +Compute rules instantiate ``cv`` independently and link through explicit +cached data products. Plot rules read those products directly; they don't +instantiate ``cv`` or enter its lazy catalogue/measurement paths. +See common.cv_init_params for the compute constructor contract. """ import os @@ -36,6 +34,22 @@ def make_cv(snakemake): return CosmologyValidation(**dict(snakemake.params["cv_init"])) +def prepare_cached_psf(cv): + """Attach FITS-only handlers; never enter the catalogue/pair-count path.""" + from shear_psf_leakage.rho_tau_stat import RhoStat, TauStat + + output = str(Path(cv.cc["paths"]["output"]) / "rho_tau_stats") + cv._rho_stat_handler = RhoStat( + output=output, treecorr_config=cv.treecorr_config, verbose=True + ) + cv._tau_stat_handler = TauStat( + catalogs=cv._rho_stat_handler.catalogs, + output=output, + treecorr_config=cv.treecorr_config, + verbose=True, + ) + + def take_reporting_parts(cv, paths): """Hand ``cv`` the reporting ξ± parts at ``paths``, one per version, in order.""" from sp_validation import sacc_io @@ -53,12 +67,3 @@ def verify_outputs(snakemake): "cosmo_val rule finished but these declared outputs are absent:\n " + "\n ".join(missing) ) - - -def touch_sentinels(snakemake): - """Create any sentinel outputs (pure-plot rules with no data product).""" - for out in snakemake.output: - path = Path(out) - if path.suffix == ".done": - path.parent.mkdir(parents=True, exist_ok=True) - path.touch() diff --git a/workflow/scripts/cv_weights.py b/workflow/scripts/cv_weights.py index 0316fe1e..293e91d6 100644 --- a/workflow/scripts/cv_weights.py +++ b/workflow/scripts/cv_weights.py @@ -1,12 +1,12 @@ -"""Rule cv_weights: weighted shear-weight histograms across versions. +"""Cache weight-weighted bin sums and bin edges; no figures.""" -Leaf rule. Reads each version's shear weights; writes weight_hist.png at a -fixed path under the output dir (declared output). -""" +import numpy as np +from cv_runner import make_cv, verify_outputs -from cv_runner import _unbuffer_streams, make_cv, verify_outputs - -_unbuffer_streams() cv = make_cv(snakemake) -cv.plot_weights() +data = {"versions": np.asarray(cv.versions)} +for i, (counts, edges) in enumerate(cv.weight_histograms().values()): + data[f"counts_{i}"] = counts + data[f"edges_{i}"] = edges +np.savez(snakemake.output.histograms, **data) verify_outputs(snakemake) diff --git a/workflow/scripts/run_rho_tau.py b/workflow/scripts/run_rho_tau.py index d04da3f0..21712470 100644 --- a/workflow/scripts/run_rho_tau.py +++ b/workflow/scripts/run_rho_tau.py @@ -18,6 +18,7 @@ npatch=int(params["npatch"]), catalog_config=params["cat_config"], output_dir=params["output_dir"], + cov_estimate_method=params["cov_estimate_method"], ) cv.calculate_rho_tau_stats(tomography=False) verify_outputs(snakemake) From ba5fb8e80e9c3087b439b554289d76d10b242a86 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 24 Oct 2026 12:50:06 -0400 Subject: [PATCH 03/18] Test cached renderers and fail-closed plot-only workflow; document restyling --- src/sp_validation/b_modes.py | 33 +-- src/sp_validation/tests/test_plot_products.py | 278 ++++++++++++++++++ src/sp_validation/tests/test_plot_workflow.py | 145 +++++++++ workflow/README.md | 33 +++ workflow/scripts/cv_rho_tau_fits.py | 6 + 5 files changed, 472 insertions(+), 23 deletions(-) create mode 100644 src/sp_validation/tests/test_plot_products.py create mode 100644 src/sp_validation/tests/test_plot_workflow.py diff --git a/src/sp_validation/b_modes.py b/src/sp_validation/b_modes.py index ea9a80ca..d4cc548f 100644 --- a/src/sp_validation/b_modes.py +++ b/src/sp_validation/b_modes.py @@ -921,28 +921,14 @@ def plot_cosebis_scale_cut_heatmap( left_edges, right_edges = edges nbins = len(left_edges) - # Initialize matrices - snrs = [np.sqrt(result["chi2_E"]) for result in cosebis_results.values()] - ptes = [result["pte_B"] for result in cosebis_results.values()] - - rows = [start for start in range(nbins) for stop in range(start, nbins)] - columns = [stop for start in range(nbins) for stop in range(start, nbins)] - + # Place only the cuts actually scanned, using their binning metadata. + # A sparse scan must not be assigned to the first N triangular cells. snr_matrix = np.full((nbins, nbins), np.nan) pte_matrix = np.full((nbins, nbins), np.nan) - for row, column, snr, pte in zip(rows, columns, snrs, ptes): - snr_matrix[row, column] = snr - pte_matrix[row, column] = pte - - # Fill matrices from COSEBIs results - for i in range(len(left_edges)): - for j in range(i, len(right_edges)): - scale_cut = (left_edges[i], right_edges[j]) - result = cosebis_results.get(scale_cut) - - if result is not None: - snr_matrix[i, j] = np.sqrt(result["chi2_E"]) - pte_matrix[i, j] = result["pte_B"] + for scale_cut, result in cosebis_results.items(): + start, stop = bins_from_scale_cut(left_edges, right_edges, scale_cut) + snr_matrix[start, stop - 1] = np.sqrt(result["chi2_E"]) + pte_matrix[start, stop - 1] = result["pte_B"] # Create 1x2 subplot layout fig, axs = plt.subplots(1, 2, figsize=(15, 7)) @@ -1348,7 +1334,6 @@ def load_cosebis_results(path): with np.load(path, allow_pickle=False) as data: results = {} for lo, hi in data["scale_cuts"]: - suffix = f"_{lo}_{hi}" # Match the writer's scalar formatting (including integer cuts). suffix = next( k.removeprefix("En") @@ -1388,8 +1373,10 @@ def save_cosebis_results(results, output_path, fiducial_scale_cut=None, edges=No output_path : str Output .npz file path fiducial_scale_cut : tuple, optional - If results has multiple scale cuts, save only this one. - If None and results has multiple scale cuts, saves all. + If given, also save untagged entries for this fiducial cut (used by + summary consumers). All scanned cuts are always retained. + edges : tuple of arrays, optional + Integration-bin left/right edges used by the scale-cut heatmap. """ is_multi = all(isinstance(k, tuple) for k in results) diff --git a/src/sp_validation/tests/test_plot_products.py b/src/sp_validation/tests/test_plot_products.py new file mode 100644 index 00000000..6fc05039 --- /dev/null +++ b/src/sp_validation/tests/test_plot_products.py @@ -0,0 +1,278 @@ +"""Small, catalogue-free products exercise the compute/render boundary.""" + +from pathlib import Path +from types import SimpleNamespace + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import pytest + +from sp_validation import b_modes +from sp_validation.cosmo_val import diagnostic_plots as plots +from sp_validation.cosmo_val.catalog_characterization import ( + CatalogCharacterizationMixin, +) + +STYLES = {"synthetic": {"colour": "C0", "marker": "o"}} + + +@pytest.fixture(autouse=True) +def no_catalogues(monkeypatch): + from sp_validation import io + + def forbidden(*args, **kwargs): + raise AssertionError("A renderer attempted a catalogue read") + + monkeypatch.setattr(io, "open_entry", forbidden) + yield + plt.close("all") + + +def assert_figures(*paths): + for path in paths: + assert Path(path).stat().st_size > 100 + + +def psf_product(): + rng = np.random.default_rng(42) + data = { + "theta": np.geomspace(1, 100, 4), + "flat_samples": rng.normal(size=(1000, 3)), + } + for component in ("plus", "minus"): + samples = np.exp(rng.normal(-10, 0.1, (100, 4))) + data.update( + { + f"samples_{component}": samples, + f"mean_{component}": samples.mean(axis=0), + f"var_{component}": samples.var(axis=0), + f"quantiles_{component}": np.quantile(samples, [0.16, 0.84], axis=0), + } + ) + return data + + +def test_cached_psf_handlers_do_not_initialize_catalogues(tmp_path): + import importlib.util + + root = Path(__file__).resolve().parents[3] + spec = importlib.util.spec_from_file_location( + "plot_test_runner", root / "workflow/scripts/cv_runner.py" + ) + runner = importlib.util.module_from_spec(spec) + spec.loader.exec_module(runner) + cv = SimpleNamespace( + cc={"paths": {"output": str(tmp_path)}}, + treecorr_config={"min_sep": 1.0, "max_sep": 100.0, "nbins": 4}, + ) + runner.prepare_cached_psf(cv) + assert cv._rho_stat_handler.catalogs._output == str(tmp_path / "rho_tau_stats") + assert cv._tau_stat_handler.catalogs is cv._rho_stat_handler.catalogs + + +def test_objectwise_compute_uses_regression_not_plot_driver(tmp_path, monkeypatch): + from contextlib import nullcontext + + from sp_validation.cosmo_val.psf_systematics import PSFSystematicsMixin, leakage + + params = {k: SimpleNamespace(value=0.01, stderr=0.001) for k in ("a11", "a22")} + calls = [] + + def regress(*args, **kwargs): + calls.append(kwargs) + return params + + monkeypatch.setattr(leakage, "corr_2d", regress) + monkeypatch.setattr(leakage, "save_to_file", lambda *args: None) + obj = SimpleNamespace( + get_out_base=lambda *args, **kwargs: str(tmp_path / "fit"), + temporarily_read_data=lambda **kwargs: nullcontext(), + get_ellipticity_weights=lambda: ([np.zeros(4), np.zeros(4)], np.ones(4)), + _dat={"p1": np.zeros(4), "p2": np.zeros(4)}, + _params={"e1_PSF_col": "p1", "e2_PSF_col": "p2", "verbose": False}, + _stats_file=None, + ) + cv = PSFSystematicsMixin() + cv.print_cyan = lambda *args: None + result = cv._objectwise_leakage_bin(obj, "synthetic", "all", True, "lin") + assert len(calls) == 1 + assert result["aii_mean"].nominal_value == pytest.approx(0.01) + + +def test_weights_roundtrip_and_overwrite(tmp_path): + path = tmp_path / "weights.npz" + w = np.array([0.1, 0.2, 0.2, 0.5, 0.9]) + + # This is exactly the compute seam; no matplotlib histogram in the renderer. + class CV(CatalogCharacterizationMixin): + versions = ["synthetic"] + results = { + "synthetic": SimpleNamespace( + temporarily_read_data=lambda: __import__("contextlib").nullcontext() + ) + } + + def _read_shear_cols(self, *args): + return w + + counts, edges = CV().weight_histograms(nbins=3)["synthetic"] + assert counts.sum() == pytest.approx(w.sum()) + np.savez(path, counts=counts, edges=edges) + with np.load(path, allow_pickle=False) as cache: + hist = {"synthetic": (cache["counts"], cache["edges"])} + output = tmp_path / "weights.png" + output.write_text("stale") + plots.plot_weights(hist, STYLES, output) + assert_figures(output) + + +def test_footprint_every_region_from_one_map(tmp_path, monkeypatch): + from sp_validation.cosmo_val.diagnostic_products import ( + footprint_counts, + load_footprint, + save_footprint, + ) + + count_map = footprint_counts([150.0, 150.0, 20.0], [40.0, 40.0, 30.0], 4, 16) + assert count_map[count_map.valid_pixels].sum() == 3 + cache = tmp_path / "map.npz" + save_footprint(count_map, cache) + loaded = load_footprint(cache) + calls = [] + + def draw(self, data, region, outpath): + np.testing.assert_array_equal(data.valid_pixels, count_map.valid_pixels) + calls.append(region) + Path(outpath).write_bytes(b"rendered" * 100) + + monkeypatch.setattr(plots.FootprintPlotter, "plot_region", draw) + outputs = {r: tmp_path / f"{r}.png" for r in plots.FootprintPlotter._regions} + plots.plot_footprint(loaded, outputs) + assert calls == list(plots.FootprintPlotter._regions.values()) + assert_figures(*outputs.values()) + + +def test_scalar_diagnostics_and_covariance(tmp_path): + coeff = { + "synthetic": { + "tomo_bin_all": { + key: [0.001, 0.0001] + for key in ("aii_mean", "alpha_mean", "alpha_1", "alpha_0") + } + } + } + leakage, covariance = tmp_path / "leakage.png", tmp_path / "covariance.pdf" + plots.plot_objectwise_leakage(coeff, STYLES, leakage) + plots.plot_covariance(np.eye(4), covariance) + assert_figures(leakage, covariance) + + +def test_psf_fit_and_sealed_ratio_rendering(tmp_path): + data = psf_product() + path = tmp_path / "psf.npz" + np.savez(path, **data) + with np.load(path, allow_pickle=False) as cache: + products = {"synthetic": dict(cache)} + contours, samples, summary, ratio = [ + tmp_path / f"{x}.png" for x in ("contours", "samples", "summary", "ratio") + ] + plots.plot_rho_tau_fits(products, STYLES, contours, samples, summary) + gg = SimpleNamespace( + meanr=data["theta"], + xip=np.ones(4) * 1e-4, + xim=np.ones(4) * 1e-4, + varxip=np.ones(4) * 1e-12, + varxim=np.ones(4) * 1e-12, + ) + plots.plot_ratio_xi_sys_xi(products, {"synthetic": gg}, STYLES, ratio) + xi_paths = [tmp_path / f"xi{i}.png" for i in range(2)] + plots.plot_2pcf({"synthetic": gg}, STYLES, xi_paths) + assert_figures(contours, samples, summary, ratio, *xi_paths) + + +@pytest.mark.parametrize("kind,indices", [("rho", range(6)), ("tau", (0, 2, 5))]) +def test_psf_overlay_arrays(tmp_path, kind, indices): + table = {"theta": np.geomspace(1, 100, 4)} + for i in indices: + table[f"{kind}_{i}_p"] = np.ones(4) * 1e-6 + table[f"var{kind}_{i}_p"] = np.ones(4) * 1e-14 + output = tmp_path / f"{kind}.png" + plots.plot_psf_stats( + {"synthetic": table}, + STYLES, + output, + kind, + [np.eye(12) * 1e-14] if kind == "tau" else [], + ) + assert_figures(output) + + +def test_pure_eb_plot_cache_is_complete(tmp_path): + n = 3 + edges = np.geomspace(1, 100, n + 1) + theta = np.sqrt(edges[:-1] * edges[1:]) + result = dict( + theta=theta, + left_edges=edges[:-1], + right_edges=edges[1:], + theta_int=theta, + xip_int=np.ones(n) * 1e-5, + xim_int=np.ones(n) * 1e-5, + xip=np.ones(n) * 1e-5, + xim=np.ones(n) * 1e-5, + var_xip=np.ones(n) * 1e-12, + var_xim=np.ones(n) * 1e-12, + cov=np.eye(6 * n) * 1e-12, + npatch=None, + ) + result.update({k: np.ones(n) * 1e-6 for k in b_modes._EB_KEYS}) + result = b_modes.calculate_eb_statistics(result) + path = tmp_path / "sealed_eb.npz" + b_modes.save_pure_eb_results(result, path) + loaded = b_modes.load_pure_eb_results(path) + assert loaded["npatch"] is None + for key in ("theta_int", "xip_int", "xim_int", "xip", "var_xip", "cov"): + np.testing.assert_array_equal(result[key], loaded[key]) + outputs = [tmp_path / f"eb{i}.png" for i in range(4)] + b_modes.plot_integration_vs_reporting(loaded, outputs[0], "synthetic") + b_modes.plot_pure_eb_correlations(loaded, outputs[1], "synthetic") + b_modes.plot_pte_2d_heatmaps(loaded, "synthetic", outputs[2]) + b_modes.plot_eb_covariance_matrix( + loaded["cov"], "analytic", outputs[3], "synthetic" + ) + assert_figures(*outputs) + + +def test_cosebis_preserves_all_cuts_for_plotting(tmp_path): + edges = np.array([1.0, 2.0, 4.0, 8.0]) + results = { + (lo, hi): dict( + En=np.array([1.0, 2.0]), + Bn=np.array([0.1, 0.2]), + cov=np.eye(4), + chi2_E=5.0, + chi2_B=0.05, + pte_B=0.9, + ) + for i, lo in enumerate(edges[:-1]) + for hi in edges[i + 1 :] + } + path = tmp_path / "sealed_cosebis.npz" + b_modes.save_cosebis_results(results, path, (1.0, 8.0), (edges[:-1], edges[1:])) + loaded, loaded_edges = b_modes.load_cosebis_results(path) + assert loaded.keys() == results.keys() + np.testing.assert_array_equal(loaded_edges, (edges[:-1], edges[1:])) + with np.load(path, allow_pickle=False) as cache: + assert float(cache["pte_B"]) == 0.9 # Existing summary consumer. + outputs = [tmp_path / f"cosebis{i}.png" for i in range(3)] + b_modes.plot_cosebis_modes(loaded[1.0, 8.0], "synthetic", outputs[0]) + b_modes.plot_cosebis_covariance_matrix( + loaded[1.0, 8.0], "synthetic", "analytic", outputs[1] + ) + b_modes.plot_cosebis_scale_cut_heatmap( + loaded, loaded_edges, "synthetic", outputs[2] + ) + assert_figures(*outputs) diff --git a/src/sp_validation/tests/test_plot_workflow.py b/src/sp_validation/tests/test_plot_workflow.py new file mode 100644 index 00000000..42c25401 --- /dev/null +++ b/src/sp_validation/tests/test_plot_workflow.py @@ -0,0 +1,145 @@ +"""Dry-run the real plotting DAG with caches present and no catalogues. + +Set SPV_TEST_SNAKEMAKE to a host Snakemake executable if it isn't in the image. +The test never launches a workflow job or opens a real science product. +""" + +import importlib.util +import os +import shutil +import subprocess +from pathlib import Path + +import pytest +import yaml + +ROOT = Path(__file__).resolve().parents[3] +PLOT_RULES = [ + "cosmo_val_plots", + "cv_plot_rho_stats", + "cv_plot_tau_stats", + "cv_plot_rho_tau_fits", + "cv_plot_footprints", + "cv_plot_objectwise_leakage", + "cv_plot_weights", + "cv_plot_2pcf", + "cv_ratio_xi_sys_xi", + "cv_plot_pure_eb", + "cv_plot_cosebis", + "cv_plot_pseudo_cl", + "covariance_plot", +] + + +@pytest.mark.parametrize("pseudo_cl", [False, True]) +def test_cached_plot_dag_has_no_compute_jobs(tmp_path, monkeypatch, pseudo_cl): + executable = os.environ.get("SPV_TEST_SNAKEMAKE") or shutil.which("snakemake") + if not executable: + pytest.skip("Snakemake not installed; set SPV_TEST_SNAKEMAKE") + monkeypatch.chdir(tmp_path) + monkeypatch.setenv("COSMO_VAL", str(tmp_path / "products")) + monkeypatch.setenv("COSMO_INFERENCE", str(tmp_path / "inference")) + monkeypatch.setenv("PYTHONUNBUFFERED", "1") + # common.configure expects this fixed relative cosmology file at parse time. + cosmology = tmp_path / "results/cosmology/planck18.json" + cosmology.parent.mkdir(parents=True) + shutil.copy(ROOT / "papers/cosmo_val/results/cosmology/planck18.json", cosmology) + config = yaml.safe_load((ROOT / "papers/cosmo_val/config/config.yaml").read_text()) + config["versions"] = [ + "SP_v1.4.6.3" + ] # Public fixture entry; no real input is opened. + config["cosmo_val"]["include_pseudo_cl"] = pseudo_cl + configfile = tmp_path / "config.yaml" + configfile.write_text(yaml.safe_dump(config)) + + spec = importlib.util.spec_from_file_location( + "plot_test_common", ROOT / "workflow/common.py" + ) + common = importlib.util.module_from_spec(spec) + spec.loader.exec_module(common) + # Set the path-helper context without importing Snakemake in the image; + # the subprocess below runs the real configure() in the host environment. + common.CATALOG_CONFIG = yaml.safe_load(Path(common.CAT_CONFIG).read_text()) + common.FIDUCIAL = config["fiducial"] + common.DEFAULT_MASK_SUFFIX = ( + "_masked" if config["covariance"].get("default_masked", False) else "" + ) + namespace = dict(vars(common), config=config) + # Execute only the ordinary Python path helpers, before the first rule. + twopoint = (ROOT / "workflow/rules/twopoint.smk").read_text().split("rule xi:")[0] + exec(twopoint, namespace) + diagnostics = ( + (ROOT / "workflow/rules/cosmo_val.smk") + .read_text() + .split("rule cv_plot_rho_stats:")[0] + ) + exec(diagnostics, namespace) + ns = namespace + version = config["versions"][0] + products = [ + *ns["CV_PSF_PRODUCTS"], + *ns["CV_FOOTPRINT_MAPS"], + str(common.COSMO_VAL / "weight_hist.npz"), + str(common.COSMO_VAL / "leakage_coefficients.json"), + common.cv_rho_stats(version, ns["CV_FIDUCIAL"]), + common.cv_tau_stats(version, ns["CV_FIDUCIAL"]), + *ns["cv_tau_covariances"](), + ns["cv_xi_sacc"](version, "reporting"), + ns["cv_pure_eb_npz"](version), + ns["cv_cosebis_npz"](version), + ns["cv_xi_cov"](version), + ns["cv_xi_cov_integration"](version), + ] + if pseudo_cl: + products += [ + ns["cv_pseudo_cl_analysis_sacc"](version), + ns["cv_pseudo_cl_cov"](version), + ] + for product in products: + path = Path(product) + path.parent.mkdir(parents=True, exist_ok=True) + path.touch() + + command = [ + executable, + "-s", + str(ROOT / "workflow/Snakefile"), + "--configfile", + str(configfile), + "--cores", + "1", + "--dry-run", + "--forceall", + "--workflow-profile", + "none", + "--rerun-triggers", + "mtime", + "--allowed-rules", + *PLOT_RULES, + "--", + "cosmo_val_plots", + ] + result = subprocess.run(command, text=True, capture_output=True, timeout=120) + log = result.stdout + result.stderr + assert result.returncode == 0, log + assert "cv_plot_weights" in log + assert "cv_plot_footprints" in log + assert "cv_plot_pure_eb" in log + assert "cv_plot_cosebis" in log + assert ("rule cv_plot_pseudo_cl:" in log) == pseudo_cl + for forbidden in ( + "cv_weights", + "cv_footprints", + "cv_rho_tau_fits", + "cv_cosebis", + "cv_pure_eb", + "xi", + "rho_tau_stats", + "covariance_run", + ): + assert f"rule {forbidden}:" not in log + # A plot-only restyle must fail closed, not silently rebuild from a catalogue. + (common.COSMO_VAL / "weight_hist.npz").unlink() + missing = subprocess.run(command, text=True, capture_output=True, timeout=120) + assert missing.returncode != 0 + assert "MissingInputException" in missing.stdout + missing.stderr diff --git a/workflow/README.md b/workflow/README.md index e3e15c6c..8336d02c 100644 --- a/workflow/README.md +++ b/workflow/README.md @@ -159,6 +159,39 @@ ln -s /n17data/cdaley/unions/code/sp_validation/cosmo_inference/data/mask **Caveat:** `rerun-triggers: code` watches rule bodies and `script:` files, not `src/`. Editing a module under `src/` does not by itself mark outputs stale — force with `-F` or `--forcerun `. +The `cosmo_val` plot rules are an exception: their plotting modules are explicit inputs, so editing those modules invalidates figures without invalidating computations. + +### Re-render the cosmo_val figures from cache + +`cosmo_val_all` runs the computations and every companion plot. +`cosmo_val_plots` gathers only figures, including the reporting/integration covariance correlation plots. +Plots read the sealed SACC/NPZ products, sparse footprint counts, weight histograms and PSF/leakage summaries; no plot rule reads a catalogue or fits a model. +Every figure is a declared output and is overwritten when its plot rule runs. + +After a full run, use the same config and output roots to render changed plots: + +```bash +snakemake --profile workflow/profiles/default -s workflow/Snakefile \ + --configfile --cores 4 cosmo_val_plots +``` + +To force every figure to be restyled, while **refusing to run any computation even if a cache is missing**, restrict the allowed rules: + +```bash +plot_rules=(cosmo_val_plots cv_plot_rho_stats cv_plot_tau_stats \ + cv_plot_rho_tau_fits cv_plot_footprints cv_plot_objectwise_leakage \ + cv_plot_weights cv_plot_2pcf cv_ratio_xi_sys_xi cv_plot_pure_eb \ + cv_plot_cosebis cv_plot_pseudo_cl covariance_plot) +snakemake --profile workflow/profiles/default -s workflow/Snakefile \ + --configfile --cores 4 --forceall \ + --allowed-rules "${plot_rules[@]}" -- cosmo_val_plots +``` + +Run this on an allocation, using the host's container/bind settings as usual. +Without the rule restriction, requesting a missing data product can schedule its producer. +A pre-split run needs one compute pass to create the new caches: old pure-E/B NPZs lack the comparison vectors and old COSEBIs NPZs contain only the fiducial cut. +That pass derives cosmological products only from the already-sealed integration SACC; no unsealed vector is added for plotting. +The objectwise rule now calls the spin-2 regression directly rather than the external all-in-one exploratory driver, so it produces the suite's leakage-coefficient comparison rather than unrelated 1D-fit/scatter side effects. To reproduce a run from the image alone, opt out: diff --git a/workflow/scripts/cv_rho_tau_fits.py b/workflow/scripts/cv_rho_tau_fits.py index d2e49378..bcdbff82 100644 --- a/workflow/scripts/cv_rho_tau_fits.py +++ b/workflow/scripts/cv_rho_tau_fits.py @@ -12,6 +12,12 @@ data = dict( theta=cv.psf_fitter.rho_stat_handler.rho_stats["theta"], flat_samples=flat_samples, + theta_edges=np.geomspace( + cv.treecorr_config["min_sep"], + cv.treecorr_config["max_sep"], + cv.treecorr_config["nbins"] + 1, + ), + theta_unit="arcmin", ) for component, samples in (("plus", plus), ("minus", minus)): data[f"samples_{component}"] = samples From 95ce64aff9ab768f8b9374289d2bf9f7533b9ad4 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 24 Oct 2026 13:28:55 -0400 Subject: [PATCH 04/18] Require inherited blind provenance in cosmological NPZ consumers --- src/sp_validation/b_modes.py | 69 +++++++++++-------- src/sp_validation/cosmo_val/cosebis.py | 5 +- src/sp_validation/cosmo_val/pure_eb.py | 6 +- src/sp_validation/sacc_io.py | 30 ++++++-- src/sp_validation/tests/test_b_modes.py | 6 +- .../tests/test_blind_consumers.py | 16 +++++ src/sp_validation/tests/test_plot_products.py | 16 +++-- workflow/rules/cosmo_val.smk | 2 + workflow/scripts/cv_cosebis.py | 4 +- workflow/scripts/cv_plot_cosebis.py | 2 +- workflow/scripts/cv_plot_pure_eb.py | 2 +- workflow/scripts/cv_pure_eb.py | 2 +- workflow/scripts/cv_summarize_bmodes.py | 4 +- 13 files changed, 117 insertions(+), 47 deletions(-) diff --git a/src/sp_validation/b_modes.py b/src/sp_validation/b_modes.py index d4cc548f..70e6f2d8 100644 --- a/src/sp_validation/b_modes.py +++ b/src/sp_validation/b_modes.py @@ -1285,7 +1285,7 @@ def plot_cosebis_modes(results, version, output_path, fiducial_scale_cut=None): plt.savefig(output_path, dpi=300, bbox_inches="tight") -def save_pure_eb_results(results, output_path): +def save_pure_eb_results(results, output_path, *, derived_from): """ Save pure E/B mode data vectors and covariance to .npz. @@ -1312,40 +1312,45 @@ def save_pure_eb_results(results, output_path): if results["npatch"] is not None: save_dict["npatch"] = np.array(results["npatch"]) - np.savez(output_path, **save_dict) + from . import sacc_io + + sacc_io.save_npz(save_dict, output_path, derived_from=derived_from) print(f"Saved pure E/B results to {output_path}") -def load_pure_eb_results(path): - """Read a sealed plot cache without transforms or catalogue access.""" - with np.load(path, allow_pickle=False) as data: - results = {k: data[k] for k in data.files if not k.startswith("pte_matrices_")} - results["pte_matrices"] = { - k.removeprefix("pte_matrices_"): data[k] - for k in data.files - if k.startswith("pte_matrices_") - } +def load_pure_eb_results(path, *, blind): + """Read a sealed plot cache only after checking the expected blind.""" + from . import sacc_io + + data = sacc_io.load_npz(path, blind=blind) + results = {k: v for k, v in data.items() if not k.startswith("pte_matrices_")} + results["pte_matrices"] = { + k.removeprefix("pte_matrices_"): v + for k, v in data.items() + if k.startswith("pte_matrices_") + } results["npatch"] = int(results["npatch"]) if "npatch" in results else None return results -def load_cosebis_results(path): - """Read all cached scale cuts and their integration-bin edges.""" - with np.load(path, allow_pickle=False) as data: - results = {} - for lo, hi in data["scale_cuts"]: - # Match the writer's scalar formatting (including integer cuts). - suffix = next( - k.removeprefix("En") - for k in data.files - if k.startswith("En_") - and tuple(map(float, k[3:].split("_"))) == (lo, hi) - ) - results[lo, hi] = { - k: data[k + suffix] - for k in ("En", "Bn", "cov", "chi2_E", "chi2_B", "pte_B") - } - edges = data["left_edges"], data["right_edges"] +def load_cosebis_results(path, *, blind): + """Read cached scale cuts after checking their inherited blind stamp.""" + from . import sacc_io + + data = sacc_io.load_npz(path, blind=blind) + results = {} + for lo, hi in data["scale_cuts"]: + # Match the writer's scalar formatting (including integer cuts). + suffix = next( + k.removeprefix("En") + for k in data + if k.startswith("En_") and tuple(map(float, k[3:].split("_"))) == (lo, hi) + ) + results[lo, hi] = { + k: data[k + suffix] + for k in ("En", "Bn", "cov", "chi2_E", "chi2_B", "pte_B") + } + edges = data["left_edges"], data["right_edges"] return results, edges @@ -1361,7 +1366,9 @@ def _cosebis_result_to_dict(r, suffix=""): } -def save_cosebis_results(results, output_path, fiducial_scale_cut=None, edges=None): +def save_cosebis_results( + results, output_path, fiducial_scale_cut=None, edges=None, *, derived_from +): """ Save COSEBIs data vectors and covariance to .npz. @@ -1403,7 +1410,9 @@ def save_cosebis_results(results, output_path, fiducial_scale_cut=None, edges=No if edges is not None: save_dict["left_edges"], save_dict["right_edges"] = edges - np.savez(output_path, **save_dict) + from . import sacc_io + + sacc_io.save_npz(save_dict, output_path, derived_from=derived_from) print(f"Saved COSEBIs results to {output_path}") diff --git a/src/sp_validation/cosmo_val/cosebis.py b/src/sp_validation/cosmo_val/cosebis.py index bea1df16..b2714e6e 100644 --- a/src/sp_validation/cosmo_val/cosebis.py +++ b/src/sp_validation/cosmo_val/cosebis.py @@ -310,6 +310,9 @@ def plot_cosebis( ) save_cosebis_results( - bin_results, bin_stub + "_data.npz", fiducial_scale_cut + bin_results, + bin_stub + "_data.npz", + fiducial_scale_cut, + derived_from=[self.xi_parts[version, "integration"]], ) self._cosebis_results[version] = results_tomo diff --git a/src/sp_validation/cosmo_val/pure_eb.py b/src/sp_validation/cosmo_val/pure_eb.py index 793fa826..96a80faa 100644 --- a/src/sp_validation/cosmo_val/pure_eb.py +++ b/src/sp_validation/cosmo_val/pure_eb.py @@ -287,5 +287,9 @@ def plot_pure_eb( ) # Save data products and store on instance - save_pure_eb_results(version_results, out_stub + "_data.npz") + save_pure_eb_results( + version_results, + out_stub + "_data.npz", + derived_from=[self.xi_parts[version, "integration"]], + ) self._pure_eb_results[version][bin_key] = version_results diff --git a/src/sp_validation/sacc_io.py b/src/sp_validation/sacc_io.py index c4bdc5c7..f20f3b3f 100644 --- a/src/sp_validation/sacc_io.py +++ b/src/sp_validation/sacc_io.py @@ -1010,6 +1010,31 @@ def seal(s, blind, theory=None): return out +def derived_stamp(parts): + """The common stamp inherited by a derivation, for SACC and NPZ alike.""" + stamps = {stamp(p) for p in parts} + if len(stamps) != 1: + raise ValueError(f"input parts carry different blinds: {sorted(stamps)}") + return stamps.pop() + + +def save_npz(arrays, path, *, derived_from): + """Save a cosmological array cache with the same inheritance as SACC.""" + np.savez(path, **{**arrays, STAMP_KEY: derived_stamp(derived_from)}) + + +def load_npz(path, *, blind): + """Check explicit provenance before reading any cached signal arrays.""" + from types import SimpleNamespace + + with np.load(path, allow_pickle=False) as data: + if STAMP_KEY not in data: + raise ValueError("cosmological NPZ is missing its blind stamp") + part = SimpleNamespace(metadata={STAMP_KEY: str(data[STAMP_KEY].item())}) + expect(part, blind) + return {k: data[k] for k in data.files if k != STAMP_KEY} + + def save(s, path, *, blind=None, theory=None, derived_from=None): """Write ``s`` to ``path`` (FITS) and return what was written. @@ -1028,10 +1053,7 @@ def save(s, path, *, blind=None, theory=None, derived_from=None): signal never leaves that function. """ if derived_from is not None: - stamps = {stamp(p) for p in derived_from} - if len(stamps) != 1: - raise ValueError(f"input parts carry different blinds: {sorted(stamps)}") - (name,) = stamps + name = derived_stamp(derived_from) out = s.copy() out.metadata[STAMP_KEY] = name if blind is not None: diff --git a/src/sp_validation/tests/test_b_modes.py b/src/sp_validation/tests/test_b_modes.py index 2ed8a341..01c30aae 100644 --- a/src/sp_validation/tests/test_b_modes.py +++ b/src/sp_validation/tests/test_b_modes.py @@ -715,7 +715,11 @@ def test_pure_eb_npz_carries_what_the_summary_reads(tmp_path): results = b_modes.calculate_eb_statistics(results) out = tmp_path / "pure_eb_data.npz" - b_modes.save_pure_eb_results(results, str(out)) + from types import SimpleNamespace + + b_modes.save_pure_eb_results( + results, str(out), derived_from=[SimpleNamespace(metadata={"blind": "none"})] + ) saved = np.load(out) for stat in ("xip_B", "xim_B", "combined"): diff --git a/src/sp_validation/tests/test_blind_consumers.py b/src/sp_validation/tests/test_blind_consumers.py index 3d3b7a8d..57e73261 100644 --- a/src/sp_validation/tests/test_blind_consumers.py +++ b/src/sp_validation/tests/test_blind_consumers.py @@ -11,6 +11,22 @@ SCRIPTS = Path(__file__).resolve().parents[3] / "workflow" / "scripts" +@pytest.mark.parametrize("loader", ["load_pure_eb_results", "load_cosebis_results"]) +@pytest.mark.parametrize("stamp", [None, "wrong"]) +def test_npz_consumers_refuse_missing_or_wrong_blind(tmp_path, loader, stamp): + import numpy as np + + from sp_validation import b_modes + + path = tmp_path / "cosmological.npz" + data = dict(scale_cuts=np.empty((0, 2)), left_edges=[], right_edges=[]) + if stamp is not None: + data["blind"] = stamp + np.savez(path, **data) + with pytest.raises(ValueError, match="stamp"): + getattr(b_modes, loader)(path, blind="expected") + + @pytest.fixture def part_path(tmp_path): part = sacc_io.new_sacc({0: ([0.1, 1.0], [1.0, 1.0])}) diff --git a/src/sp_validation/tests/test_plot_products.py b/src/sp_validation/tests/test_plot_products.py index 6fc05039..c5b7e2dc 100644 --- a/src/sp_validation/tests/test_plot_products.py +++ b/src/sp_validation/tests/test_plot_products.py @@ -231,8 +231,11 @@ def test_pure_eb_plot_cache_is_complete(tmp_path): result.update({k: np.ones(n) * 1e-6 for k in b_modes._EB_KEYS}) result = b_modes.calculate_eb_statistics(result) path = tmp_path / "sealed_eb.npz" - b_modes.save_pure_eb_results(result, path) - loaded = b_modes.load_pure_eb_results(path) + parent = SimpleNamespace(metadata={"blind": "toy"}) + b_modes.save_pure_eb_results(result, path, derived_from=[parent]) + loaded = b_modes.load_pure_eb_results(path, blind="toy") + with np.load(path, allow_pickle=False) as cache: + assert cache["blind"].item() == "toy" assert loaded["npatch"] is None for key in ("theta_int", "xip_int", "xim_int", "xip", "var_xip", "cov"): np.testing.assert_array_equal(result[key], loaded[key]) @@ -261,8 +264,13 @@ def test_cosebis_preserves_all_cuts_for_plotting(tmp_path): for hi in edges[i + 1 :] } path = tmp_path / "sealed_cosebis.npz" - b_modes.save_cosebis_results(results, path, (1.0, 8.0), (edges[:-1], edges[1:])) - loaded, loaded_edges = b_modes.load_cosebis_results(path) + parent = SimpleNamespace(metadata={"blind": "toy"}) + b_modes.save_cosebis_results( + results, path, (1.0, 8.0), (edges[:-1], edges[1:]), derived_from=[parent] + ) + loaded, loaded_edges = b_modes.load_cosebis_results(path, blind="toy") + with np.load(path, allow_pickle=False) as cache: + assert cache["blind"].item() == "toy" assert loaded.keys() == results.keys() np.testing.assert_array_equal(loaded_edges, (edges[:-1], edges[1:])) with np.load(path, allow_pickle=False) as cache: diff --git a/workflow/rules/cosmo_val.smk b/workflow/rules/cosmo_val.smk index eb37511e..6a9a4158 100644 --- a/workflow/rules/cosmo_val.smk +++ b/workflow/rules/cosmo_val.smk @@ -540,6 +540,7 @@ rule cv_plot_pure_eb: **cv_pure_eb_figures("{version}"), params: version="{version}", + blind=lambda w: blind_of(w.version), fiducial_scale_cut=CV["fiducial_scale_cut"], script: "../scripts/cv_plot_pure_eb.py" @@ -554,6 +555,7 @@ rule cv_plot_cosebis: **cv_cosebis_figures("{version}"), params: version="{version}", + blind=lambda w: blind_of(w.version), fiducial_scale_cut=CV["fiducial_scale_cut"], script: "../scripts/cv_plot_cosebis.py" diff --git a/workflow/scripts/cv_cosebis.py b/workflow/scripts/cv_cosebis.py index 97308d4d..1f58a304 100644 --- a/workflow/scripts/cv_cosebis.py +++ b/workflow/scripts/cv_cosebis.py @@ -41,7 +41,9 @@ fiducial_key = find_conservative_scale_cut_key(results, fiducial_scale_cut) fiducial = results[fiducial_key] -save_cosebis_results(results, snakemake.output["npz"], fiducial_scale_cut, edges) +save_cosebis_results( + results, snakemake.output["npz"], fiducial_scale_cut, edges, derived_from=[part] +) s = cosebis_to_sacc({0: sacc_io.get_nz(part, 0)}, part.metadata, fiducial, fiducial_key) sacc_io.save(s, snakemake.output["sacc"], derived_from=[part]) diff --git a/workflow/scripts/cv_plot_cosebis.py b/workflow/scripts/cv_plot_cosebis.py index 7e3f10da..0d2085a5 100644 --- a/workflow/scripts/cv_plot_cosebis.py +++ b/workflow/scripts/cv_plot_cosebis.py @@ -9,7 +9,7 @@ ) p = snakemake.params -results, edges = load_cosebis_results(snakemake.input.npz) +results, edges = load_cosebis_results(snakemake.input.npz, blind=p.blind) cut = tuple(p.fiducial_scale_cut) fiducial = results[find_conservative_scale_cut_key(results, cut)] plot_cosebis_modes(fiducial, p.version, snakemake.output.figure_modes, cut) diff --git a/workflow/scripts/cv_plot_pure_eb.py b/workflow/scripts/cv_plot_pure_eb.py index fec9a0e5..b0c2b9a1 100644 --- a/workflow/scripts/cv_plot_pure_eb.py +++ b/workflow/scripts/cv_plot_pure_eb.py @@ -10,7 +10,7 @@ ) p = snakemake.params -results = load_pure_eb_results(snakemake.input.npz) +results = load_pure_eb_results(snakemake.input.npz, blind=p.blind) kwargs = dict( fiducial_xip_scale_cut=tuple(p.fiducial_scale_cut), fiducial_xim_scale_cut=tuple(p.fiducial_scale_cut), diff --git a/workflow/scripts/cv_pure_eb.py b/workflow/scripts/cv_pure_eb.py index cf561e08..50a65fbc 100644 --- a/workflow/scripts/cv_pure_eb.py +++ b/workflow/scripts/cv_pure_eb.py @@ -41,7 +41,7 @@ ) results = calculate_eb_statistics(results) -save_pure_eb_results(results, snakemake.output["npz"]) +save_pure_eb_results(results, snakemake.output["npz"], derived_from=[part]) sacc_io.save( pure_eb_to_sacc( diff --git a/workflow/scripts/cv_summarize_bmodes.py b/workflow/scripts/cv_summarize_bmodes.py index 8175e0d6..f91648d2 100644 --- a/workflow/scripts/cv_summarize_bmodes.py +++ b/workflow/scripts/cv_summarize_bmodes.py @@ -26,7 +26,7 @@ for i, version in enumerate(p["versions"]): row = {} - pure_eb = np.load(snakemake.input["pure_eb"][i]) + pure_eb = sacc_io.load_npz(snakemake.input["pure_eb"][i], blind=p["blinds"][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"): @@ -41,7 +41,7 @@ # The COSEBIs .npz is written at the fiducial cut, so its PTE is the one # this table wants. - cosebis = np.load(snakemake.input["cosebis"][i]) + cosebis = sacc_io.load_npz(snakemake.input["cosebis"][i], blind=p["blinds"][i]) row["COSEBIS"] = float(cosebis["pte_B"]) cov_methods.add("COSEBIs: analytic (OneCovariance ξ± through the COSEBIs kernel)") From d23257fe7f16510bb60522496a3aed385caaa388 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 24 Oct 2026 13:40:09 -0400 Subject: [PATCH 05/18] Persist COSEBIs selected bins and place heatmaps from those indices --- src/sp_validation/b_modes.py | 14 +++-- src/sp_validation/tests/test_plot_products.py | 55 ++++++++++++++++++- 2 files changed, 62 insertions(+), 7 deletions(-) diff --git a/src/sp_validation/b_modes.py b/src/sp_validation/b_modes.py index 70e6f2d8..0e57210c 100644 --- a/src/sp_validation/b_modes.py +++ b/src/sp_validation/b_modes.py @@ -501,6 +501,7 @@ def cosebis_scan_from_xi( "pte_B": 1 - stats.chi2.cdf(chi2_B, nmodes), "scale_cut": (min_theta, max_theta), "inds": inds, + "bin_indices": np.array([start_bin, stop_bin]), } print(f"COSEBIs results [{min_theta:.1f}-{max_theta:.1f} arcmin]:") @@ -925,8 +926,8 @@ def plot_cosebis_scale_cut_heatmap( # A sparse scan must not be assigned to the first N triangular cells. snr_matrix = np.full((nbins, nbins), np.nan) pte_matrix = np.full((nbins, nbins), np.nan) - for scale_cut, result in cosebis_results.items(): - start, stop = bins_from_scale_cut(left_edges, right_edges, scale_cut) + for result in cosebis_results.values(): + start, stop = result["bin_indices"] snr_matrix[start, stop - 1] = np.sqrt(result["chi2_E"]) pte_matrix[start, stop - 1] = result["pte_B"] @@ -983,10 +984,10 @@ def plot_cosebis_scale_cut_heatmap( # Add fiducial scale cut cross-hatching if provided if fiducial_scale_cut is not None: - min_scale, max_scale = fiducial_scale_cut - start_bin, stop_bin = bins_from_edges( - left_edges, right_edges, min_scale, max_scale + fiducial_key = find_conservative_scale_cut_key( + cosebis_results, fiducial_scale_cut ) + start_bin, stop_bin = cosebis_results[fiducial_key]["bin_indices"] if stop_bin > start_bin and start_bin < nbins and stop_bin > 0: # Add cross-hatching at the fiducial scale cut matrix element @@ -1348,7 +1349,7 @@ def load_cosebis_results(path, *, blind): ) results[lo, hi] = { k: data[k + suffix] - for k in ("En", "Bn", "cov", "chi2_E", "chi2_B", "pte_B") + for k in ("En", "Bn", "cov", "chi2_E", "chi2_B", "pte_B", "bin_indices") } edges = data["left_edges"], data["right_edges"] return results, edges @@ -1363,6 +1364,7 @@ def _cosebis_result_to_dict(r, suffix=""): f"chi2_E{suffix}": np.array(r["chi2_E"]), f"chi2_B{suffix}": np.array(r["chi2_B"]), f"pte_B{suffix}": np.array(r["pte_B"]), + f"bin_indices{suffix}": np.asarray(r["bin_indices"]), } diff --git a/src/sp_validation/tests/test_plot_products.py b/src/sp_validation/tests/test_plot_products.py index c5b7e2dc..2b109c49 100644 --- a/src/sp_validation/tests/test_plot_products.py +++ b/src/sp_validation/tests/test_plot_products.py @@ -249,6 +249,58 @@ def test_pure_eb_plot_cache_is_complete(tmp_path): assert_figures(*outputs) +def test_cosebis_transform_caches_selected_bins(tmp_path, monkeypatch): + from cosmo_numba.B_modes import cosebis + + class SyntheticKernel: + def __init__(self, **kwargs): + pass + + def cosebis_from_xipm(self, theta, xip, xim, **kwargs): + return np.ones(2), np.zeros(2) + + def cosebis_covariance_from_xipm_covariance(self, theta, cov): + return np.eye(4) + + monkeypatch.setattr(cosebis, "COSEBIS", SyntheticKernel) + edges = b_modes.log_bin_edges(0.08, 300, 1000) + theta = np.sqrt(edges[0] * edges[1]) + results = b_modes.cosebis_scan_from_xi( + theta, + np.zeros(1000), + np.zeros(1000), + np.eye(2000), + *edges, + nmodes=2, + scale_cuts=[(12, 83)], + ) + np.testing.assert_array_equal(results[12, 83]["bin_indices"], [609, 843]) + path = tmp_path / "cuts.npz" + b_modes.save_cosebis_results( + results, + path, + (12, 83), + edges, + derived_from=[SimpleNamespace(metadata={"blind": "toy"})], + ) + loaded, _ = b_modes.load_cosebis_results(path, blind="toy") + np.testing.assert_array_equal(loaded[12, 83]["bin_indices"], [609, 843]) + + +def test_cosebis_heatmap_uses_stored_inward_bins(tmp_path, monkeypatch): + # Nearest edges would put this at (0, 2); the transform selected (1, 2). + result = {(1.1, 7.9): dict(chi2_E=4.0, pte_B=0.7, bin_indices=[1, 2])} + edges = (np.array([1.0, 2.0, 4.0]), np.array([2.0, 4.0, 8.0])) + b_modes.plot_cosebis_scale_cut_heatmap( + result, edges, "synthetic", tmp_path / "cut.png", (1.1, 7.9) + ) + ax = plt.gcf().axes[1] + matrix = ax.images[0].get_array() + assert matrix[1, 1] == pytest.approx(0.7) + assert np.ma.count(matrix) == 1 + assert ax.patches[0].get_xy() == (1, 1) + + def test_cosebis_preserves_all_cuts_for_plotting(tmp_path): edges = np.array([1.0, 2.0, 4.0, 8.0]) results = { @@ -259,9 +311,10 @@ def test_cosebis_preserves_all_cuts_for_plotting(tmp_path): chi2_E=5.0, chi2_B=0.05, pte_B=0.9, + bin_indices=np.array([i, j]), ) for i, lo in enumerate(edges[:-1]) - for hi in edges[i + 1 :] + for j, hi in enumerate(edges[i + 1 :], start=i + 1) } path = tmp_path / "sealed_cosebis.npz" parent = SimpleNamespace(metadata={"blind": "toy"}) From ca55e1a460e47d6db8d6927fd46ec29374b60954 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 24 Oct 2026 13:53:05 -0400 Subject: [PATCH 06/18] Require every plot cache directly even when figures already exist --- src/sp_validation/tests/test_plot_workflow.py | 46 +++++++++++++++++-- workflow/rules/cosmo_val.smk | 21 +++++++++ 2 files changed, 62 insertions(+), 5 deletions(-) diff --git a/src/sp_validation/tests/test_plot_workflow.py b/src/sp_validation/tests/test_plot_workflow.py index 42c25401..0aee0c3c 100644 --- a/src/sp_validation/tests/test_plot_workflow.py +++ b/src/sp_validation/tests/test_plot_workflow.py @@ -138,8 +138,44 @@ def test_cached_plot_dag_has_no_compute_jobs(tmp_path, monkeypatch, pseudo_cl): "covariance_run", ): assert f"rule {forbidden}:" not in log - # A plot-only restyle must fail closed, not silently rebuild from a catalogue. - (common.COSMO_VAL / "weight_hist.npz").unlink() - missing = subprocess.run(command, text=True, capture_output=True, timeout=120) - assert missing.returncode != 0 - assert "MissingInputException" in missing.stdout + missing.stderr + # Existing PNGs must never hide a missing cache when allowed-rules prunes + # a plot job. Populate *all* figures, then remove each cache in turn. + figures = [ + *ns["CV_FOOTPRINT_FIGURES"], + *ns["cv_pure_eb_figures"](version).values(), + *ns["cv_cosebis_figures"](version).values(), + *[ + str(common.COSMO_VAL / name) + for name in ( + "weight_hist.png", + "leakage_coefficients.png", + "ratio_xi_sys_xi.png", + "xi_pm_tomography_False.png", + "xi_pm_theta_tomography_False.png", + "rho_tau_stats/rho_stats.png", + "rho_tau_stats/tau_stats.png", + "rho_tau_stats/contours_tau_stat.png", + "rho_tau_stats/xi_psf_sys_tomography_False_sample.png", + "rho_tau_stats/xi_psf_sys_tomography_False_mean_and_std.png", + ) + ], + ns["cv_xi_cov"](version).removesuffix(".dat") + "_corr_plot.pdf", + ns["cv_xi_cov_integration"](version).removesuffix(".dat") + "_corr_plot.pdf", + ] + if pseudo_cl: + figures += [ + str(common.COSMO_VAL / f"cell_{name}.png") for name in ("ee", "eb", "bb") + ] + for figure in figures: + path = Path(figure) + path.parent.mkdir(parents=True, exist_ok=True) + path.write_text("existing figure") + for product in products: + path = Path(product) + path.unlink() + missing = subprocess.run(command, text=True, capture_output=True, timeout=120) + log = missing.stdout + missing.stderr + assert missing.returncode != 0, f"Missing cache {product} was ignored:\n{log}" + assert "MissingInputException" in log + assert "existing figure" == Path(figures[0]).read_text() + path.touch() diff --git a/workflow/rules/cosmo_val.smk b/workflow/rules/cosmo_val.smk index 6a9a4158..2f0194e3 100644 --- a/workflow/rules/cosmo_val.smk +++ b/workflow/rules/cosmo_val.smk @@ -664,6 +664,27 @@ rule cosmo_val_plots: list(cv_pseudo_cl_figures().values()) if CV.get("include_pseudo_cl", False) else [], [cv_xi_cov(v).removesuffix(".dat") + "_corr_plot.pdf" for v in CV_VERSIONS], [cv_xi_cov_integration(v).removesuffix(".dat") + "_corr_plot.pdf" for v in CV_VERSIONS], + # Direct dependencies are intentional. With --allowed-rules Snakemake + # can prune a plot whose figure exists even when its cache is missing. + # The aggregate must validate caches independently of those plot jobs. + caches=[ + *rules.cv_plot_rho_stats.input.stats, + *rules.cv_plot_tau_stats.input.stats, + *rules.cv_plot_tau_stats.input.covariances, + *rules.cv_plot_rho_tau_fits.input.products, + *rules.cv_plot_footprints.input.maps, + rules.cv_plot_objectwise_leakage.input.coefficients, + rules.cv_plot_weights.input.histograms, + *rules.cv_plot_2pcf.input.xi, + *[cv_pure_eb_npz(v) for v in CV_VERSIONS], + *[cv_cosebis_npz(v) for v in CV_VERSIONS], + *[cv_xi_cov(v) for v in CV_VERSIONS], + *[cv_xi_cov_integration(v) for v in CV_VERSIONS], + *([cv_pseudo_cl_analysis_sacc(v) for v in CV_VERSIONS] + if CV.get("include_pseudo_cl", False) else []), + *([cv_pseudo_cl_cov(v) for v in CV_VERSIONS] + if CV.get("include_pseudo_cl", False) else []), + ], rule cosmo_val_all: From b33c06e661475cc98df521a1f3cb5de32531652e Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 24 Oct 2026 14:04:29 -0400 Subject: [PATCH 07/18] Remove retired reporting-parts seam and test NPZ consumers instead --- .../tests/test_blind_consumers.py | 22 +++---------------- workflow/scripts/cv_runner.py | 9 -------- 2 files changed, 3 insertions(+), 28 deletions(-) diff --git a/src/sp_validation/tests/test_blind_consumers.py b/src/sp_validation/tests/test_blind_consumers.py index 57e73261..2cb4fcb7 100644 --- a/src/sp_validation/tests/test_blind_consumers.py +++ b/src/sp_validation/tests/test_blind_consumers.py @@ -13,7 +13,8 @@ @pytest.mark.parametrize("loader", ["load_pure_eb_results", "load_cosebis_results"]) @pytest.mark.parametrize("stamp", [None, "wrong"]) -def test_npz_consumers_refuse_missing_or_wrong_blind(tmp_path, loader, stamp): +@pytest.mark.parametrize("expected", ["none", "expected"]) +def test_npz_consumers_refuse_missing_or_wrong_blind(tmp_path, loader, stamp, expected): import numpy as np from sp_validation import b_modes @@ -24,7 +25,7 @@ def test_npz_consumers_refuse_missing_or_wrong_blind(tmp_path, loader, stamp): data["blind"] = stamp np.savez(path, **data) with pytest.raises(ValueError, match="stamp"): - getattr(b_modes, loader)(path, blind="expected") + getattr(b_modes, loader)(path, blind=expected) @pytest.fixture @@ -36,23 +37,6 @@ def part_path(tmp_path): return path -@pytest.mark.parametrize("expected", ["none", "toy"]) -def test_reporting_parts(part_path, monkeypatch, expected): - monkeypatch.syspath_prepend(str(SCRIPTS)) - from cv_runner import take_reporting_parts - - cv = SimpleNamespace( - versions=["v"], xi_parts={}, blind=lambda v: blinding.Blind(expected) - ) - if expected == "toy": - with pytest.raises(ValueError, match="stamped none.*toy"): - take_reporting_parts(cv, [part_path]) - assert not cv.xi_parts - else: - take_reporting_parts(cv, [part_path]) - assert sacc_io.stamp(cv.xi_parts["v", "reporting"]) == "none" - - @pytest.mark.parametrize( "script, input_key, reader", [ diff --git a/workflow/scripts/cv_runner.py b/workflow/scripts/cv_runner.py index 14430166..dda72778 100644 --- a/workflow/scripts/cv_runner.py +++ b/workflow/scripts/cv_runner.py @@ -50,15 +50,6 @@ def prepare_cached_psf(cv): ) -def take_reporting_parts(cv, paths): - """Hand ``cv`` the reporting ξ± parts at ``paths``, one per version, in order.""" - from sp_validation import sacc_io - - for version, path in zip(cv.versions, paths, strict=True): - part = sacc_io.load(path, blind=cv.blind(version)) - cv.xi_parts[version, "reporting"] = part - - def verify_outputs(snakemake): """Fail loudly if any declared output is missing after the method ran.""" missing = [str(o) for o in snakemake.output if not Path(o).exists()] From 4cde5961a7fedce2584f85a2c8277a3d0cecd1a5 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 24 Oct 2026 14:19:03 -0400 Subject: [PATCH 08/18] Keep notebook COSEBIs caches on their selected integration grid --- src/sp_validation/cosmo_val/cosebis.py | 9 ++-- src/sp_validation/tests/test_plot_products.py | 44 +++++++++++++++++++ 2 files changed, 49 insertions(+), 4 deletions(-) diff --git a/src/sp_validation/cosmo_val/cosebis.py b/src/sp_validation/cosmo_val/cosebis.py index b2714e6e..efce2cfc 100644 --- a/src/sp_validation/cosmo_val/cosebis.py +++ b/src/sp_validation/cosmo_val/cosebis.py @@ -297,13 +297,13 @@ def plot_cosebis( plot_results, version, var_method, bin_stub + "_covariance.png" ) + # Stored bin indices refer to the integration grid, including when + # the requested cuts were generated from a coarser reporting grid. + edges = log_bin_edges(min_sep_int, max_sep_int, nbins_int) if multiple_scale_cuts and len(bin_results) > 1: - binning = self._binning(min_sep, max_sep, nbins) plot_cosebis_scale_cut_heatmap( bin_results, - log_bin_edges( - binning["min_sep"], binning["max_sep"], binning["nbins"] - ), + edges, version, bin_stub + "_scalecut_ptes.png", fiducial_scale_cut=fiducial_scale_cut, @@ -313,6 +313,7 @@ def plot_cosebis( bin_results, bin_stub + "_data.npz", fiducial_scale_cut, + edges, derived_from=[self.xi_parts[version, "integration"]], ) self._cosebis_results[version] = results_tomo diff --git a/src/sp_validation/tests/test_plot_products.py b/src/sp_validation/tests/test_plot_products.py index 2b109c49..9c547c6e 100644 --- a/src/sp_validation/tests/test_plot_products.py +++ b/src/sp_validation/tests/test_plot_products.py @@ -301,6 +301,50 @@ def test_cosebis_heatmap_uses_stored_inward_bins(tmp_path, monkeypatch): assert ax.patches[0].get_xy() == (1, 1) +def test_notebook_cosebis_cache_and_heatmap_use_integration_edges( + tmp_path, monkeypatch +): + from sp_validation.cosmo_val import cosebis as module + + class CV(module.CosebisMixin): + versions = ["synthetic"] + cc = {"paths": {"output": str(tmp_path)}} + treecorr_config = {"npatch": 1} + xi_parts = { + ("synthetic", "integration"): SimpleNamespace(metadata={"blind": "toy"}) + } + _cosebis_results = {} + + def _binning(self, *args): + return dict(min_sep=1.0, max_sep=8.0, nbins=3) + + edges = b_modes.log_bin_edges(1.0, 8.0, 10) + results = { + cut: dict( + En=np.ones(2), + Bn=np.zeros(2), + cov=np.eye(4), + chi2_E=2.0, + chi2_B=0.0, + pte_B=1.0, + bin_indices=indices, + ) + for cut, indices in [((1.0, 8.0), [0, 10]), ((1.0, 4.0), [0, 6])] + } + monkeypatch.setattr(module, "plot_cosebis_modes", lambda *args, **kwargs: None) + monkeypatch.setattr(module, "plot_cosebis_covariance_matrix", lambda *args: None) + + def heatmap(data, used_edges, *args, **kwargs): + np.testing.assert_array_equal(used_edges, edges) + + monkeypatch.setattr(module, "plot_cosebis_scale_cut_heatmap", heatmap) + CV().plot_cosebis(results=results, min_sep_int=1.0, max_sep_int=8.0, nbins_int=10) + _, cached_edges = b_modes.load_cosebis_results( + next(tmp_path.glob("*.npz")), blind="toy" + ) + np.testing.assert_array_equal(cached_edges, edges) + + def test_cosebis_preserves_all_cuts_for_plotting(tmp_path): edges = np.array([1.0, 2.0, 4.0, 8.0]) results = { From 8aa855f5f1ee11317ae0cc423f9ac39e3209f5c8 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sun, 25 Oct 2026 07:16:54 -0400 Subject: [PATCH 09/18] fix(plots): compact cached COSEBIs cuts and unclip mode amplitudes --- src/sp_validation/b_modes.py | 57 ++++++++----------- src/sp_validation/tests/test_plot_products.py | 7 ++- 2 files changed, 30 insertions(+), 34 deletions(-) diff --git a/src/sp_validation/b_modes.py b/src/sp_validation/b_modes.py index 0e57210c..f9c3b0f8 100644 --- a/src/sp_validation/b_modes.py +++ b/src/sp_validation/b_modes.py @@ -920,16 +920,19 @@ def plot_cosebis_scale_cut_heatmap( (min_scale, max_scale) for cross-hatching """ left_edges, right_edges = edges - nbins = len(left_edges) + starts = sorted({int(r["bin_indices"][0]) for r in cosebis_results.values()}) + stops = sorted({int(r["bin_indices"][1]) for r in cosebis_results.values()}) + nx, ny = len(starts), len(stops) # Place only the cuts actually scanned, using their binning metadata. # A sparse scan must not be assigned to the first N triangular cells. - snr_matrix = np.full((nbins, nbins), np.nan) - pte_matrix = np.full((nbins, nbins), np.nan) + snr_matrix = np.full((nx, ny), np.nan) + pte_matrix = np.full((nx, ny), np.nan) for result in cosebis_results.values(): start, stop = result["bin_indices"] - snr_matrix[start, stop - 1] = np.sqrt(result["chi2_E"]) - pte_matrix[start, stop - 1] = result["pte_B"] + i, j = starts.index(start), stops.index(stop) + snr_matrix[i, j] = np.sqrt(result["chi2_E"]) + pte_matrix[i, j] = result["pte_B"] # Create 1x2 subplot layout fig, axs = plt.subplots(1, 2, figsize=(15, 7)) @@ -951,11 +954,11 @@ def plot_cosebis_scale_cut_heatmap( cmap=mako_cmap, vmin=0, vmax=snr_max, - extent=[0, nbins, 0, nbins], + extent=[0, nx, 0, ny], ) axs[0].set_title(r"E-mode SNR: $\sqrt{\chi^2_E}$") - axs[0].set_xlabel("Lower scale cut") - axs[0].set_ylabel("Upper scale cut") + axs[0].set_xlabel("Lower scale cut [arcmin]") + axs[0].set_ylabel("Upper scale cut [arcmin]") # Right plot: B-mode PTE with contours pte_plot_data = pte_matrix.T @@ -966,18 +969,10 @@ def plot_cosebis_scale_cut_heatmap( cmap=vlag_cmap, vmin=0, vmax=1, - extent=[0, nbins, 0, nbins], + extent=[0, nx, 0, ny], ) - # Add contours to PTE plot - cs = axs[1].contour( - pte_plot_data, - levels=[0.05, 0.95], - colors="black", - linewidths=1, - extent=[0, nbins, 0, nbins], - ) - axs[1].clabel(cs, inline=True, fontsize=8) + # Sparse discrete cuts do not support interpolated significance contours. axs[1].set_title("B-mode PTE") axs[1].set_xlabel("Lower scale cut") @@ -989,10 +984,10 @@ def plot_cosebis_scale_cut_heatmap( ) start_bin, stop_bin = cosebis_results[fiducial_key]["bin_indices"] - if stop_bin > start_bin and start_bin < nbins and stop_bin > 0: - # Add cross-hatching at the fiducial scale cut matrix element - rect_x = start_bin - rect_y = stop_bin - 1 + if start_bin in starts and stop_bin in stops: + # Mark the cached selected cut. + rect_x = starts.index(start_bin) + rect_y = stops.index(stop_bin) rect_width = 1 rect_height = 1 @@ -1012,15 +1007,10 @@ def plot_cosebis_scale_cut_heatmap( ) # Set angular scale ticks - tick_indices = np.arange(0, nbins) - x_tick_labels = [ - f"{left_edges[i]:.1f}" for i in tick_indices if i < len(left_edges) - ] - y_tick_labels = [ - f"{right_edges[i]:.1f}" for i in tick_indices if i < len(right_edges) - ] - x_tick_positions = tick_indices + 0.5 - y_tick_positions = tick_indices + 0.5 + x_tick_labels = [f"{left_edges[i]:.2f}" for i in starts] + y_tick_labels = [f"{right_edges[i - 1]:.2f}" for i in stops] + x_tick_positions = np.arange(nx) + 0.5 + y_tick_positions = np.arange(ny) + 0.5 for ax in axs: ax.set_xticks(x_tick_positions[: len(x_tick_labels)]) @@ -1270,7 +1260,10 @@ def plot_cosebis_modes(results, version, output_path, fiducial_scale_cut=None): plt.legend() plt.xlabel("n (mode)") plt.ylabel("E_n, B_n") - plt.ylim(-0.5e-10, 3e-10) + errors = np.sqrt(np.diag(results["cov"])) + values = np.concatenate([results["En"], results["Bn"]]) + limit = 1.1 * np.nanmax(np.abs(values) + errors) + plt.ylim(-limit, limit) # Add scale cut information to title - use actual scale cut from results scale_info = "" diff --git a/src/sp_validation/tests/test_plot_products.py b/src/sp_validation/tests/test_plot_products.py index 9c547c6e..f510b933 100644 --- a/src/sp_validation/tests/test_plot_products.py +++ b/src/sp_validation/tests/test_plot_products.py @@ -296,9 +296,12 @@ def test_cosebis_heatmap_uses_stored_inward_bins(tmp_path, monkeypatch): ) ax = plt.gcf().axes[1] matrix = ax.images[0].get_array() - assert matrix[1, 1] == pytest.approx(0.7) + assert matrix.shape == (1, 1) + assert matrix[0, 0] == pytest.approx(0.7) assert np.ma.count(matrix) == 1 - assert ax.patches[0].get_xy() == (1, 1) + assert ax.patches[0].get_xy() == (0, 0) + assert ax.get_xticklabels()[0].get_text() == "2.00" + assert ax.get_yticklabels()[0].get_text() == "4.00" def test_notebook_cosebis_cache_and_heatmap_use_integration_edges( From c9ade52ba15f8a0154d41b92032bc5d2e6a07847 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sun, 25 Oct 2026 07:18:30 -0400 Subject: [PATCH 10/18] fix(plots): derive symmetric pure E/B bounds including errors --- src/sp_validation/b_modes.py | 12 +++++++++++- src/sp_validation/tests/test_plot_products.py | 4 ++++ 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/src/sp_validation/b_modes.py b/src/sp_validation/b_modes.py index f9c3b0f8..5853b0c2 100644 --- a/src/sp_validation/b_modes.py +++ b/src/sp_validation/b_modes.py @@ -851,7 +851,17 @@ def plot_pure_eb_correlations( ax.set(xscale="log", xlabel=r"$\theta$ [arcmin]") ax.axhline(0, alpha=0.3, color="k", linestyle="--", linewidth=0.5) ax.legend(loc="upper left") - ax.set_ylim(-0.5, 2) + # Include every displayed mode and its error bars, symmetric about zero. + bounds = [] + for xi_type, var_attr in (("xip", "var_xip"), ("xim", "var_xim")): + bounds.extend(theta * (np.abs(results[xi_type]) + + np.sqrt(results[var_attr])) / scale_factor) + for mode in ("E", "B", "amb"): + key = f"{xi_type}_{mode}" + bounds.extend(theta * (np.abs(results[key]) + + results[f"std_{key}"]) / scale_factor) + limit = 1.1 * np.nanmax(bounds) + axs[0].set_ylim(-limit, limit) # Save the axis limits after data plotting but before adding gray regions original_xlims = [ax.get_xlim() for ax in axs] diff --git a/src/sp_validation/tests/test_plot_products.py b/src/sp_validation/tests/test_plot_products.py index f510b933..774e9e16 100644 --- a/src/sp_validation/tests/test_plot_products.py +++ b/src/sp_validation/tests/test_plot_products.py @@ -242,6 +242,10 @@ def test_pure_eb_plot_cache_is_complete(tmp_path): outputs = [tmp_path / f"eb{i}.png" for i in range(4)] b_modes.plot_integration_vs_reporting(loaded, outputs[0], "synthetic") b_modes.plot_pure_eb_correlations(loaded, outputs[1], "synthetic") + lo, hi = plt.gcf().axes[0].get_ylim() + assert lo == -hi + assert hi > np.max(theta * (loaded["xip"] + + np.sqrt(loaded["var_xip"])) / 1e-4) b_modes.plot_pte_2d_heatmaps(loaded, "synthetic", outputs[2]) b_modes.plot_eb_covariance_matrix( loaded["cov"], "analytic", outputs[3], "synthetic" From e50308779bc4db6c9b1e2d0dd1d055e970d3f028 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sun, 25 Oct 2026 07:19:07 -0400 Subject: [PATCH 11/18] fix(plots): label footprint galaxy counts and pixel NSIDE --- src/sp_validation/cosmo_val/diagnostic_plots.py | 5 ++++- src/sp_validation/tests/test_plot_products.py | 3 ++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/sp_validation/cosmo_val/diagnostic_plots.py b/src/sp_validation/cosmo_val/diagnostic_plots.py index 59ba45ef..fc16feb4 100644 --- a/src/sp_validation/cosmo_val/diagnostic_plots.py +++ b/src/sp_validation/cosmo_val/diagnostic_plots.py @@ -13,7 +13,10 @@ def plot_footprint(count_map, outputs): """Render every requested region from one cached HealSparse count map.""" plotter = FootprintPlotter() for region, path in outputs.items(): - plotter.plot_region(count_map, plotter._regions[region], outpath=str(path)) + plotter.plot_region( + count_map, plotter._regions[region], outpath=str(path), + colorbar_label=f"Galaxies per pixel (NSIDE={count_map.nside_sparse})", + ) plt.close() diff --git a/src/sp_validation/tests/test_plot_products.py b/src/sp_validation/tests/test_plot_products.py index 774e9e16..a406fa3a 100644 --- a/src/sp_validation/tests/test_plot_products.py +++ b/src/sp_validation/tests/test_plot_products.py @@ -143,7 +143,8 @@ def test_footprint_every_region_from_one_map(tmp_path, monkeypatch): loaded = load_footprint(cache) calls = [] - def draw(self, data, region, outpath): + def draw(self, data, region, outpath, colorbar_label): + assert colorbar_label == "Galaxies per pixel (NSIDE=16)" np.testing.assert_array_equal(data.valid_pixels, count_map.valid_pixels) calls.append(region) Path(outpath).write_bytes(b"rendered" * 100) From 49e0a92862ab640d117ec826ce4bb23420839df6 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sun, 25 Oct 2026 07:24:58 -0400 Subject: [PATCH 12/18] fix(plots): constrain xi panel layout to separate axis labels --- .../cosmo_val/diagnostic_plots.py | 2 +- .../tests/test_plotfix_layout.py | 20 +++++++++++++++++++ 2 files changed, 21 insertions(+), 1 deletion(-) create mode 100644 src/sp_validation/tests/test_plotfix_layout.py diff --git a/src/sp_validation/cosmo_val/diagnostic_plots.py b/src/sp_validation/cosmo_val/diagnostic_plots.py index fc16feb4..c00f140b 100644 --- a/src/sp_validation/cosmo_val/diagnostic_plots.py +++ b/src/sp_validation/cosmo_val/diagnostic_plots.py @@ -183,7 +183,7 @@ def plot_ratio_xi_sys_xi( def plot_2pcf(correlations, styles, outputs): """Draw sealed reporting xi±, both plain and theta-weighted.""" for times_theta, output in enumerate(outputs): - fig, axes = plt.subplots(1, 2, figsize=(14, 6)) + fig, axes = plt.subplots(1, 2, figsize=(14, 6), layout="constrained") for ver, gg in correlations.items(): factor = gg.meanr if times_theta else 1 for ax, xi, var, sign in zip( diff --git a/src/sp_validation/tests/test_plotfix_layout.py b/src/sp_validation/tests/test_plotfix_layout.py new file mode 100644 index 00000000..39da5ec8 --- /dev/null +++ b/src/sp_validation/tests/test_plotfix_layout.py @@ -0,0 +1,20 @@ +"""Presentation regressions with synthetic arrays only.""" +from types import SimpleNamespace + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np + +from sp_validation.cosmo_val import diagnostic_plots as plots + + +def test_xi_panels_use_constrained_layout(tmp_path, monkeypatch): + figures = [] + monkeypatch.setattr(plt, "close", lambda fig: figures.append(fig)) + gg = SimpleNamespace(meanr=np.array([1., 10.]), xip=np.array([1., -1.]), + xim=np.array([-1., 1.]), varxip=np.ones(2), + varxim=np.ones(2)) + plots.plot_2pcf({"toy": gg}, {"toy": {"marker": "o", "colour": "C0"}}, + [tmp_path / "xi.png", tmp_path / "theta_xi.png"]) + assert all(fig.get_constrained_layout() for fig in figures if hasattr(fig, "axes")) From d0a4973ddb039a48d1557c64036d2f74c6f43c48 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sun, 25 Oct 2026 07:31:25 -0400 Subject: [PATCH 13/18] fix(plots): retain signed PSF contamination and uncertainty bands --- .../cosmo_val/diagnostic_plots.py | 4 +-- src/sp_validation/tests/test_plotfix_psf.py | 26 +++++++++++++++++++ 2 files changed, 28 insertions(+), 2 deletions(-) create mode 100644 src/sp_validation/tests/test_plotfix_psf.py diff --git a/src/sp_validation/cosmo_val/diagnostic_plots.py b/src/sp_validation/cosmo_val/diagnostic_plots.py index c00f140b..e9a3185e 100644 --- a/src/sp_validation/cosmo_val/diagnostic_plots.py +++ b/src/sp_validation/cosmo_val/diagnostic_plots.py @@ -109,7 +109,7 @@ def plot_rho_tau_fits(products, styles, contours, samples, summary): close=True, ) for path, draw_samples in ((samples, True), (summary, False)): - fig, axes = plt.subplots(1, 2, figsize=(14, 6)) + fig, axes = plt.subplots(1, 2, figsize=(14, 6), layout="constrained") for ver, data in products.items(): theta, colour = data["theta"], styles[ver]["colour"] for ax, component, sign in zip(axes, ("plus", "minus"), ("+", "-")): @@ -128,7 +128,7 @@ def plot_rho_tau_fits(products, styles, contours, samples, summary): ax.fill_between(theta, lo, hi, color=colour, alpha=0.3) ax.set( xscale="log", - yscale="log", + yscale="linear", xlabel=r"$\theta$ [arcmin]", ylabel=rf"$\xi^{{\rm PSF,sys}}_{sign}(\theta)$", ) diff --git a/src/sp_validation/tests/test_plotfix_psf.py b/src/sp_validation/tests/test_plotfix_psf.py new file mode 100644 index 00000000..172c3980 --- /dev/null +++ b/src/sp_validation/tests/test_plotfix_psf.py @@ -0,0 +1,26 @@ +"""Signed PSF contamination must remain visible.""" +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np + +from sp_validation.cosmo_val import diagnostic_plots as plots + + +def test_signed_psf_summary(tmp_path, monkeypatch): + from shear_psf_leakage import plots as external + monkeypatch.setattr(external, "plot_contours", lambda *a, **kw: None) + figures = [] + monkeypatch.setattr(plt, "close", lambda fig: figures.append(fig)) + data = {"theta": np.array([1., 10.]), "flat_samples": np.ones((10, 3))} + for component in ("plus", "minus"): + data[f"samples_{component}"] = np.array([[1., -1.]]) + data[f"mean_{component}"] = np.array([1., -1.]) + data[f"quantiles_{component}"] = np.array([[0., -2.], [2., 0.]]) + plots.plot_rho_tau_fits({"toy": data}, {"toy": {"colour": "C0"}}, + tmp_path / "contours.png", tmp_path / "samples.png", + tmp_path / "summary.png") + for fig in figures: + if hasattr(fig, "axes"): + assert all(ax.get_yscale() == "linear" for ax in fig.axes) + assert all(ax.get_ylim()[0] < 0 for ax in fig.axes) From 81b5e8ed6092368a9002c3c38dda4c911ab58698 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sun, 25 Oct 2026 07:33:53 -0400 Subject: [PATCH 14/18] fix(plots): use signed symlog pseudo-Cl panel --- src/sp_validation/cosmo_val/pseudo_cl.py | 10 ++++++++-- src/sp_validation/tests/test_plotfix_cell.py | 19 +++++++++++++++++++ 2 files changed, 27 insertions(+), 2 deletions(-) create mode 100644 src/sp_validation/tests/test_plotfix_cell.py diff --git a/src/sp_validation/cosmo_val/pseudo_cl.py b/src/sp_validation/cosmo_val/pseudo_cl.py index 361369ac..bfd3b720 100644 --- a/src/sp_validation/cosmo_val/pseudo_cl.py +++ b/src/sp_validation/cosmo_val/pseudo_cl.py @@ -51,7 +51,7 @@ def plot_pseudo_cl_spectrum(datasets, spectrum, output_path): ``datasets`` maps a version to ``{"ell", "cl", "cov", "style"}``, where ``style`` carries the ``marker`` and ``colour`` the version is drawn with. """ - fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(8, 8)) + fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(8, 8), layout="constrained") minor_ticks = [i * 10 for i in range(1, 10)] + [i * 100 for i in range(1, 21)] for panel, scaled in ((ax[0], True), (ax[1], False)): @@ -77,7 +77,13 @@ def plot_pseudo_cl_spectrum(datasets, spectrum, output_path): panel.xaxis.set_ticks(minor_ticks, minor=True) ax[1].set_xlabel(r"$\ell$") - ax[1].set_yscale("log") + # A signed scale preserves negative spectra and both ends of error bars. + amplitudes = np.concatenate([np.abs(d["cl"]) for d in datasets.values()]) + positive = amplitudes[np.isfinite(amplitudes) & (amplitudes > 0)] + threshold = np.median(positive) if positive.size else 1.0 + ax[1].set_yscale("symlog", linthresh=threshold) + ax[1].yaxis.get_major_locator().set_params(numticks=7) + ax[1].axhline(0, color="grey", linewidth=0.5) plt.suptitle(f"Pseudo-Cl {spectrum} (Gaussian covariance)") plt.legend() plt.savefig(output_path) diff --git a/src/sp_validation/tests/test_plotfix_cell.py b/src/sp_validation/tests/test_plotfix_cell.py new file mode 100644 index 00000000..abba51b6 --- /dev/null +++ b/src/sp_validation/tests/test_plotfix_cell.py @@ -0,0 +1,19 @@ +"""Pseudo-Cl plots must preserve signs, including error bars.""" +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np + +from sp_validation.cosmo_val.pseudo_cl import plot_pseudo_cl_spectrum + + +def test_signed_cell(tmp_path, monkeypatch): + figures = [] + monkeypatch.setattr(plt, "close", lambda fig: figures.append(fig)) + plot_pseudo_cl_spectrum({"toy": {"ell": np.array([100., 400.]), + "cl": np.array([-1., 1.]), "cov": np.eye(2) * 4}}, + "EB", tmp_path / "eb.png") + fig = next(fig for fig in reversed(figures) if hasattr(fig, "axes")) + assert fig.axes[1].get_yscale() == "symlog" + lo, hi = fig.axes[1].get_ylim() + assert lo < -3 and hi > 3 From c9422bd7a4cb6e7b1c826cdbab09ae13e8af25fc Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sun, 25 Oct 2026 08:24:57 -0400 Subject: [PATCH 15/18] fix(plots): omit symlog labels crowded against zero --- src/sp_validation/cosmo_val/pseudo_cl.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/sp_validation/cosmo_val/pseudo_cl.py b/src/sp_validation/cosmo_val/pseudo_cl.py index bfd3b720..6ddc3507 100644 --- a/src/sp_validation/cosmo_val/pseudo_cl.py +++ b/src/sp_validation/cosmo_val/pseudo_cl.py @@ -83,6 +83,8 @@ def plot_pseudo_cl_spectrum(datasets, spectrum, output_path): threshold = np.median(positive) if positive.size else 1.0 ax[1].set_yscale("symlog", linthresh=threshold) ax[1].yaxis.get_major_locator().set_params(numticks=7) + ticks = ax[1].get_yticks() + ax[1].set_yticks([t for t in ticks if t == 0 or abs(t) >= 100 * threshold]) ax[1].axhline(0, color="grey", linewidth=0.5) plt.suptitle(f"Pseudo-Cl {spectrum} (Gaussian covariance)") plt.legend() From d2bb3d6ce7bb4fe8d4bbfdfbca3e718bcbc29161 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sun, 25 Oct 2026 08:35:44 -0400 Subject: [PATCH 16/18] fix(plots): separate rho and tau panel labels --- src/sp_validation/cosmo_val/diagnostic_plots.py | 1 + src/sp_validation/tests/test_plotfix_layout.py | 12 ++++++++++++ 2 files changed, 13 insertions(+) diff --git a/src/sp_validation/cosmo_val/diagnostic_plots.py b/src/sp_validation/cosmo_val/diagnostic_plots.py index e9a3185e..30d19bd9 100644 --- a/src/sp_validation/cosmo_val/diagnostic_plots.py +++ b/src/sp_validation/cosmo_val/diagnostic_plots.py @@ -216,6 +216,7 @@ def plot_psf_stats(tables, styles, output, kind, covariances=()): 3, figsize=(15, 10 if kind == "rho" else 6), squeeze=False, + layout="constrained", ) for i, stat in enumerate(indices): ax = axes.flat[i] diff --git a/src/sp_validation/tests/test_plotfix_layout.py b/src/sp_validation/tests/test_plotfix_layout.py index 39da5ec8..e2c875dc 100644 --- a/src/sp_validation/tests/test_plotfix_layout.py +++ b/src/sp_validation/tests/test_plotfix_layout.py @@ -18,3 +18,15 @@ def test_xi_panels_use_constrained_layout(tmp_path, monkeypatch): plots.plot_2pcf({"toy": gg}, {"toy": {"marker": "o", "colour": "C0"}}, [tmp_path / "xi.png", tmp_path / "theta_xi.png"]) assert all(fig.get_constrained_layout() for fig in figures if hasattr(fig, "axes")) + + +def test_tau_panels_use_constrained_layout(tmp_path, monkeypatch): + figures = [] + monkeypatch.setattr(plt, "close", lambda fig: figures.append(fig)) + table = {"theta": np.array([1., 10.])} + for i in (0, 2, 5): + table[f"tau_{i}_p"] = np.zeros(2) + table[f"vartau_{i}_p"] = np.ones(2) + plots.plot_psf_stats({"toy": table}, {"toy": {"colour": "C0"}}, + tmp_path / "tau.png", "tau") + assert all(fig.get_constrained_layout() for fig in figures if hasattr(fig, "axes")) From 9f6670961d55b41073cb24470cfd886be79c2857 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Mon, 26 Oct 2026 02:37:40 -0400 Subject: [PATCH 17/18] style: format plot renderers and presentation regression tests --- src/sp_validation/b_modes.py | 12 ++++++--- .../cosmo_val/diagnostic_plots.py | 4 ++- src/sp_validation/tests/test_plot_products.py | 3 +-- src/sp_validation/tests/test_plotfix_cell.py | 16 +++++++++--- .../tests/test_plotfix_layout.py | 26 +++++++++++++------ src/sp_validation/tests/test_plotfix_psf.py | 21 ++++++++++----- 6 files changed, 57 insertions(+), 25 deletions(-) diff --git a/src/sp_validation/b_modes.py b/src/sp_validation/b_modes.py index 5853b0c2..48349043 100644 --- a/src/sp_validation/b_modes.py +++ b/src/sp_validation/b_modes.py @@ -854,12 +854,16 @@ def plot_pure_eb_correlations( # Include every displayed mode and its error bars, symmetric about zero. bounds = [] for xi_type, var_attr in (("xip", "var_xip"), ("xim", "var_xim")): - bounds.extend(theta * (np.abs(results[xi_type]) - + np.sqrt(results[var_attr])) / scale_factor) + bounds.extend( + theta + * (np.abs(results[xi_type]) + np.sqrt(results[var_attr])) + / scale_factor + ) for mode in ("E", "B", "amb"): key = f"{xi_type}_{mode}" - bounds.extend(theta * (np.abs(results[key]) - + results[f"std_{key}"]) / scale_factor) + bounds.extend( + theta * (np.abs(results[key]) + results[f"std_{key}"]) / scale_factor + ) limit = 1.1 * np.nanmax(bounds) axs[0].set_ylim(-limit, limit) diff --git a/src/sp_validation/cosmo_val/diagnostic_plots.py b/src/sp_validation/cosmo_val/diagnostic_plots.py index 30d19bd9..6634c03a 100644 --- a/src/sp_validation/cosmo_val/diagnostic_plots.py +++ b/src/sp_validation/cosmo_val/diagnostic_plots.py @@ -14,7 +14,9 @@ def plot_footprint(count_map, outputs): plotter = FootprintPlotter() for region, path in outputs.items(): plotter.plot_region( - count_map, plotter._regions[region], outpath=str(path), + count_map, + plotter._regions[region], + outpath=str(path), colorbar_label=f"Galaxies per pixel (NSIDE={count_map.nside_sparse})", ) plt.close() diff --git a/src/sp_validation/tests/test_plot_products.py b/src/sp_validation/tests/test_plot_products.py index a406fa3a..2b72f2e0 100644 --- a/src/sp_validation/tests/test_plot_products.py +++ b/src/sp_validation/tests/test_plot_products.py @@ -245,8 +245,7 @@ def test_pure_eb_plot_cache_is_complete(tmp_path): b_modes.plot_pure_eb_correlations(loaded, outputs[1], "synthetic") lo, hi = plt.gcf().axes[0].get_ylim() assert lo == -hi - assert hi > np.max(theta * (loaded["xip"] - + np.sqrt(loaded["var_xip"])) / 1e-4) + assert hi > np.max(theta * (loaded["xip"] + np.sqrt(loaded["var_xip"])) / 1e-4) b_modes.plot_pte_2d_heatmaps(loaded, "synthetic", outputs[2]) b_modes.plot_eb_covariance_matrix( loaded["cov"], "analytic", outputs[3], "synthetic" diff --git a/src/sp_validation/tests/test_plotfix_cell.py b/src/sp_validation/tests/test_plotfix_cell.py index abba51b6..e652d159 100644 --- a/src/sp_validation/tests/test_plotfix_cell.py +++ b/src/sp_validation/tests/test_plotfix_cell.py @@ -1,5 +1,7 @@ """Pseudo-Cl plots must preserve signs, including error bars.""" + import matplotlib + matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np @@ -10,9 +12,17 @@ def test_signed_cell(tmp_path, monkeypatch): figures = [] monkeypatch.setattr(plt, "close", lambda fig: figures.append(fig)) - plot_pseudo_cl_spectrum({"toy": {"ell": np.array([100., 400.]), - "cl": np.array([-1., 1.]), "cov": np.eye(2) * 4}}, - "EB", tmp_path / "eb.png") + plot_pseudo_cl_spectrum( + { + "toy": { + "ell": np.array([100.0, 400.0]), + "cl": np.array([-1.0, 1.0]), + "cov": np.eye(2) * 4, + } + }, + "EB", + tmp_path / "eb.png", + ) fig = next(fig for fig in reversed(figures) if hasattr(fig, "axes")) assert fig.axes[1].get_yscale() == "symlog" lo, hi = fig.axes[1].get_ylim() diff --git a/src/sp_validation/tests/test_plotfix_layout.py b/src/sp_validation/tests/test_plotfix_layout.py index e2c875dc..2610b5e3 100644 --- a/src/sp_validation/tests/test_plotfix_layout.py +++ b/src/sp_validation/tests/test_plotfix_layout.py @@ -1,7 +1,9 @@ """Presentation regressions with synthetic arrays only.""" + from types import SimpleNamespace import matplotlib + matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np @@ -12,21 +14,29 @@ def test_xi_panels_use_constrained_layout(tmp_path, monkeypatch): figures = [] monkeypatch.setattr(plt, "close", lambda fig: figures.append(fig)) - gg = SimpleNamespace(meanr=np.array([1., 10.]), xip=np.array([1., -1.]), - xim=np.array([-1., 1.]), varxip=np.ones(2), - varxim=np.ones(2)) - plots.plot_2pcf({"toy": gg}, {"toy": {"marker": "o", "colour": "C0"}}, - [tmp_path / "xi.png", tmp_path / "theta_xi.png"]) + gg = SimpleNamespace( + meanr=np.array([1.0, 10.0]), + xip=np.array([1.0, -1.0]), + xim=np.array([-1.0, 1.0]), + varxip=np.ones(2), + varxim=np.ones(2), + ) + plots.plot_2pcf( + {"toy": gg}, + {"toy": {"marker": "o", "colour": "C0"}}, + [tmp_path / "xi.png", tmp_path / "theta_xi.png"], + ) assert all(fig.get_constrained_layout() for fig in figures if hasattr(fig, "axes")) def test_tau_panels_use_constrained_layout(tmp_path, monkeypatch): figures = [] monkeypatch.setattr(plt, "close", lambda fig: figures.append(fig)) - table = {"theta": np.array([1., 10.])} + table = {"theta": np.array([1.0, 10.0])} for i in (0, 2, 5): table[f"tau_{i}_p"] = np.zeros(2) table[f"vartau_{i}_p"] = np.ones(2) - plots.plot_psf_stats({"toy": table}, {"toy": {"colour": "C0"}}, - tmp_path / "tau.png", "tau") + plots.plot_psf_stats( + {"toy": table}, {"toy": {"colour": "C0"}}, tmp_path / "tau.png", "tau" + ) assert all(fig.get_constrained_layout() for fig in figures if hasattr(fig, "axes")) diff --git a/src/sp_validation/tests/test_plotfix_psf.py b/src/sp_validation/tests/test_plotfix_psf.py index 172c3980..5b6c3359 100644 --- a/src/sp_validation/tests/test_plotfix_psf.py +++ b/src/sp_validation/tests/test_plotfix_psf.py @@ -1,5 +1,7 @@ """Signed PSF contamination must remain visible.""" + import matplotlib + matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np @@ -9,17 +11,22 @@ def test_signed_psf_summary(tmp_path, monkeypatch): from shear_psf_leakage import plots as external + monkeypatch.setattr(external, "plot_contours", lambda *a, **kw: None) figures = [] monkeypatch.setattr(plt, "close", lambda fig: figures.append(fig)) - data = {"theta": np.array([1., 10.]), "flat_samples": np.ones((10, 3))} + data = {"theta": np.array([1.0, 10.0]), "flat_samples": np.ones((10, 3))} for component in ("plus", "minus"): - data[f"samples_{component}"] = np.array([[1., -1.]]) - data[f"mean_{component}"] = np.array([1., -1.]) - data[f"quantiles_{component}"] = np.array([[0., -2.], [2., 0.]]) - plots.plot_rho_tau_fits({"toy": data}, {"toy": {"colour": "C0"}}, - tmp_path / "contours.png", tmp_path / "samples.png", - tmp_path / "summary.png") + data[f"samples_{component}"] = np.array([[1.0, -1.0]]) + data[f"mean_{component}"] = np.array([1.0, -1.0]) + data[f"quantiles_{component}"] = np.array([[0.0, -2.0], [2.0, 0.0]]) + plots.plot_rho_tau_fits( + {"toy": data}, + {"toy": {"colour": "C0"}}, + tmp_path / "contours.png", + tmp_path / "samples.png", + tmp_path / "summary.png", + ) for fig in figures: if hasattr(fig, "axes"): assert all(ax.get_yscale() == "linear" for ax in fig.axes) From e4805219585f2ea6adaa7f75493a723ce39aae3a Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 28 Oct 2026 19:59:43 +0100 Subject: [PATCH 18/18] fix: retain default rho tau covariance method Co-Authored-By: Claude Opus 5.5 --- workflow/scripts/run_rho_tau.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/workflow/scripts/run_rho_tau.py b/workflow/scripts/run_rho_tau.py index 21712470..995fbdc1 100644 --- a/workflow/scripts/run_rho_tau.py +++ b/workflow/scripts/run_rho_tau.py @@ -18,7 +18,7 @@ npatch=int(params["npatch"]), catalog_config=params["cat_config"], output_dir=params["output_dir"], - cov_estimate_method=params["cov_estimate_method"], + cov_estimate_method=params.get("cov_estimate_method", "th"), ) cv.calculate_rho_tau_stats(tomography=False) verify_outputs(snakemake)