API(maps): add deconvolve() method to mapper classes

heracles/maps/_healpix.py

@@ -233,3 +233,24 @@ def transform(
             update_metadata(alms[1], **md)
 
         return alms
+
+    def deconvolve(self, alm: ArrayLike, *, inplace: bool = False) -> ArrayLike:
+        """
+        Remove HEALPix pixel window function from *alm*.
+        """
+
+        lmax = hp.Alm.getlmax(alm.size)
+
+        md = alm.dtype.metadata or {}
+        spin = md.get("spin", 0)
+
+        if spin == 0:
+            pw = hp.pixwin(self.__nside, lmax=lmax)
+        elif spin == 2:
+            _, pw = hp.pixwin(self.__nside, lmax=lmax, pol=True)
+            pw[:2] = 1.0
+        else:
+            msg = f"unsupported spin for deconvolve: {spin}"
+            raise ValueError(msg)
+
+        return hp.almxfl(alm, 1 / pw, inplace=inplace)

heracles/maps/_mapper.py

@@ -126,3 +126,9 @@ def transform(
         """
         The spherical harmonic transform for this mapper.
         """
+
+    @abstractmethod
+    def deconvolve(self, alm: ArrayLike, *, inplace: bool = False) -> ArrayLike:
+        """
+        Remove this mapper's convolution kernel from *alm*.
+        """

heracles/maps/_mapping.py

@@ -147,6 +147,7 @@ def transform_maps(
     maps: Mapping[tuple[Any, Any], NDArray],
     *,
     lmax: int | Mapping[Any, int] | None = None,
+    deconvolve: bool = True,
     out: MutableMapping[tuple[Any, Any], NDArray] | None = None,
     progress: bool = False,
     **kwargs,
@@ -187,6 +188,13 @@ def transform_maps(
 
             alms = _mapper.transform(m, _lmax)
 
+            if deconvolve:
+                if isinstance(alms, tuple):
+                    for alm in alms:
+                        _mapper.deconvolve(alm, inplace=True)
+                else:
+                    _mapper.deconvolve(alms, inplace=True)
+
             if isinstance(alms, tuple):
                 out[f"{k}_E", i] = alms[0]
                 out[f"{k}_B", i] = alms[1]

heracles/twopoint.py

@@ -182,80 +182,6 @@ def debias_cls(cls, bias=None, *, inplace=False):
     return out
 
 
-def depixelate_cls(cls, *, inplace=False):
-    """remove discretisation kernel from cls"""
-
-    logger.info("depixelate %d cl(s)%s", len(cls), " in place" if inplace else "")
-    t = time.monotonic()
-
-    # keep a cache of convolution kernels (i.e. pixel window functions)
-    fls = {
-        "healpix": {},
-    }
-
-    # the output toc dict
-    out = cls if inplace else TocDict()
-
-    # remove effect of convolution for each cl in turn
-    for key in cls:
-        logger.info("depixelate %s x %s cl for bins %s, %s", *key)
-
-        cl = cls[key]
-        md = cl.dtype.metadata or {}
-
-        if not inplace:
-            cl = cl.copy()
-            update_metadata(cl, **md)
-
-        lmax = len(cl) - 1
-
-        spins = [md.get("spin_1", 0), md.get("spin_2", 0)]
-        kernels = [md.get("kernel_1"), md.get("kernel_2")]
-
-        # minimum l for corrections
-        lmin = max(map(abs, spins))
-
-        # deconvolve the kernels of the first and second map
-        for i, spin, kernel in zip([1, 2], spins, kernels):
-            logger.info("- spin-%s %s kernel", spin, kernel)
-            if kernel is None:
-                fl = None
-            elif kernel == "healpix":
-                nside = md[f"nside_{i}"]
-                if (nside, lmax, spin) not in fls[kernel]:
-                    fl0, fl2 = hp.pixwin(nside, lmax=lmax, pol=True)
-                    fls[kernel][nside, lmax, 0] = fl0
-                    fls[kernel][nside, lmax, 2] = fl2
-                fl = fls[kernel].get((nside, lmax, spin))
-                if fl is None:
-                    logger.warning(
-                        "no HEALPix kernel for NSIDE = %s, LMAX = %s, SPIN = %s",
-                        nside,
-                        lmax,
-                        spin,
-                    )
-            else:
-                msg = f"unknown kernel: {kernel}"
-                raise ValueError(msg)
-            if fl is not None:
-                if cl.dtype.names is None:
-                    cl[lmin:] /= fl[lmin:]
-                else:
-                    cl["CL"][lmin:] /= fl[lmin:]
-
-        # store depixelated cl in output set
-        out[key] = cl
-
-    logger.info(
-        "depixelated %d cl(s) in %s",
-        len(out),
-        timedelta(seconds=(time.monotonic() - t)),
-    )
-
-    # return the toc dict of depixelated cls
-    return out
-
-
 def mixing_matrices(cls, *, l1max=None, l2max=None, l3max=None):
     """compute mixing matrices from a set of cls"""
 
@@ -317,64 +243,6 @@ def mixing_matrices(cls, *, l1max=None, l2max=None, l3max=None):
     return mms
 
 
-def pixelate_mms_healpix(mms, nside, *, inplace=False):
-    """apply HEALPix pixel window function to mms"""
-
-    logger.info("pixelate %d mm(s)%s", len(mms), " in place" if inplace else "")
-    logger.info("kernel: HEALPix, NSIDE=%d", nside)
-    t = time.monotonic()
-
-    # pixel window functions
-    lmax = 4 * nside
-    fl0, fl2 = hp.pixwin(nside, lmax=lmax, pol=True)
-
-    # will be multiplied over rows
-    fl0 = fl0[:, np.newaxis]
-    fl2 = fl2[:, np.newaxis]
-
-    # the output toc dict
-    out = mms if inplace else TocDict()
-
-    # apply discretisation kernel from cl to each mm in turn
-    for key in mms:
-        logger.info("pixelate %s mm for bins %s, %s", *key)
-
-        mm = mms[key]
-        if not inplace:
-            mm = mm.copy()
-
-        n = np.shape(mm)[-2]
-        if n >= lmax:
-            logger.error(
-                "no HEALPix pixel window function for NSIDE=%d and LMAX=%d",
-                nside,
-                n - 1,
-            )
-
-        name = key[0]
-        if name == "00":
-            mm *= fl0[:n] * fl0[:n]
-        elif name in ["0+", "+0"]:
-            mm *= fl0[:n] * fl2[:n]
-        elif name in ["++", "--", "+-", "-+"]:
-            mm *= fl2[:n] * fl2[:n]
-        else:
-            logger.warning("unknown mixing matrix, assuming spin-0")
-            mm *= fl0[:n] * fl0[:n]
-
-        # store pixelated mm in output set
-        out[key] = mm
-
-    logger.info(
-        "pixelated %d mm(s) in %s",
-        len(out),
-        timedelta(seconds=(time.monotonic() - t)),
-    )
-
-    # return the toc dict of modified cls
-    return out
-
-
 def bin2pt(arr, bins, name, *, weights=None):
     """Compute binned two-point data."""
 

pyproject.toml

@@ -93,7 +93,6 @@ select = [
     "ICN",
     "INP",
     "ISC",
-    "PD",
     "PIE",
     "PYI",
     "Q",

tests/test_twopoint.py

@@ -164,35 +164,6 @@ def test_mixing_matrices(rng):
     assert mms["XY", 0, 1].shape == (lmax + 1, lmax + 1)
 
 
-def test_pixelate_mms_healpix():
-    import healpy as hp
-
-    from heracles.twopoint import pixelate_mms_healpix
-
-    nside = 512
-    lmax = 1000
-
-    fl0, fl2 = hp.pixwin(nside, lmax=lmax, pol=True)
-
-    mms = {
-        ("00", 0, 0): np.eye(lmax + 1),
-        ("0+", 0, 0): np.eye(lmax + 1),
-        ("++", 0, 0): np.eye(lmax + 1),
-        ("--", 0, 0): np.eye(lmax + 1),
-        ("+-", 0, 0): np.eye(lmax + 1),
-        ("ab", 0, 0): np.eye(lmax + 1),
-    }
-
-    pixelate_mms_healpix(mms, nside, inplace=True)
-
-    assert np.all(mms["00", 0, 0] == np.diag(fl0 * fl0))
-    assert np.all(mms["0+", 0, 0] == np.diag(fl0 * fl2))
-    assert np.all(mms["++", 0, 0] == np.diag(fl2 * fl2))
-    assert np.all(mms["--", 0, 0] == np.diag(fl2 * fl2))
-    assert np.all(mms["+-", 0, 0] == np.diag(fl2 * fl2))
-    assert np.all(mms["ab", 0, 0] == np.diag(fl0 * fl0))
-
-
 @pytest.mark.parametrize("weights", [None, "l(l+1)", "2l+1", "<rand>"])
 def test_binned_cls(rng, weights):
     from heracles.twopoint import binned_cls