From 0152aa22824636b8a715bfca4d4a60a6ef223a2e Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Thu, 1 Oct 2026 04:37:56 +0200 Subject: [PATCH] Fix one-multipole shift in pseudo-Cl covariance theory NaMaster arrays for fields built with lmax=b_lmax run over l = 0..b_lmax with index equal to l. The iNKA covariance built its fiducial C_l (and pixel window) on l = 1..lmax, so the theory at l landed at index l-1: every input multipole was shifted down by one (41% at l=2, 14% at l=3, ~2-3% at l=5-20, 1% at l=100 for a Planck18 single-bin spectrum). Add pseudo_cl.pixelised_theory_cl, which evaluates a theory C_l on NaMaster's l = 0..b_lmax grid (zero at l < 2) times pw^2(l), and use it in calculate_pseudo_cl_eb_cov with the shared pseudo_cl_geometry. It also takes the single-bin "W1xW1" entry from cs_util's get_theo_c_ell, which returns a dict keyed by bin pair; the covariance method raised TypeError without it. Label the coupled-spectrum ell row in the GLASS mock pseudo-Cl helpers as 0..lmax-1, matching compute_coupled_cell. Tests pin the l <-> index convention against NaMaster's own coupling and smoke-test the covariance method end to end on the synthetic catalog. Co-Authored-By: Claude Opus 5.5 --- src/sp_validation/cosmo_val/pseudo_cl.py | 25 +++++--------- src/sp_validation/glass_mock.py | 4 +-- src/sp_validation/pseudo_cl.py | 36 ++++++++++++++++++-- src/sp_validation/tests/test_pseudo_cl.py | 40 +++++++++++++++++++++++ 4 files changed, 83 insertions(+), 22 deletions(-) diff --git a/src/sp_validation/cosmo_val/pseudo_cl.py b/src/sp_validation/cosmo_val/pseudo_cl.py index cf74efc6..321acfbd 100644 --- a/src/sp_validation/cosmo_val/pseudo_cl.py +++ b/src/sp_validation/cosmo_val/pseudo_cl.py @@ -24,6 +24,8 @@ get_pseudo_cls_catalog, get_pseudo_cls_map, make_namaster_bin, + pixelised_theory_cl, + pseudo_cl_geometry, ) from ..rho_tau import get_params_rho_tau from ..statistics import chi2_and_pte, cov_from_one_covariance @@ -157,36 +159,25 @@ def calculate_pseudo_cl_eb_cov(self): self.print_cyan(f"Extracting the fiducial power spectrum for {ver}") - lmax = 2 * self.nside - ell = np.arange(1, lmax + 1) - pw = hp.pixwin(nside, lmax=lmax) - if pw.shape[0] != len(ell) + 1: - raise ValueError( - "Unexpected pixwin length for lmax=" - f"{lmax}: got {pw.shape[0]}, expected {len(ell) + 1}" - ) - pw = pw[1 : len(ell) + 1] + lmin, lmax, b_lmax = pseudo_cl_geometry(nside) # Load redshift distribution and calculate theory C_ell path_redshift_distr = self.cc[ver]["shear"]["redshift_path"] z, dndz = np.loadtxt(path_redshift_distr, unpack=True) - fiducial_cl = ( - get_theo_c_ell( + fiducial_cl = pixelised_theory_cl( + lambda ell: get_theo_c_ell( ell=ell, z=z, nz=dndz, backend="ccl", cosmo=self.cosmo, - ) - * pw**2 + )["W1xW1"], + nside, + b_lmax, ) self.print_cyan("Getting a binning, n_gal_map, field and workspace.") - lmin = 8 - lmax = 2 * self.nside - b_lmax = lmax - 1 - b = self.get_namaster_bin(lmin, lmax, b_lmax) # Load data and create shear and noise maps diff --git a/src/sp_validation/glass_mock.py b/src/sp_validation/glass_mock.py index de8312c7..43e8df9a 100644 --- a/src/sp_validation/glass_mock.py +++ b/src/sp_validation/glass_mock.py @@ -423,7 +423,7 @@ def compute_two_point_cl(cat, nside=1024, lmin=8, n_bins=32): cl_coupled = nmt.compute_coupled_cell(f_all, f_all) cl_all = wsp.decouple_cell(cl_coupled) - cl_coupled = np.concatenate([np.arange(1, lmax + 1)[np.newaxis, :], cl_coupled]) + cl_coupled = np.concatenate([np.arange(lmax)[np.newaxis, :], cl_coupled]) cl_all = np.concatenate([ell_eff[np.newaxis, ...], cl_all]) return cl_coupled, cl_all @@ -478,7 +478,7 @@ def compute_two_point_cl_map(cat, nside=1024, lmin=8, n_bins=32): cl_coupled = nmt.compute_coupled_cell(f_all, f_all) cl_all = wsp.decouple_cell(cl_coupled) - cl_coupled = np.concatenate([np.arange(1, lmax + 1)[np.newaxis, :], cl_coupled]) + cl_coupled = np.concatenate([np.arange(lmax)[np.newaxis, :], cl_coupled]) cl_all = np.concatenate([ell_eff[np.newaxis, ...], cl_all]) return cl_coupled, cl_all diff --git a/src/sp_validation/pseudo_cl.py b/src/sp_validation/pseudo_cl.py index 34cbc68d..a76581a0 100644 --- a/src/sp_validation/pseudo_cl.py +++ b/src/sp_validation/pseudo_cl.py @@ -4,9 +4,9 @@ mirroring ``b_modes.py`` and ``rho_tau.py``: the orchestrator mixin in ``sp_validation.cosmo_val.pseudo_cl`` (and analysis scripts directly) call these free functions. Everything here is pure computation -- NaMaster binning, -weighted galaxy number-density maps, random-rotation noise debiasing, and the -map-/catalog-based pseudo-Cl estimators. Depends on pymaster (NaMaster) and -healpy. +theory C_ℓ on NaMaster's multipole grid, weighted galaxy number-density maps, +random-rotation noise debiasing, and the map-/catalog-based pseudo-Cl +estimators. Depends on pymaster (NaMaster) and healpy. The harmonic geometry the estimators use is fixed by ``nside``: ``lmin = 8``, ``lmax = 2 * nside``, ``b_lmax = lmax - 1``. ``pseudo_cl_geometry`` returns that @@ -31,6 +31,36 @@ def pseudo_cl_geometry(nside): return LMIN, lmax, lmax - 1 +def pixelised_theory_cl(cl_of_ell, nside, b_lmax): + """Theory C_ℓ on NaMaster's multipole grid, times the pixel window pw²(ℓ). + + Arrays NaMaster takes or returns for fields built with ``lmax=b_lmax`` + (coupled spectra, ``couple_cell``, ``gaussian_covariance`` inputs, the + coupling matrix) run over ℓ = 0..b_lmax, with index equal to ℓ. This + returns a C_ℓ on that grid: ``cl_of_ell`` evaluated at ℓ ≥ 2 and multiplied + by the HEALPix pixel window squared at ``nside``; entries ℓ = 0, 1 are zero + (spin-2 fields carry no monopole or dipole). + + Parameters + ---------- + cl_of_ell : callable + Maps an integer multipole array to the C_ℓ at those multipoles. + nside : int + HEALPix resolution of the maps the spectrum is measured on. + b_lmax : int + Maximum multipole of the NaMaster fields. + + Returns + ------- + np.ndarray + Shape ``(b_lmax + 1,)``; element ℓ is pw²(ℓ) C_ℓ. + """ + ell = np.arange(b_lmax + 1) + cl = np.zeros(ell.size) + cl[2:] = cl_of_ell(ell[2:]) + return cl * hp.pixwin(nside, lmax=b_lmax) ** 2 + + def make_namaster_bin( lmin, lmax, b_lmax, binning, *, ell_step=10, n_ell_bins=32, power=0.5 ): diff --git a/src/sp_validation/tests/test_pseudo_cl.py b/src/sp_validation/tests/test_pseudo_cl.py index ca6bbb1f..31cad7ea 100644 --- a/src/sp_validation/tests/test_pseudo_cl.py +++ b/src/sp_validation/tests/test_pseudo_cl.py @@ -612,3 +612,43 @@ def test_calculate_pseudo_cl_out_path_rejects_multiversion(cv): cv.versions = [cv._test_version, "SecondVersion"] with pytest.raises(ValueError, match="one part to one path"): cv.calculate_pseudo_cl(out_path=cv._output_path("pseudo_cl_x.sacc")) + + +def test_pixelised_theory_cl_indexes_by_ell(): + """Element ℓ of the theory vector is pw²(ℓ) C_ℓ on NaMaster's ℓ = 0..b_lmax grid. + + The vector feeds ``couple_cell`` / ``gaussian_covariance`` directly, so it + must match NaMaster's own multipole axis: same length as a coupled spectrum, + and on the full sky (coupling = identity for ℓ ≥ 2) coupling returns it + unchanged at every ℓ. + """ + from sp_validation.pseudo_cl import pixelised_theory_cl, pseudo_cl_geometry + + nside = 16 + _lmin, _lmax, b_lmax = pseudo_cl_geometry(nside) + cl = pixelised_theory_cl(lambda ell: 1.0 / ell**2, nside, b_lmax) + + ell = np.arange(b_lmax + 1) + pw2 = healpy.pixwin(nside, lmax=b_lmax) ** 2 + npt.assert_array_equal(cl[:2], 0.0) + npt.assert_allclose(cl[2:], pw2[2:] / ell[2:] ** 2, rtol=1e-14) + + mask = np.ones(healpy.nside2npix(nside)) + f = pymaster.NmtField(mask=mask, maps=[np.zeros_like(mask)] * 2, lmax=b_lmax) + b = pymaster.NmtBin.from_lmax_linear(b_lmax, 4) + wsp = pymaster.NmtWorkspace.from_fields(f, f, b) + assert pymaster.compute_coupled_cell(f, f).shape == (4, cl.size) + + zero = np.zeros_like(cl) + coupled = wsp.couple_cell(np.array([cl, zero, zero, zero])) + npt.assert_allclose(coupled[0, 2:], cl[2:], rtol=1e-3) + + +def test_calculate_pseudo_cl_eb_cov_runs(cv): + """The iNKA covariance runs end to end and returns a symmetric, positive EE block.""" + cv.calculate_pseudo_cl_eb_cov() + cov = cv._pseudo_cls[cv._test_version]["cov"] + ee = cov["COVAR_EE_EE"].data + assert ee.shape == (N_ELL_BINS, N_ELL_BINS) + npt.assert_allclose(ee, ee.T, rtol=1e-6) + assert np.all(np.diag(ee) > 0)