From 148b699c6f49dc705b09d8c7256f2a86810b01b7 Mon Sep 17 00:00:00 2001 From: Alastair Hewitt Date: Tue, 21 Oct 2025 00:35:47 -0400 Subject: [PATCH 1/9] Handle binning mismatches in cross-spectrum script --- scripts/compute_cross_spectrum.py | 37 +++++++++++++++++++++++++++---- 1 file changed, 33 insertions(+), 4 deletions(-) diff --git a/scripts/compute_cross_spectrum.py b/scripts/compute_cross_spectrum.py index 1cfee4f..0213134 100644 --- a/scripts/compute_cross_spectrum.py +++ b/scripts/compute_cross_spectrum.py @@ -27,7 +27,7 @@ def _load_bin_info( fallback_nlb: int | None, ) -> _BinInfo: try: - bins, _meta = load_bins_from_prereg(prereg, nside=nside) + bins, bins_meta = load_bins_from_prereg(prereg, nside=nside) except FileNotFoundError: bins = None except Exception as exc: # pragma: no cover - defensive fallback for CLI use @@ -43,12 +43,23 @@ def _load_bin_info( ell_lists.append(np.asarray(group, dtype=int)) except Exception: # pragma: no cover - guard against API mismatches ell_lists = [] - return _BinInfo(ell_effective=ell_eff, ell_lists=ell_lists) + if ell_eff.size == nbins: + return _BinInfo(ell_effective=ell_eff, ell_lists=ell_lists) + + summary_line( + "prereg bin definition does not match input spectra; " + "falling back to CLI-provided binning" + ) + if fallback_lmin is None and isinstance(bins_meta, dict): + fallback_lmin = bins_meta.get("lmin") + if fallback_nlb is None and isinstance(bins_meta, dict): + fallback_nlb = bins_meta.get("nlb") if fallback_nlb is None: raise ValueError("Prereg bins unavailable; please provide --nlb for fallback binning") if fallback_lmin is None: - raise ValueError("Prereg bins unavailable; please provide --lmin for fallback binning") + summary_line("no fallback lmin supplied; assuming lmin=0 for binning") + fallback_lmin = 0 lmin = int(fallback_lmin) step = int(fallback_nlb) @@ -137,12 +148,30 @@ def main(argv: Iterable[str] | None = None) -> int: nbins = cl_a.size nside = int(order_a.get("nside", order_b.get("nside", 0))) + def _extract_scalar(value: np.ndarray | None) -> int | None: + if value is None: + return None + try: + arr = np.asarray(value) + except Exception: + return None + if arr.size == 0: + return None + try: + return int(arr.ravel()[0]) + except Exception: + return None + + fallback_nlb = args.nlb + if fallback_nlb is None: + fallback_nlb = _extract_scalar(order_a.get("nlb")) or _extract_scalar(order_b.get("nlb")) + info = _load_bin_info( args.prereg, nside=nside if nside > 0 else 256, nbins=nbins, fallback_lmin=args.lmin, - fallback_nlb=args.nlb, + fallback_nlb=fallback_nlb, ) ell_eff = info.ell_effective From 0da51836689a56ac6b350d4b906618f0685e8c17 Mon Sep 17 00:00:00 2001 From: Alastair Hewitt Date: Tue, 21 Oct 2025 00:47:56 -0400 Subject: [PATCH 2/9] Require explicit fallback binning parameters --- README.md | 8 ++++++-- scripts/compute_cross_spectrum.py | 7 +------ tests/test_compute_cross_spectrum.py | 29 ++++++++++++++++++++++++++++ 3 files changed, 36 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index 017b247..2eaf7aa 100644 --- a/README.md +++ b/README.md @@ -216,12 +216,16 @@ earlier. --order-a artifacts/order_A_to_B_full.npz \ --order-b artifacts/order_B_to_A_full.npz \ --theory data/theory/tk_planck2018.npz \ + --lmin 30 \ + --nlb 30 \ --cov artifacts/cov_delta_full.npy \ --out artifacts/cross_tk_full.npz \ --summary artifacts/cross_summary_full.json ``` - Check `artifacts/cross_summary_full.json` to confirm the mean and - maximum |z| are consistent with a null detection. + The explicit `--lmin`/`--nlb` values ensure the theory is binned with + the same geometry as the commutator runs when the preregistration + metadata diverges. Check `artifacts/cross_summary_full.json` to confirm + the mean and maximum |z| are consistent with a null detection. 7. **Generate publication figures and a textual digest.** ```bash diff --git a/scripts/compute_cross_spectrum.py b/scripts/compute_cross_spectrum.py index 0213134..1dfef82 100644 --- a/scripts/compute_cross_spectrum.py +++ b/scripts/compute_cross_spectrum.py @@ -50,16 +50,11 @@ def _load_bin_info( "prereg bin definition does not match input spectra; " "falling back to CLI-provided binning" ) - if fallback_lmin is None and isinstance(bins_meta, dict): - fallback_lmin = bins_meta.get("lmin") - if fallback_nlb is None and isinstance(bins_meta, dict): - fallback_nlb = bins_meta.get("nlb") if fallback_nlb is None: raise ValueError("Prereg bins unavailable; please provide --nlb for fallback binning") if fallback_lmin is None: - summary_line("no fallback lmin supplied; assuming lmin=0 for binning") - fallback_lmin = 0 + raise ValueError("Prereg bins unavailable; please provide --lmin for fallback binning") lmin = int(fallback_lmin) step = int(fallback_nlb) diff --git a/tests/test_compute_cross_spectrum.py b/tests/test_compute_cross_spectrum.py index e9a7523..fa1bf76 100644 --- a/tests/test_compute_cross_spectrum.py +++ b/tests/test_compute_cross_spectrum.py @@ -85,3 +85,32 @@ def test_compute_cross_spectrum_outputs(tmp_path, monkeypatch): summary_payload = json.loads(summary.read_text()) assert summary_payload["nbins"] == 3 assert messages and "cross-spectrum" in messages[-1] + + +def test_compute_cross_spectrum_requires_lmin_when_falling_back(tmp_path): + module = _load_script_module( + "compute_cross_spectrum", Path("scripts/compute_cross_spectrum.py") + ) + + cl = np.array([1.0, 2.0, 3.0], dtype=float) + order = tmp_path / "order.npz" + np.savez(order, cl=cl, nside=256, nlb=10) + + theory_path = tmp_path / "theory.npz" + np.savez(theory_path, ell=np.arange(0, 40, dtype=float), cl_tk=0.05 * np.arange(0, 40)) + + with pytest.raises(ValueError, match="provide --lmin"): + module.main( + [ + "--order-a", + str(order), + "--order-b", + str(order), + "--theory", + str(theory_path), + "--prereg", + str(tmp_path / "missing_prereg.yaml"), + "--nlb", + "10", + ] + ) From 258198816c2de24b73bd0e824c29aa76fa9c3986 Mon Sep 17 00:00:00 2001 From: Alastair Hewitt Date: Tue, 21 Oct 2025 00:59:07 -0400 Subject: [PATCH 3/9] Embed bin metadata in order outputs --- README.md | 6 +++ scripts/compute_cross_spectrum.py | 79 ++++++++++++++++++++++------ scripts/run_order_A_to_B.py | 47 ++++++++++++++++- scripts/run_order_B_to_A.py | 47 ++++++++++++++++- tests/test_compute_cross_spectrum.py | 44 +++++++++++++--- 5 files changed, 197 insertions(+), 26 deletions(-) diff --git a/README.md b/README.md index 2eaf7aa..2dd2ea3 100644 --- a/README.md +++ b/README.md @@ -227,6 +227,12 @@ earlier. metadata diverges. Check `artifacts/cross_summary_full.json` to confirm the mean and maximum |z| are consistent with a null detection. + The ordering scripts now record the binning limits alongside the + spectra, so the cross-spectrum CLI can infer `lmin`/`nlb` directly from + those artifacts. Passing the explicit values remains recommended for a + reproducible command log, especially when exchanging files generated by + older commits that predate the embedded metadata. + 7. **Generate publication figures and a textual digest.** ```bash micromamba run -n comet python scripts/summarize_results.py \ diff --git a/scripts/compute_cross_spectrum.py b/scripts/compute_cross_spectrum.py index 1dfef82..29d410a 100644 --- a/scripts/compute_cross_spectrum.py +++ b/scripts/compute_cross_spectrum.py @@ -2,7 +2,7 @@ import argparse import json -from collections.abc import Iterable, Sequence +from collections.abc import Iterable, Mapping, Sequence from dataclasses import dataclass from pathlib import Path @@ -115,6 +115,49 @@ def _bin_theory( return np.zeros(0, dtype=float) +def _extract_scalar(value: object) -> int | None: + if value is None: + return None + if isinstance(value, (int, np.integer)): + return int(value) + if isinstance(value, (float, np.floating)): + if np.isnan(value): + return None + return int(value) + try: + arr = np.asarray(value) + except Exception: + return None + if arr.size == 0: + return None + try: + return int(arr.reshape(-1)[0]) + except Exception: + return None + + +def _load_order_metadata(order_path: Path) -> Mapping[str, object]: + meta_path = order_path.with_suffix(".json") + if not meta_path.exists(): + return {} + try: + return json.loads(meta_path.read_text()) + except Exception as exc: # pragma: no cover - metadata is optional for CLI use + summary_line(f"failed to parse {meta_path.name}: {exc}; ignoring metadata") + return {} + + +def _extract_from_metadata(meta: Mapping[str, object], key: str) -> object | None: + if not isinstance(meta, Mapping): + return None + if key in meta: + return meta[key] + bins = meta.get("bins") + if isinstance(bins, Mapping) and key in bins: + return bins[key] + return None + + def main(argv: Iterable[str] | None = None) -> int: ap = argparse.ArgumentParser(description="Build the T×κ cross-spectrum science product") ap.add_argument("--order-a", default="artifacts/order_A_to_B.npz") @@ -143,29 +186,33 @@ def main(argv: Iterable[str] | None = None) -> int: nbins = cl_a.size nside = int(order_a.get("nside", order_b.get("nside", 0))) - def _extract_scalar(value: np.ndarray | None) -> int | None: - if value is None: - return None - try: - arr = np.asarray(value) - except Exception: - return None - if arr.size == 0: - return None - try: - return int(arr.ravel()[0]) - except Exception: - return None - fallback_nlb = args.nlb if fallback_nlb is None: fallback_nlb = _extract_scalar(order_a.get("nlb")) or _extract_scalar(order_b.get("nlb")) + meta_a: Mapping[str, object] | None = None + meta_b: Mapping[str, object] | None = None + if fallback_nlb is None or args.lmin is None: + meta_a = _load_order_metadata(Path(args.order_a)) + meta_b = _load_order_metadata(Path(args.order_b)) + + if fallback_nlb is None: + fallback_nlb = _extract_scalar(_extract_from_metadata(meta_a or {}, "nlb")) or _extract_scalar( + _extract_from_metadata(meta_b or {}, "nlb") + ) + + fallback_lmin = args.lmin + if fallback_lmin is None: + fallback_lmin = _extract_scalar(order_a.get("lmin")) or _extract_scalar(order_b.get("lmin")) + if fallback_lmin is None and meta_a is not None and meta_b is not None: + fallback_lmin = _extract_scalar(_extract_from_metadata(meta_a, "lmin")) or _extract_scalar( + _extract_from_metadata(meta_b, "lmin") + ) info = _load_bin_info( args.prereg, nside=nside if nside > 0 else 256, nbins=nbins, - fallback_lmin=args.lmin, + fallback_lmin=fallback_lmin, fallback_nlb=fallback_nlb, ) diff --git a/scripts/run_order_A_to_B.py b/scripts/run_order_A_to_B.py index 85446e0..a70805e 100755 --- a/scripts/run_order_A_to_B.py +++ b/scripts/run_order_A_to_B.py @@ -2,6 +2,7 @@ from __future__ import annotations import argparse +from collections.abc import Mapping from pathlib import Path import healpy as hp @@ -85,6 +86,30 @@ def main(): nlb=args.nlb, ) + def _infer_bin_limits() -> tuple[int | None, int | None]: + lmin_val = getattr(bins, "lmin", None) + lmax_val = getattr(bins, "lmax", None) + ell_list = getattr(bins, "get_ell_list", None) + if callable(ell_list): + try: + groups = list(ell_list()) + except Exception: + groups = [] + for group in groups: + arr = np.asarray(group, dtype=int) + if arr.size: + lmin_val = int(arr[0]) + break + for group in reversed(groups): + arr = np.asarray(group, dtype=int) + if arr.size: + lmax_val = int(arr[-1]) + break + return ( + None if lmin_val is None else int(lmin_val), + None if lmax_val is None else int(lmax_val), + ) + if not args.disable_prereg: try: windows_cfg = load_windows_from_prereg(args.prereg) @@ -113,7 +138,25 @@ def main(): field_names=("cmb", "phi"), ) - np.savez(Path(args.out), cl=cl, nside=nside, nlb=args.nlb) + lmin_used, lmax_used = _infer_bin_limits() + payload: dict[str, object] = {"cl": cl, "nside": nside, "nlb": int(args.nlb)} + if lmin_used is not None: + payload["lmin"] = np.array(lmin_used, dtype=int) + if lmax_used is not None: + payload["lmax"] = np.array(lmax_used, dtype=int) + np.savez(Path(args.out), **payload) + + bins_summary: dict[str, object] = {} + if isinstance(bins_meta, Mapping): + bins_summary.update({k: v for k, v in bins_meta.items()}) + if lmin_used is not None and "lmin" not in bins_summary: + bins_summary["lmin"] = int(lmin_used) + if lmax_used is not None and "lmax" not in bins_summary: + bins_summary["lmax"] = int(lmax_used) + if "nlb" not in bins_summary: + bins_summary["nlb"] = int(args.nlb) + bins_summary["nbins"] = int(cl.size) + save_json( { "order": "A_to_B", @@ -122,7 +165,7 @@ def main(): "mask_threshold_sigma": args.threshold_sigma, "mask_apod_arcmin": args.apod_arcmin, "bins_source": "prereg" if bins_meta is not None else "cli", - "bins": bins_meta, + "bins": bins_summary, "windows": windows_cfg.to_metadata() if windows_cfg is not None else None, }, Path(args.out).with_suffix(".json"), diff --git a/scripts/run_order_B_to_A.py b/scripts/run_order_B_to_A.py index 0ca81e6..54c9cef 100755 --- a/scripts/run_order_B_to_A.py +++ b/scripts/run_order_B_to_A.py @@ -2,6 +2,7 @@ from __future__ import annotations import argparse +from collections.abc import Mapping from pathlib import Path import healpy as hp @@ -85,6 +86,30 @@ def main(): nlb=args.nlb, ) + def _infer_bin_limits() -> tuple[int | None, int | None]: + lmin_val = getattr(bins, "lmin", None) + lmax_val = getattr(bins, "lmax", None) + ell_list = getattr(bins, "get_ell_list", None) + if callable(ell_list): + try: + groups = list(ell_list()) + except Exception: + groups = [] + for group in groups: + arr = np.asarray(group, dtype=int) + if arr.size: + lmin_val = int(arr[0]) + break + for group in reversed(groups): + arr = np.asarray(group, dtype=int) + if arr.size: + lmax_val = int(arr[-1]) + break + return ( + None if lmin_val is None else int(lmin_val), + None if lmax_val is None else int(lmax_val), + ) + if not args.disable_prereg: try: windows_cfg = load_windows_from_prereg(args.prereg) @@ -113,7 +138,25 @@ def main(): field_names=("phi", "cmb"), ) - np.savez(Path(args.out), cl=cl, nside=nside, nlb=args.nlb) + lmin_used, lmax_used = _infer_bin_limits() + payload: dict[str, object] = {"cl": cl, "nside": nside, "nlb": int(args.nlb)} + if lmin_used is not None: + payload["lmin"] = np.array(lmin_used, dtype=int) + if lmax_used is not None: + payload["lmax"] = np.array(lmax_used, dtype=int) + np.savez(Path(args.out), **payload) + + bins_summary: dict[str, object] = {} + if isinstance(bins_meta, Mapping): + bins_summary.update({k: v for k, v in bins_meta.items()}) + if lmin_used is not None and "lmin" not in bins_summary: + bins_summary["lmin"] = int(lmin_used) + if lmax_used is not None and "lmax" not in bins_summary: + bins_summary["lmax"] = int(lmax_used) + if "nlb" not in bins_summary: + bins_summary["nlb"] = int(args.nlb) + bins_summary["nbins"] = int(cl.size) + save_json( { "order": "B_to_A", @@ -122,7 +165,7 @@ def main(): "mask_threshold_sigma": args.threshold_sigma, "mask_apod_arcmin": args.apod_arcmin, "bins_source": "prereg" if bins_meta is not None else "cli", - "bins": bins_meta, + "bins": bins_summary, "windows": windows_cfg.to_metadata() if windows_cfg is not None else None, }, Path(args.out).with_suffix(".json"), diff --git a/tests/test_compute_cross_spectrum.py b/tests/test_compute_cross_spectrum.py index fa1bf76..fbb7be9 100644 --- a/tests/test_compute_cross_spectrum.py +++ b/tests/test_compute_cross_spectrum.py @@ -28,8 +28,8 @@ def test_compute_cross_spectrum_outputs(tmp_path, monkeypatch): cl_b = np.array([1.2, 1.8, 3.2], dtype=float) order_a = tmp_path / "order_a.npz" order_b = tmp_path / "order_b.npz" - np.savez(order_a, cl=cl_a, nside=256) - np.savez(order_b, cl=cl_b, nside=256) + np.savez(order_a, cl=cl_a, nside=256, nlb=5, lmin=10) + np.savez(order_b, cl=cl_b, nside=256, nlb=5, lmin=10) diag = np.array([0.25, 0.16, 0.09], dtype=float) cov = tmp_path / "cov.npy" @@ -60,10 +60,6 @@ def test_compute_cross_spectrum_outputs(tmp_path, monkeypatch): str(theory_path), "--prereg", str(tmp_path / "missing_prereg.yaml"), - "--lmin", - "10", - "--nlb", - "5", "--cov", str(cov), "--out", @@ -114,3 +110,39 @@ def test_compute_cross_spectrum_requires_lmin_when_falling_back(tmp_path): "10", ] ) + + +def test_compute_cross_spectrum_reads_lmin_from_json_metadata(tmp_path): + module = _load_script_module( + "compute_cross_spectrum", Path("scripts/compute_cross_spectrum.py") + ) + + cl = np.array([1.0, 2.0, 3.0], dtype=float) + order = tmp_path / "order.npz" + np.savez(order, cl=cl, nside=256, nlb=5) + (order.with_suffix(".json")).write_text(json.dumps({"bins": {"lmin": 0, "nlb": 5}})) + + theory_path = tmp_path / "theory.npz" + np.savez(theory_path, ell=np.arange(0, 40, dtype=float), cl_tk=0.05 * np.arange(0, 40)) + + out = tmp_path / "cross.npz" + summary = tmp_path / "summary.json" + module.main( + [ + "--order-a", + str(order), + "--order-b", + str(order), + "--theory", + str(theory_path), + "--prereg", + str(tmp_path / "missing_prereg.yaml"), + "--out", + str(out), + "--summary", + str(summary), + ] + ) + + data = np.load(out) + np.testing.assert_allclose(data["ell"], np.array([2.0, 7.0, 12.0])) From d0d9c46271f1feb3af6368ca57792d763cdc2c27 Mon Sep 17 00:00:00 2001 From: Alastair Hewitt Date: Tue, 21 Oct 2025 18:58:28 -0400 Subject: [PATCH 4/9] Persist default lmin in order artifacts --- scripts/run_order_A_to_B.py | 14 ++++++++++++++ scripts/run_order_B_to_A.py | 13 +++++++++++++ 2 files changed, 27 insertions(+) diff --git a/scripts/run_order_A_to_B.py b/scripts/run_order_A_to_B.py index a70805e..5332a9c 100755 --- a/scripts/run_order_A_to_B.py +++ b/scripts/run_order_A_to_B.py @@ -139,6 +139,20 @@ def _infer_bin_limits() -> tuple[int | None, int | None]: ) lmin_used, lmax_used = _infer_bin_limits() + if lmin_used is None: + if args.lmin is not None: + lmin_used = int(args.lmin) + elif isinstance(bins_meta, Mapping) and "lmin" in bins_meta: + try: + lmin_used = int(bins_meta["lmin"]) + except Exception: + lmin_used = None + if lmin_used is None: + # NaMaster's linear bin helper defaults to starting at ℓ=0 when no + # explicit lower bound is supplied. Persist that implicit origin so + # downstream stages can always reconstruct the band edges when + # preregistered metadata is unavailable. + lmin_used = 0 payload: dict[str, object] = {"cl": cl, "nside": nside, "nlb": int(args.nlb)} if lmin_used is not None: payload["lmin"] = np.array(lmin_used, dtype=int) diff --git a/scripts/run_order_B_to_A.py b/scripts/run_order_B_to_A.py index 54c9cef..a6a3d02 100755 --- a/scripts/run_order_B_to_A.py +++ b/scripts/run_order_B_to_A.py @@ -139,6 +139,19 @@ def _infer_bin_limits() -> tuple[int | None, int | None]: ) lmin_used, lmax_used = _infer_bin_limits() + if lmin_used is None: + if args.lmin is not None: + lmin_used = int(args.lmin) + elif isinstance(bins_meta, Mapping) and "lmin" in bins_meta: + try: + lmin_used = int(bins_meta["lmin"]) + except Exception: + lmin_used = None + if lmin_used is None: + # Mirror the implicit NaMaster behaviour of starting at ℓ=0 so the + # cross-spectrum step can fall back to linear binning even when + # preregistration metadata is absent. + lmin_used = 0 payload: dict[str, object] = {"cl": cl, "nside": nside, "nlb": int(args.nlb)} if lmin_used is not None: payload["lmin"] = np.array(lmin_used, dtype=int) From 7170719e73faea1a5676ed33e9ae113abe36e40f Mon Sep 17 00:00:00 2001 From: Alastair Hewitt Date: Tue, 21 Oct 2025 19:36:00 -0400 Subject: [PATCH 5/9] Record bin geometry in order outputs --- scripts/compute_cross_spectrum.py | 110 ++++++++++++++++++++------- scripts/run_order_A_to_B.py | 74 +++++++++++++----- scripts/run_order_B_to_A.py | 74 +++++++++++++----- tests/test_compute_cross_spectrum.py | 67 +++++++++++++++- 4 files changed, 261 insertions(+), 64 deletions(-) diff --git a/scripts/compute_cross_spectrum.py b/scripts/compute_cross_spectrum.py index 29d410a..dae934c 100644 --- a/scripts/compute_cross_spectrum.py +++ b/scripts/compute_cross_spectrum.py @@ -18,6 +18,58 @@ class _BinInfo: ell_lists: list[np.ndarray] +def _bin_info_from_orders( + orders: Sequence[np.lib.npyio.NpzFile], + nbins: int, +) -> _BinInfo | None: + for order in orders: + ell_eff_raw = order.get("ell_effective") + ell_left = order.get("ell_left_edges") + ell_right = order.get("ell_right_edges") + + ell_lists: list[np.ndarray] = [] + lists_valid = True + if ell_left is not None and ell_right is not None: + left = np.asarray(ell_left, dtype=int) + right = np.asarray(ell_right, dtype=int) + if left.shape == right.shape and left.shape[0] == nbins: + for lo, hi in zip(left, right): + lo_i = int(lo) + hi_i = int(hi) + if hi_i > lo_i >= 0: + ell_lists.append(np.arange(lo_i, hi_i)) + else: + lists_valid = False + break + else: + lists_valid = False + else: + lists_valid = False + + if not lists_valid or len(ell_lists) != nbins: + ell_lists = [] + + ell_eff = None + if ell_eff_raw is not None: + ell_eff = np.asarray(ell_eff_raw, dtype=float) + if ell_eff.size != nbins: + ell_eff = None + + if ell_eff is None and ell_lists and len(ell_lists) == nbins: + ell_eff = np.asarray( + [ + 0.5 * (float(group[0]) + float(group[-1])) if group.size else 0.0 + for group in ell_lists + ], + dtype=float, + ) + + if ell_eff is not None and ell_eff.size == nbins: + return _BinInfo(ell_effective=ell_eff, ell_lists=ell_lists) + + return None + + def _load_bin_info( prereg: Path, *, @@ -186,36 +238,38 @@ def main(argv: Iterable[str] | None = None) -> int: nbins = cl_a.size nside = int(order_a.get("nside", order_b.get("nside", 0))) - fallback_nlb = args.nlb - if fallback_nlb is None: - fallback_nlb = _extract_scalar(order_a.get("nlb")) or _extract_scalar(order_b.get("nlb")) - meta_a: Mapping[str, object] | None = None - meta_b: Mapping[str, object] | None = None - if fallback_nlb is None or args.lmin is None: - meta_a = _load_order_metadata(Path(args.order_a)) - meta_b = _load_order_metadata(Path(args.order_b)) - - if fallback_nlb is None: - fallback_nlb = _extract_scalar(_extract_from_metadata(meta_a or {}, "nlb")) or _extract_scalar( - _extract_from_metadata(meta_b or {}, "nlb") + info = _bin_info_from_orders((order_a, order_b), nbins) + if info is None: + fallback_nlb = args.nlb + if fallback_nlb is None: + fallback_nlb = _extract_scalar(order_a.get("nlb")) or _extract_scalar(order_b.get("nlb")) + meta_a: Mapping[str, object] | None = None + meta_b: Mapping[str, object] | None = None + if fallback_nlb is None or args.lmin is None: + meta_a = _load_order_metadata(Path(args.order_a)) + meta_b = _load_order_metadata(Path(args.order_b)) + + if fallback_nlb is None: + fallback_nlb = _extract_scalar( + _extract_from_metadata(meta_a or {}, "nlb") + ) or _extract_scalar(_extract_from_metadata(meta_b or {}, "nlb")) + + fallback_lmin = args.lmin + if fallback_lmin is None: + fallback_lmin = _extract_scalar(order_a.get("lmin")) or _extract_scalar(order_b.get("lmin")) + if fallback_lmin is None and meta_a is not None and meta_b is not None: + fallback_lmin = _extract_scalar( + _extract_from_metadata(meta_a, "lmin") + ) or _extract_scalar(_extract_from_metadata(meta_b, "lmin")) + + info = _load_bin_info( + args.prereg, + nside=nside if nside > 0 else 256, + nbins=nbins, + fallback_lmin=fallback_lmin, + fallback_nlb=fallback_nlb, ) - fallback_lmin = args.lmin - if fallback_lmin is None: - fallback_lmin = _extract_scalar(order_a.get("lmin")) or _extract_scalar(order_b.get("lmin")) - if fallback_lmin is None and meta_a is not None and meta_b is not None: - fallback_lmin = _extract_scalar(_extract_from_metadata(meta_a, "lmin")) or _extract_scalar( - _extract_from_metadata(meta_b, "lmin") - ) - - info = _load_bin_info( - args.prereg, - nside=nside if nside > 0 else 256, - nbins=nbins, - fallback_lmin=fallback_lmin, - fallback_nlb=fallback_nlb, - ) - ell_eff = info.ell_effective cl_data = 0.5 * (cl_a + cl_b) diff --git a/scripts/run_order_A_to_B.py b/scripts/run_order_A_to_B.py index 5332a9c..470c179 100755 --- a/scripts/run_order_A_to_B.py +++ b/scripts/run_order_A_to_B.py @@ -69,6 +69,8 @@ def main(): bins = None bins_meta = None + bins_ell_groups: list[np.ndarray] = [] + ell_effective: np.ndarray | None = None windows_cfg: WindowConfig | None = None if not args.disable_prereg: try: @@ -86,25 +88,31 @@ def main(): nlb=args.nlb, ) + if hasattr(bins, "get_effective_ells"): + try: + ell_effective = np.asarray(bins.get_effective_ells(), dtype=float) + except Exception: + ell_effective = None + + if hasattr(bins, "get_ell_list"): + try: + bins_ell_groups = [np.asarray(group, dtype=int) for group in bins.get_ell_list()] + except Exception: + bins_ell_groups = [] + def _infer_bin_limits() -> tuple[int | None, int | None]: lmin_val = getattr(bins, "lmin", None) lmax_val = getattr(bins, "lmax", None) - ell_list = getattr(bins, "get_ell_list", None) - if callable(ell_list): - try: - groups = list(ell_list()) - except Exception: - groups = [] - for group in groups: - arr = np.asarray(group, dtype=int) - if arr.size: - lmin_val = int(arr[0]) - break - for group in reversed(groups): - arr = np.asarray(group, dtype=int) - if arr.size: - lmax_val = int(arr[-1]) - break + for group in bins_ell_groups: + arr = np.asarray(group, dtype=int) + if arr.size: + lmin_val = int(arr[0]) + break + for group in reversed(bins_ell_groups): + arr = np.asarray(group, dtype=int) + if arr.size: + lmax_val = int(arr[-1]) + break return ( None if lmin_val is None else int(lmin_val), None if lmax_val is None else int(lmax_val), @@ -153,11 +161,33 @@ def _infer_bin_limits() -> tuple[int | None, int | None]: # downstream stages can always reconstruct the band edges when # preregistered metadata is unavailable. lmin_used = 0 - payload: dict[str, object] = {"cl": cl, "nside": nside, "nlb": int(args.nlb)} + bin_left_edges: np.ndarray | None = None + bin_right_edges: np.ndarray | None = None + if bins_ell_groups and len(bins_ell_groups) == cl.size: + bin_left_edges = np.full(cl.size, -1, dtype=int) + bin_right_edges = np.full(cl.size, -1, dtype=int) + for idx, group in enumerate(bins_ell_groups): + if group.size: + bin_left_edges[idx] = int(group[0]) + bin_right_edges[idx] = int(group[-1]) + 1 + + nlb_value = int(args.nlb) + if isinstance(bins_meta, Mapping) and "nlb" in bins_meta: + try: + nlb_value = int(bins_meta["nlb"]) + except Exception: + nlb_value = int(args.nlb) + + payload: dict[str, object] = {"cl": cl, "nside": nside, "nlb": np.array(nlb_value, dtype=int)} if lmin_used is not None: payload["lmin"] = np.array(lmin_used, dtype=int) if lmax_used is not None: payload["lmax"] = np.array(lmax_used, dtype=int) + if ell_effective is not None and ell_effective.size == cl.size: + payload["ell_effective"] = ell_effective.astype(float) + if bin_left_edges is not None and bin_right_edges is not None: + payload["ell_left_edges"] = bin_left_edges + payload["ell_right_edges"] = bin_right_edges np.savez(Path(args.out), **payload) bins_summary: dict[str, object] = {} @@ -170,6 +200,16 @@ def _infer_bin_limits() -> tuple[int | None, int | None]: if "nlb" not in bins_summary: bins_summary["nlb"] = int(args.nlb) bins_summary["nbins"] = int(cl.size) + if ell_effective is not None and ell_effective.size == cl.size and "ell_effective" not in bins_summary: + bins_summary["ell_effective"] = ell_effective.astype(float).tolist() + if ( + bin_left_edges is not None + and bin_right_edges is not None + and "ell_left_edges" not in bins_summary + and "ell_right_edges" not in bins_summary + ): + bins_summary["ell_left_edges"] = bin_left_edges.astype(int).tolist() + bins_summary["ell_right_edges"] = bin_right_edges.astype(int).tolist() save_json( { diff --git a/scripts/run_order_B_to_A.py b/scripts/run_order_B_to_A.py index a6a3d02..fd0f70a 100755 --- a/scripts/run_order_B_to_A.py +++ b/scripts/run_order_B_to_A.py @@ -69,6 +69,8 @@ def main(): bins = None bins_meta = None + bins_ell_groups: list[np.ndarray] = [] + ell_effective: np.ndarray | None = None windows_cfg: WindowConfig | None = None if not args.disable_prereg: try: @@ -86,25 +88,31 @@ def main(): nlb=args.nlb, ) + if hasattr(bins, "get_effective_ells"): + try: + ell_effective = np.asarray(bins.get_effective_ells(), dtype=float) + except Exception: + ell_effective = None + + if hasattr(bins, "get_ell_list"): + try: + bins_ell_groups = [np.asarray(group, dtype=int) for group in bins.get_ell_list()] + except Exception: + bins_ell_groups = [] + def _infer_bin_limits() -> tuple[int | None, int | None]: lmin_val = getattr(bins, "lmin", None) lmax_val = getattr(bins, "lmax", None) - ell_list = getattr(bins, "get_ell_list", None) - if callable(ell_list): - try: - groups = list(ell_list()) - except Exception: - groups = [] - for group in groups: - arr = np.asarray(group, dtype=int) - if arr.size: - lmin_val = int(arr[0]) - break - for group in reversed(groups): - arr = np.asarray(group, dtype=int) - if arr.size: - lmax_val = int(arr[-1]) - break + for group in bins_ell_groups: + arr = np.asarray(group, dtype=int) + if arr.size: + lmin_val = int(arr[0]) + break + for group in reversed(bins_ell_groups): + arr = np.asarray(group, dtype=int) + if arr.size: + lmax_val = int(arr[-1]) + break return ( None if lmin_val is None else int(lmin_val), None if lmax_val is None else int(lmax_val), @@ -152,11 +160,33 @@ def _infer_bin_limits() -> tuple[int | None, int | None]: # cross-spectrum step can fall back to linear binning even when # preregistration metadata is absent. lmin_used = 0 - payload: dict[str, object] = {"cl": cl, "nside": nside, "nlb": int(args.nlb)} + bin_left_edges: np.ndarray | None = None + bin_right_edges: np.ndarray | None = None + if bins_ell_groups and len(bins_ell_groups) == cl.size: + bin_left_edges = np.full(cl.size, -1, dtype=int) + bin_right_edges = np.full(cl.size, -1, dtype=int) + for idx, group in enumerate(bins_ell_groups): + if group.size: + bin_left_edges[idx] = int(group[0]) + bin_right_edges[idx] = int(group[-1]) + 1 + + nlb_value = int(args.nlb) + if isinstance(bins_meta, Mapping) and "nlb" in bins_meta: + try: + nlb_value = int(bins_meta["nlb"]) + except Exception: + nlb_value = int(args.nlb) + + payload: dict[str, object] = {"cl": cl, "nside": nside, "nlb": np.array(nlb_value, dtype=int)} if lmin_used is not None: payload["lmin"] = np.array(lmin_used, dtype=int) if lmax_used is not None: payload["lmax"] = np.array(lmax_used, dtype=int) + if ell_effective is not None and ell_effective.size == cl.size: + payload["ell_effective"] = ell_effective.astype(float) + if bin_left_edges is not None and bin_right_edges is not None: + payload["ell_left_edges"] = bin_left_edges + payload["ell_right_edges"] = bin_right_edges np.savez(Path(args.out), **payload) bins_summary: dict[str, object] = {} @@ -169,6 +199,16 @@ def _infer_bin_limits() -> tuple[int | None, int | None]: if "nlb" not in bins_summary: bins_summary["nlb"] = int(args.nlb) bins_summary["nbins"] = int(cl.size) + if ell_effective is not None and ell_effective.size == cl.size and "ell_effective" not in bins_summary: + bins_summary["ell_effective"] = ell_effective.astype(float).tolist() + if ( + bin_left_edges is not None + and bin_right_edges is not None + and "ell_left_edges" not in bins_summary + and "ell_right_edges" not in bins_summary + ): + bins_summary["ell_left_edges"] = bin_left_edges.astype(int).tolist() + bins_summary["ell_right_edges"] = bin_right_edges.astype(int).tolist() save_json( { diff --git a/tests/test_compute_cross_spectrum.py b/tests/test_compute_cross_spectrum.py index fbb7be9..9d1adbe 100644 --- a/tests/test_compute_cross_spectrum.py +++ b/tests/test_compute_cross_spectrum.py @@ -26,10 +26,31 @@ def test_compute_cross_spectrum_outputs(tmp_path, monkeypatch): cl_a = np.array([1.0, 2.0, 3.0], dtype=float) cl_b = np.array([1.2, 1.8, 3.2], dtype=float) + ell_left = np.array([10, 15, 20], dtype=int) + ell_right = np.array([15, 20, 25], dtype=int) + ell_eff = np.array([12.0, 17.0, 22.0], dtype=float) order_a = tmp_path / "order_a.npz" order_b = tmp_path / "order_b.npz" - np.savez(order_a, cl=cl_a, nside=256, nlb=5, lmin=10) - np.savez(order_b, cl=cl_b, nside=256, nlb=5, lmin=10) + np.savez( + order_a, + cl=cl_a, + nside=256, + nlb=5, + lmin=10, + ell_left_edges=ell_left, + ell_right_edges=ell_right, + ell_effective=ell_eff, + ) + np.savez( + order_b, + cl=cl_b, + nside=256, + nlb=5, + lmin=10, + ell_left_edges=ell_left, + ell_right_edges=ell_right, + ell_effective=ell_eff, + ) diag = np.array([0.25, 0.16, 0.09], dtype=float) cov = tmp_path / "cov.npy" @@ -146,3 +167,45 @@ def test_compute_cross_spectrum_reads_lmin_from_json_metadata(tmp_path): data = np.load(out) np.testing.assert_allclose(data["ell"], np.array([2.0, 7.0, 12.0])) + + +def test_compute_cross_spectrum_recovers_geometry_from_edges(tmp_path): + module = _load_script_module( + "compute_cross_spectrum", Path("scripts/compute_cross_spectrum.py") + ) + + cl = np.array([1.0, 2.0, 3.0], dtype=float) + order = tmp_path / "order.npz" + np.savez( + order, + cl=cl, + nside=256, + nlb=5, + ell_left_edges=np.array([8, 18, 28], dtype=int), + ell_right_edges=np.array([18, 28, 35], dtype=int), + ) + + theory_path = tmp_path / "theory.npz" + np.savez(theory_path, ell=np.arange(0, 50, dtype=float), cl_tk=0.05 * np.arange(0, 50)) + + out = tmp_path / "cross.npz" + summary = tmp_path / "summary.json" + module.main( + [ + "--order-a", + str(order), + "--order-b", + str(order), + "--theory", + str(theory_path), + "--prereg", + str(tmp_path / "missing_prereg.yaml"), + "--out", + str(out), + "--summary", + str(summary), + ] + ) + + data = np.load(out) + np.testing.assert_allclose(data["ell"], np.array([12.5, 22.5, 31.0])) From e41d52519b5ba9b30eaf9e84a2ebf519bc09f1a3 Mon Sep 17 00:00:00 2001 From: Alastair Hewitt Date: Tue, 21 Oct 2025 20:57:18 -0400 Subject: [PATCH 6/9] Ensure theory fixtures include required spectra --- tests/test_compute_cross_spectrum.py | 27 ++++++++++++++++++++++++--- 1 file changed, 24 insertions(+), 3 deletions(-) diff --git a/tests/test_compute_cross_spectrum.py b/tests/test_compute_cross_spectrum.py index 9d1adbe..bb278e4 100644 --- a/tests/test_compute_cross_spectrum.py +++ b/tests/test_compute_cross_spectrum.py @@ -114,7 +114,14 @@ def test_compute_cross_spectrum_requires_lmin_when_falling_back(tmp_path): np.savez(order, cl=cl, nside=256, nlb=10) theory_path = tmp_path / "theory.npz" - np.savez(theory_path, ell=np.arange(0, 40, dtype=float), cl_tk=0.05 * np.arange(0, 40)) + ell = np.arange(0, 40, dtype=float) + np.savez( + theory_path, + ell=ell, + cl_tt=0.2 * ell, + cl_kk=0.1 * ell, + cl_tk=0.05 * ell, + ) with pytest.raises(ValueError, match="provide --lmin"): module.main( @@ -144,7 +151,14 @@ def test_compute_cross_spectrum_reads_lmin_from_json_metadata(tmp_path): (order.with_suffix(".json")).write_text(json.dumps({"bins": {"lmin": 0, "nlb": 5}})) theory_path = tmp_path / "theory.npz" - np.savez(theory_path, ell=np.arange(0, 40, dtype=float), cl_tk=0.05 * np.arange(0, 40)) + ell = np.arange(0, 40, dtype=float) + np.savez( + theory_path, + ell=ell, + cl_tt=0.2 * ell, + cl_kk=0.1 * ell, + cl_tk=0.05 * ell, + ) out = tmp_path / "cross.npz" summary = tmp_path / "summary.json" @@ -186,7 +200,14 @@ def test_compute_cross_spectrum_recovers_geometry_from_edges(tmp_path): ) theory_path = tmp_path / "theory.npz" - np.savez(theory_path, ell=np.arange(0, 50, dtype=float), cl_tk=0.05 * np.arange(0, 50)) + ell = np.arange(0, 50, dtype=float) + np.savez( + theory_path, + ell=ell, + cl_tt=0.2 * ell, + cl_kk=0.1 * ell, + cl_tk=0.05 * ell, + ) out = tmp_path / "cross.npz" summary = tmp_path / "summary.json" From b5da3fb6538c277fbbc54900b4da1ce67c5c0b6f Mon Sep 17 00:00:00 2001 From: Alastair Hewitt Date: Tue, 21 Oct 2025 21:08:00 -0400 Subject: [PATCH 7/9] Fix cross-spectrum tests to use local covariance --- tests/test_compute_cross_spectrum.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/tests/test_compute_cross_spectrum.py b/tests/test_compute_cross_spectrum.py index bb278e4..6fa7013 100644 --- a/tests/test_compute_cross_spectrum.py +++ b/tests/test_compute_cross_spectrum.py @@ -160,6 +160,9 @@ def test_compute_cross_spectrum_reads_lmin_from_json_metadata(tmp_path): cl_tk=0.05 * ell, ) + cov = tmp_path / "cov.npy" + np.save(cov, np.eye(3, dtype=float)) + out = tmp_path / "cross.npz" summary = tmp_path / "summary.json" module.main( @@ -172,6 +175,8 @@ def test_compute_cross_spectrum_reads_lmin_from_json_metadata(tmp_path): str(theory_path), "--prereg", str(tmp_path / "missing_prereg.yaml"), + "--cov", + str(cov), "--out", str(out), "--summary", @@ -209,6 +214,9 @@ def test_compute_cross_spectrum_recovers_geometry_from_edges(tmp_path): cl_tk=0.05 * ell, ) + cov = tmp_path / "cov.npy" + np.save(cov, np.eye(3, dtype=float)) + out = tmp_path / "cross.npz" summary = tmp_path / "summary.json" module.main( @@ -221,6 +229,8 @@ def test_compute_cross_spectrum_recovers_geometry_from_edges(tmp_path): str(theory_path), "--prereg", str(tmp_path / "missing_prereg.yaml"), + "--cov", + str(cov), "--out", str(out), "--summary", From 81fe3f4bfabbff1e7fe21a2b5e5991a5279f5c16 Mon Sep 17 00:00:00 2001 From: Alastair Hewitt Date: Tue, 21 Oct 2025 21:09:03 -0400 Subject: [PATCH 8/9] fix lint --- scripts/compute_cross_spectrum.py | 8 ++++++-- scripts/run_order_A_to_B.py | 6 +++++- scripts/run_order_B_to_A.py | 6 +++++- 3 files changed, 16 insertions(+), 4 deletions(-) diff --git a/scripts/compute_cross_spectrum.py b/scripts/compute_cross_spectrum.py index dae934c..303a36f 100644 --- a/scripts/compute_cross_spectrum.py +++ b/scripts/compute_cross_spectrum.py @@ -242,7 +242,9 @@ def main(argv: Iterable[str] | None = None) -> int: if info is None: fallback_nlb = args.nlb if fallback_nlb is None: - fallback_nlb = _extract_scalar(order_a.get("nlb")) or _extract_scalar(order_b.get("nlb")) + fallback_nlb = _extract_scalar(order_a.get("nlb")) or _extract_scalar( + order_b.get("nlb") + ) meta_a: Mapping[str, object] | None = None meta_b: Mapping[str, object] | None = None if fallback_nlb is None or args.lmin is None: @@ -256,7 +258,9 @@ def main(argv: Iterable[str] | None = None) -> int: fallback_lmin = args.lmin if fallback_lmin is None: - fallback_lmin = _extract_scalar(order_a.get("lmin")) or _extract_scalar(order_b.get("lmin")) + fallback_lmin = _extract_scalar(order_a.get("lmin")) or _extract_scalar( + order_b.get("lmin") + ) if fallback_lmin is None and meta_a is not None and meta_b is not None: fallback_lmin = _extract_scalar( _extract_from_metadata(meta_a, "lmin") diff --git a/scripts/run_order_A_to_B.py b/scripts/run_order_A_to_B.py index 470c179..6d04b1d 100755 --- a/scripts/run_order_A_to_B.py +++ b/scripts/run_order_A_to_B.py @@ -200,7 +200,11 @@ def _infer_bin_limits() -> tuple[int | None, int | None]: if "nlb" not in bins_summary: bins_summary["nlb"] = int(args.nlb) bins_summary["nbins"] = int(cl.size) - if ell_effective is not None and ell_effective.size == cl.size and "ell_effective" not in bins_summary: + if ( + ell_effective is not None + and ell_effective.size == cl.size + and "ell_effective" not in bins_summary + ): bins_summary["ell_effective"] = ell_effective.astype(float).tolist() if ( bin_left_edges is not None diff --git a/scripts/run_order_B_to_A.py b/scripts/run_order_B_to_A.py index fd0f70a..3c437aa 100755 --- a/scripts/run_order_B_to_A.py +++ b/scripts/run_order_B_to_A.py @@ -199,7 +199,11 @@ def _infer_bin_limits() -> tuple[int | None, int | None]: if "nlb" not in bins_summary: bins_summary["nlb"] = int(args.nlb) bins_summary["nbins"] = int(cl.size) - if ell_effective is not None and ell_effective.size == cl.size and "ell_effective" not in bins_summary: + if ( + ell_effective is not None + and ell_effective.size == cl.size + and "ell_effective" not in bins_summary + ): bins_summary["ell_effective"] = ell_effective.astype(float).tolist() if ( bin_left_edges is not None From 3b68fecec418d8929101b8f35293d854c4f85ea7 Mon Sep 17 00:00:00 2001 From: Alastair Hewitt Date: Tue, 21 Oct 2025 21:10:34 -0400 Subject: [PATCH 9/9] fix indent --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 2dd2ea3..f3b33f1 100644 --- a/README.md +++ b/README.md @@ -202,7 +202,7 @@ earlier. --cov artifacts/cov_delta_full.npy \ --out-delta artifacts/delta_ell_full.npy \ --out-summary artifacts/summary_full.json - ``` + ``` The resulting JSON contains the Δ vector length and the stabilized χ ("z") statistic for the null test. When reusing an older covariance that has one or two extra high-ℓ bins, add