API: make bias work for HEALPix and spherical convolution (#111)

Changes how the bias is handled in Heracles. All alms are now
deconvolved by default. For HEALPix, an ell-dependent bias is removed
from angular power spectra. Uses the new reduced mixing matrix
methodology to correctly handle the HEALPix pixel window function in the
mixing matrices (i.e. not at all, same as for spherical convolution).

Closes: #107

heracles/core.py

@@ -143,3 +143,11 @@ def update_metadata(array, *sources, **metadata):
         raise ValueError(msg)
     # set the new dtype in array
     array.dtype = dt
+
+
+def items_with_suffix(d: Mapping[str, Any], suffix: str) -> Mapping[str, Any]:
+    """
+    Return items from *d* where keys end in *suffix*.  Returns a mapping
+    where *suffix* is removed from keys.
+    """
+    return {k.removesuffix(suffix): v for k, v in d.items() if k.endswith(suffix)}

heracles/fields.py

@@ -320,7 +320,7 @@ async def __call__(
 
         # total weighted variance from online algorithm
         ngal = 0
-        wmean, w2mean, var = 0.0, 0.0, 0.0
+        wmean, var = 0.0, 0.0
 
         # go through pages in catalogue and map values
         async for page in _pages(catalog, progress):
@@ -338,7 +338,6 @@ async def __call__(
                     var += (v**2 - var).sum() / ngal
                 else:
                     wmean += (w - wmean).sum() / ngal
-                    w2mean += (w**2 - w2mean).sum() / ngal
                     var += ((w * v) ** 2 - var).sum() / ngal
 
                 del lon, lat, v, w
@@ -348,7 +347,7 @@ async def __call__(
 
         # fix mean weight if there was no column for it
         if wcol is None:
-            wmean = w2mean = 1.0
+            wmean = 1.0
 
         # compute mean visibility
         if catalog.visibility is None:
@@ -365,11 +364,8 @@ async def __call__(
         # compute bias from variance (per object)
         bias = 4 * np.pi * vbar**2 * (var / wmean**2) / ngal
 
-        # bias correction factor for intrinsic variance
-        bcor = 4 * np.pi * vbar**2 * (w2mean / wmean**2) / ngal
-
         # set metadata of array
-        update_metadata(val, self, catalog, mapper, wbar=wbar, bias=bias, bcor=bcor)
+        update_metadata(val, self, catalog, mapper, wbar=wbar, bias=bias)
 
         # return the value map
         return val
@@ -402,7 +398,7 @@ async def __call__(
 
         # total weighted variance from online algorithm
         ngal = 0
-        wmean, w2mean, var = 0.0, 0.0, 0.0
+        wmean, var = 0.0, 0.0
 
         # go through pages in catalogue and get the shear values,
         async for page in _pages(catalog, progress):
@@ -421,7 +417,6 @@ async def __call__(
                     var += (re**2 + im**2 - var).sum() / ngal
                 else:
                     wmean += (w - wmean).sum() / ngal
-                    w2mean += (w**2 - w2mean).sum() / ngal
                     var += ((w * re) ** 2 + (w * im) ** 2 - var).sum() / ngal
 
                 del lon, lat, re, im, w
@@ -430,7 +425,7 @@ async def __call__(
 
         # set mean weight if there was no column for it
         if wcol is None:
-            wmean = w2mean = 1.0
+            wmean = 1.0
 
         # compute mean visibility
         if catalog.visibility is None:
@@ -447,11 +442,8 @@ async def __call__(
         # bias from measured variance, for E/B decomposition
         bias = 2 * np.pi * vbar**2 * (var / wmean**2) / ngal
 
-        # bias correction factor for intrinsic field variance
-        bcor = 2 * np.pi * vbar**2 * (w2mean / wmean**2) / ngal
-
         # set metadata of array
-        update_metadata(val, self, catalog, mapper, wbar=wbar, bias=bias, bcor=bcor)
+        update_metadata(val, self, catalog, mapper, wbar=wbar, bias=bias)
 
         # return the shear map
         return val
@@ -513,29 +505,52 @@ async def __call__(
         # weight map
         wht = np.zeros(mapper.size, mapper.dtype)
 
+        # total weighted variance from online algorithm
+        ngal = 0
+        wmean, w2mean = 0.0, 0.0
+
         # map catalogue
         async for page in _pages(catalog, progress):
-            lon, lat = page.get(*col)
+            if wcol is not None:
+                page.delete(page[wcol] == 0)
 
-            if wcol is None:
-                w = None
-            else:
-                w = page.get(wcol)
+            if page.size:
+                lon, lat = page.get(*col)
+
+                w = page.get(wcol) if wcol is not None else None
+
+                mapper.map_values(lon, lat, [wht], None, w)
+
+                ngal += page.size
+                if w is not None:
+                    wmean += (w - wmean).sum() / ngal
+                    w2mean += (w**2 - w2mean).sum() / ngal
 
-            mapper.map_values(lon, lat, [wht], None, w)
+                del lon, lat, w
+
+            del page
+
+        # set mean weight if there was no column for it
+        if wcol is None:
+            wmean = w2mean = 1.0
 
-            del page, lon, lat, w
+        # compute mean visibility
+        if catalog.visibility is None:
+            vbar = 1
+        else:
+            vbar = np.mean(catalog.visibility)
 
-        # compute average weight in nonzero pixels
-        wbar = wht.mean()
-        if catalog.visibility is not None:
-            wbar /= np.mean(catalog.visibility)
+        # mean weight per effective mapper "pixel"
+        wbar = ngal / (4 * np.pi * vbar) * wmean * mapper.area
 
         # normalise the map
         wht /= wbar
 
-        # set metadata of arrays
-        update_metadata(wht, self, catalog, mapper, wbar=wbar)
+        # bias from weights
+        bias = 4 * np.pi * vbar**2 * (w2mean / wmean**2) / ngal
+
+        # set metadata of array
+        update_metadata(wht, self, catalog, mapper, wbar=wbar, bias=bias)
 
         # return the weight map
         return wht

heracles/maps/__init__.py

@@ -25,9 +25,10 @@
     "Mapper",
     "get_kernels",
     "map_catalogs",
+    "mapper_from_dict",
     "transform_maps",
 )
 
 from ._healpix import Healpix
-from ._mapper import Mapper, get_kernels
+from ._mapper import Mapper, get_kernels, mapper_from_dict
 from ._mapping import map_catalogs, transform_maps

heracles/maps/_healpix.py

@@ -34,9 +34,10 @@
 from ._mapper import Mapper
 
 if TYPE_CHECKING:
-    from collections.abc import Sequence
+    from collections.abc import Mapping, Sequence
+    from typing import Any, Self
 
-    from numpy.typing import ArrayLike, DTypeLike
+    from numpy.typing import ArrayLike, DTypeLike, NDArray
 
 
 def _asnative(arr):
@@ -120,12 +121,25 @@ class Healpix(Mapper, kernel="healpix"):
     available as the *nside* property.
     """
 
+    DATAPATH: str | None = None
+
+    @classmethod
+    def from_dict(cls, d: Mapping[str, Any]) -> Self:
+        """
+        Create a HEALPix mapper from a dictionary.
+        """
+        nside = d["nside"]
+        if isinstance(nside, str) and nside.isdigit():
+            nside = int(nside)
+        elif not isinstance(nside, int):
+            msg = f"value for key 'nside' is {type(nside).__name__}, expected int"
+            raise ValueError(msg)
+        return cls(nside)
+
     def __init__(
         self,
         nside: int,
         dtype: DTypeLike = np.float64,
-        *,
-        datapath: str | None = None,
     ) -> None:
         """
         Mapper for HEALPix maps with the given *nside* parameter.
@@ -134,7 +148,6 @@ def __init__(
         self.__nside = nside
         self.__npix = hp.nside2npix(nside)
         self.__dtype = np.dtype(dtype)
-        self.__datapath = datapath
         self._metadata |= {
             "nside": nside,
         }
@@ -222,7 +235,7 @@ def transform(
             lmax=lmax,
             pol=pol,
             use_pixel_weights=True,
-            datapath=self.__datapath,
+            datapath=self.DATAPATH,
         )
 
         if spin == 0:
@@ -234,23 +247,26 @@ def transform(
 
         return alms
 
-    def deconvolve(self, alm: ArrayLike, *, inplace: bool = False) -> ArrayLike:
+    def kl(self, lmax: int, spin: int = 0) -> NDArray[Any]:
         """
-        Remove HEALPix pixel window function from *alm*.
+        Return the HEALPix pixel window function.
         """
 
-        lmax = hp.Alm.getlmax(alm.size)
-
-        md = alm.dtype.metadata or {}
-        spin = md.get("spin", 0)
-
+        pw: NDArray[Any]
         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}"
+            msg = f"unsupported spin: {spin}"
             raise ValueError(msg)
 
-        return hp.almxfl(alm, 1 / pw, inplace=inplace)
+        return pw
+
+    def bl(self, lmax: int, spin: int = 0) -> NDArray[Any]:
+        """
+        Return the biasing kernel for HEALPix.
+        """
+        kl = self.kl(lmax, spin)
+        where = np.arange(lmax + 1) >= abs(spin)
+        return np.divide(1.0, kl, where=where, out=np.zeros(lmax + 1))

heracles/maps/_mapper.py

@@ -28,8 +28,9 @@
 
 if TYPE_CHECKING:
     from collections.abc import Mapping, Sequence
+    from typing import Self
 
-    from numpy.typing import ArrayLike, DTypeLike
+    from numpy.typing import ArrayLike, DTypeLike, NDArray
 
 # dictionary of kernel names and their corresponding Mapper classes
 _KERNELS: dict[str, type[Mapper]] = {}
@@ -43,39 +44,63 @@ def get_kernels() -> Mapping[str, type[Mapper]]:
     return MappingProxyType(_KERNELS)
 
 
-class MapperMeta(ABCMeta):
+def mapper_from_dict(d: Mapping[str, Any]) -> Mapper:
     """
-    Metaclass for mappers.
+    Return a mapper that matches the given metadata.
     """
 
-    _kernel: str | None = None
+    try:
+        kernel = d["kernel"]
+    except KeyError:
+        msg = "no 'kernel' in mapping"
+        raise ValueError(msg) from None
 
-    @property
-    def kernel(cls) -> str | None:
-        return cls._kernel
+    try:
+        cls = _KERNELS[kernel]
+    except KeyError:
+        msg = f"unknown kernel: {kernel}"
+        raise ValueError(msg) from None
+
+    return cls.from_dict(d)
 
 
-class Mapper(metaclass=MapperMeta):
+class Mapper(metaclass=ABCMeta):
     """
     Abstract base class for mappers.
     """
 
-    def __init_subclass__(cls, /, kernel: str, **kwargs):
+    __kernel: str
+
+    def __init_subclass__(cls, /, kernel: str, **kwargs) -> None:
         """
         Initialise mapper subclasses with a *kernel* parameter.
         """
         super().__init_subclass__(**kwargs)
-        cls._kernel = kernel
+        cls.__kernel = kernel
         _KERNELS[kernel] = cls
 
+    @classmethod
+    @abstractmethod
+    def from_dict(cls, d: Mapping[str, Any]) -> Self:
+        """
+        Create a new mapper instance from a dictionary of parameters.
+        """
+
     def __init__(self) -> None:
         """
         Initialise a new mapper instance.
         """
         self._metadata: dict[str, Any] = {
-            "kernel": self.__class__.kernel,
+            "kernel": self.__kernel,
         }
 
+    @property
+    def kernel(self) -> str:
+        """
+        Return the name of the kernel for this mapper.
+        """
+        return self.__kernel
+
     @property
     def metadata(self) -> Mapping[str, Any]:
         """
@@ -128,7 +153,13 @@ def transform(
         """
 
     @abstractmethod
-    def deconvolve(self, alm: ArrayLike, *, inplace: bool = False) -> ArrayLike:
+    def kl(self, lmax: int, spin: int = 0) -> NDArray[Any]:
+        """
+        Return the convolution kernel in harmonic space.
+        """
+
+    def bl(self, lmax: int, spin: int = 0) -> None | NDArray[Any]:
         """
-        Remove this mapper's convolution kernel from *alm*.
+        Return the biasing kernel in harmonic space.
         """
+        return None