Fix one-multipole shift in the pseudo-Cℓ Gaussian covariance theory

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

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
 

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
 ):

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)