gh-280: add mypy to codebase (#359)

Adding `mypy` to a pre-existing codebase is never easy. I initially
attempted this in #308, but have since split this up into separate
issues:
- [x] #347
- [x] #356
- [ ] #358

In this PR, `mypy` is added but in order for CI to pass, the `# type:
ignore[]` syntax is used throughout. I didn't want to tackle these here
too (see #308) as it gets quite messy.

One thing I have done (following #356) is change every empty
`npt.NDArray` to `npt.NDArray[typing.Any]` (see #330), as this actually
results in fewer errors than leaving them all blank. Ideally, we'd like
to fill in as many of the `typing.Any` as possible - they're a bit
useless by themselves. However, that is not the priority for now. Plus,
I expect typing to change when #67 is tackled.

.pre-commit-config.yaml

@@ -65,3 +65,9 @@ repos:
     rev: v1.7.3
     hooks:
       - id: actionlint
+  - repo: https://github.com/pre-commit/mirrors-mypy
+    rev: v1.11.2
+    hooks:
+      - id: mypy
+        additional_dependencies:
+          - numpy

docs/conf.py

@@ -54,7 +54,7 @@
 
 html_logo = "_static/logo.png"
 html_favicon = "_static/favicon.ico"
-html_css_files = []
+html_css_files = []  # type: ignore[var-annotated]
 
 
 # -- Intersphinx -------------------------------------------------------------

glass/core/array.py

@@ -5,14 +5,14 @@
 import numpy as np
 
 
-def broadcast_first(*arrays):
+def broadcast_first(*arrays):  # type: ignore[no-untyped-def]
     """Broadcast arrays, treating the first axis as common."""
     arrays = tuple(np.moveaxis(a, 0, -1) if np.ndim(a) else a for a in arrays)
     arrays = np.broadcast_arrays(*arrays)
     return tuple(np.moveaxis(a, -1, 0) if np.ndim(a) else a for a in arrays)
 
 
-def broadcast_leading_axes(*args):
+def broadcast_leading_axes(*args):  # type: ignore[no-untyped-def]
     """
     Broadcast all but the last N axes.
 
@@ -49,7 +49,7 @@ def broadcast_leading_axes(*args):
     return (dims, *arrs)
 
 
-def ndinterp(x, xp, fp, axis=-1, left=None, right=None, period=None):  # noqa: PLR0913
+def ndinterp(x, xp, fp, axis=-1, left=None, right=None, period=None):  # type: ignore[no-untyped-def] # noqa: PLR0913
     """Interpolate multi-dimensional array over axis."""
     return np.apply_along_axis(
         partial(np.interp, x, xp),
@@ -61,7 +61,7 @@ def ndinterp(x, xp, fp, axis=-1, left=None, right=None, period=None):  # noqa: P
     )
 
 
-def trapz_product(f, *ff, axis=-1):
+def trapz_product(f, *ff, axis=-1):  # type: ignore[no-untyped-def]
     """Trapezoidal rule for a product of functions."""
     x, _ = f
     for x_, _ in ff:
@@ -72,10 +72,10 @@ def trapz_product(f, *ff, axis=-1):
     y = np.interp(x, *f)
     for f_ in ff:
         y *= np.interp(x, *f_)
-    return np.trapz(y, x, axis=axis)
+    return np.trapz(y, x, axis=axis)  # type: ignore[attr-defined]
 
 
-def cumtrapz(f, x, dtype=None, out=None):
+def cumtrapz(f, x, dtype=None, out=None):  # type: ignore[no-untyped-def]
     """Cumulative trapezoidal rule along last axis."""
     if out is None:
         out = np.empty_like(f, dtype=dtype)

glass/ext/__init__.py

@@ -7,7 +7,7 @@
 """
 
 
-def _extend_path(path, name) -> list:
+def _extend_path(path, name) -> list:  # type: ignore[no-untyped-def, type-arg]
     import os.path
     from pkgutil import extend_path
 

glass/fields.py

@@ -37,17 +37,19 @@
 
 # types
 Size = typing.Optional[typing.Union[int, tuple[int, ...]]]
-Iternorm = tuple[typing.Optional[int], npt.NDArray, npt.NDArray]
-ClTransform = typing.Union[str, typing.Callable[[npt.NDArray], npt.NDArray]]
+Iternorm = tuple[typing.Optional[int], npt.NDArray[typing.Any], npt.NDArray[typing.Any]]
+ClTransform = typing.Union[
+    str, typing.Callable[[npt.NDArray[typing.Any]], npt.NDArray[typing.Any]]
+]
 Cls = collections.abc.Sequence[
-    typing.Union[npt.NDArray, collections.abc.Sequence[float]]
+    typing.Union[npt.NDArray[typing.Any], collections.abc.Sequence[float]]
 ]
-Alms = npt.NDArray
+Alms = npt.NDArray[typing.Any]
 
 
 def iternorm(
     k: int,
-    cov: collections.abc.Iterable[npt.NDArray],
+    cov: collections.abc.Iterable[npt.NDArray[typing.Any]],
     size: Size = None,
 ) -> collections.abc.Generator[Iternorm, None, None]:
     """Return the vector a and variance sigma^2 for iterative normal sampling."""
@@ -109,15 +111,15 @@ def iternorm(
 
 def cls2cov(
     cls: Cls, nl: int, nf: int, nc: int
-) -> collections.abc.Generator[npt.NDArray, None, None]:
+) -> collections.abc.Generator[npt.NDArray[typing.Any], None, None]:
     """Return array of cls as a covariance matrix for iterative sampling."""
     cov = np.zeros((nl, nc + 1))
     end = 0
     for j in range(nf):
         begin, end = end, end + j + 1
         for i, cl in enumerate(cls[begin:end][: nc + 1]):
             if cl is None:
-                cov[:, i] = 0
+                cov[:, i] = 0  # type: ignore[unreachable]
             else:
                 if i == 0 and np.any(np.less(cl, 0)):
                     msg = "negative values in cl"
@@ -129,7 +131,7 @@ def cls2cov(
         yield cov
 
 
-def multalm(alm: Alms, bl: npt.NDArray, *, inplace: bool = False) -> Alms:
+def multalm(alm: Alms, bl: npt.NDArray[typing.Any], *, inplace: bool = False) -> Alms:
     """Multiply alm by bl."""
     n = len(bl)
     out = np.asanyarray(alm) if inplace else np.copy(alm)
@@ -142,7 +144,7 @@ def transform_cls(cls: Cls, tfm: ClTransform, pars: tuple[typing.Any, ...] = ())
     """Transform Cls to Gaussian Cls."""
     gls = []
     for cl in cls:
-        if cl is not None and len(cl) > 0:
+        if cl is not None and len(cl) > 0:  # type: ignore[redundant-expr]
             monopole = 0.0 if cl[0] == 0 else None
             gl, info, _, _ = gaussiancl(cl, tfm, pars, monopole=monopole)
             if info == 0:
@@ -187,7 +189,7 @@ def gaussian_gls(
 
     gls = []
     for cl in cls:
-        if cl is not None and len(cl) > 0:
+        if cl is not None and len(cl) > 0:  # type: ignore[redundant-expr]
             if lmax is not None:
                 cl = cl[: lmax + 1]  # noqa: PLW2901
             if nside is not None:
@@ -216,7 +218,7 @@ def generate_gaussian(
     *,
     ncorr: int | None = None,
     rng: np.random.Generator | None = None,
-) -> collections.abc.Generator[npt.NDArray, None, None]:
+) -> collections.abc.Generator[npt.NDArray[typing.Any], None, None]:
     """
     Sample Gaussian random fields from Cls iteratively.
 
@@ -253,7 +255,7 @@ def generate_gaussian(
         ncorr = ngrf - 1
 
     # number of modes
-    n = max((len(gl) for gl in gls if gl is not None), default=0)
+    n = max((len(gl) for gl in gls if gl is not None), default=0)  # type: ignore[redundant-expr]
     if n == 0:
         msg = "all gls are empty"
         raise ValueError(msg)
@@ -302,7 +304,7 @@ def generate_lognormal(
     *,
     ncorr: int | None = None,
     rng: np.random.Generator | None = None,
-) -> collections.abc.Generator[npt.NDArray, None, None]:
+) -> collections.abc.Generator[npt.NDArray[typing.Any], None, None]:
     """Sample lognormal random fields from Gaussian Cls iteratively."""
     for i, m in enumerate(generate_gaussian(gls, nside, ncorr=ncorr, rng=rng)):
         # compute the variance of the auto-correlation
@@ -324,7 +326,7 @@ def generate_lognormal(
         yield m
 
 
-def getcl(cls, i, j, lmax=None):
+def getcl(cls, i, j, lmax=None):  # type: ignore[no-untyped-def]
     """
     Return a specific angular power spectrum from an array.
 
@@ -354,7 +356,7 @@ def getcl(cls, i, j, lmax=None):
     return cl
 
 
-def effective_cls(
+def effective_cls(  # type: ignore[no-untyped-def]
     cls, weights1, weights2=None, *, lmax=None
 ) -> npt.NDArray[np.float64]:
     r"""
@@ -408,7 +410,7 @@ def effective_cls(
     if weights2 is weights1:
         pairs = combinations_with_replacement(np.ndindex(shape1[1:]), 2)
     else:
-        pairs = product(np.ndindex(shape1[1:]), np.ndindex(shape2[1:]))
+        pairs = product(np.ndindex(shape1[1:]), np.ndindex(shape2[1:]))  # type: ignore[assignment]
 
     # create the output array: axes for all input axes plus lmax+1
     out = np.empty(shape1[1:] + shape2[1:] + (lmax + 1,))
@@ -421,7 +423,7 @@ def effective_cls(
     for j1, j2 in pairs:
         w1, w2 = weights1[c + j1], weights2[c + j2]
         cl = sum(
-            w1[i1] * w2[i2] * getcl(cls, i1, i2, lmax=lmax)
+            w1[i1] * w2[i2] * getcl(cls, i1, i2, lmax=lmax)  # type: ignore[no-untyped-call]
             for i1, i2 in np.ndindex(n, n)
         )
         out[j1 + j2] = cl

glass/galaxies.py

@@ -37,7 +37,7 @@ def redshifts(
     w: RadialWindow,
     *,
     rng: np.random.Generator | None = None,
-) -> npt.NDArray:
+) -> npt.NDArray[typing.Any]:
     """
     Sample redshifts from a radial window function.
 
@@ -67,7 +67,7 @@ def redshifts_from_nz(
     *,
     rng: np.random.Generator | None = None,
     warn: bool = True,
-) -> npt.NDArray:
+) -> npt.NDArray[typing.Any]:
     """
     Generate galaxy redshifts from a source distribution.
 
@@ -111,43 +111,43 @@ def redshifts_from_nz(
         rng = np.random.default_rng()
 
     # bring inputs' leading axes into common shape
-    dims, count, z, nz = broadcast_leading_axes((count, 0), (z, 1), (nz, 1))
+    dims, count, z, nz = broadcast_leading_axes((count, 0), (z, 1), (nz, 1))  # type: ignore[no-untyped-call]
 
     # list of results for all dimensions
-    redshifts = np.empty(count.sum())
+    redshifts = np.empty(count.sum())  # type: ignore[union-attr]
 
     # keep track of the number of sampled redshifts
     total = 0
 
     # go through extra dimensions; also works if dims is empty
     for k in np.ndindex(dims):
         # compute the CDF of each galaxy population
-        cdf = cumtrapz(nz[k], z[k], dtype=float)
+        cdf = cumtrapz(nz[k], z[k], dtype=float)  # type: ignore[call-overload, index, no-untyped-call]
         cdf /= cdf[-1]
 
         # sample redshifts and store result
-        redshifts[total : total + count[k]] = np.interp(
-            rng.uniform(0, 1, size=count[k]),
+        redshifts[total : total + count[k]] = np.interp(  # type: ignore[call-overload, index, misc, operator]
+            rng.uniform(0, 1, size=count[k]),  # type: ignore[arg-type, call-overload, index]
             cdf,
-            z[k],
+            z[k],  # type: ignore[arg-type, call-overload, index]
         )
-        total += count[k]
+        total += count[k]  # type: ignore[assignment, call-overload, index, operator]
 
     assert total == redshifts.size  # noqa: S101
 
     return redshifts
 
 
 def galaxy_shear(  # noqa: PLR0913
-    lon: npt.NDArray,
-    lat: npt.NDArray,
-    eps: npt.NDArray,
-    kappa: npt.NDArray,
-    gamma1: npt.NDArray,
-    gamma2: npt.NDArray,
+    lon: npt.NDArray[typing.Any],
+    lat: npt.NDArray[typing.Any],
+    eps: npt.NDArray[typing.Any],
+    kappa: npt.NDArray[typing.Any],
+    gamma1: npt.NDArray[typing.Any],
+    gamma2: npt.NDArray[typing.Any],
     *,
     reduced_shear: bool = True,
-) -> npt.NDArray:
+) -> npt.NDArray[typing.Any]:
     """
     Observed galaxy shears from weak lensing.
 
@@ -212,7 +212,7 @@ def gaussian_phz(
     lower: npt.ArrayLike | None = None,
     upper: npt.ArrayLike | None = None,
     rng: np.random.Generator | None = None,
-) -> npt.NDArray:
+) -> npt.NDArray[typing.Any]:
     r"""
     Photometric redshifts assuming a Gaussian error.
 
@@ -264,26 +264,26 @@ def gaussian_phz(
     sigma = np.add(1, z) * sigma_0
     dims = np.shape(sigma)
 
-    zphot = rng.normal(z, sigma)
+    zphot = rng.normal(z, sigma)  # type: ignore[arg-type]
 
     if lower is None:
         lower = 0.0
     if upper is None:
         upper = np.inf
 
-    if not np.all(lower < upper):
+    if not np.all(lower < upper):  # type: ignore[operator]
         msg = "requires lower < upper"
         raise ValueError(msg)
 
     if not dims:
-        while zphot < lower or zphot > upper:
-            zphot = rng.normal(z, sigma)
+        while zphot < lower or zphot > upper:  # type: ignore[operator]
+            zphot = rng.normal(z, sigma)  # type: ignore[arg-type]
     else:
         z = np.broadcast_to(z, dims)
-        trunc = np.where((zphot < lower) | (zphot > upper))[0]
+        trunc = np.where((zphot < lower) | (zphot > upper))[0]  # type: ignore[operator]
         while trunc.size:
-            znew = rng.normal(z[trunc], sigma[trunc])
-            zphot[trunc] = znew
-            trunc = trunc[(znew < lower) | (znew > upper)]
+            znew = rng.normal(z[trunc], sigma[trunc])  # type: ignore[arg-type]
+            zphot[trunc] = znew  # type: ignore[index]
+            trunc = trunc[(znew < lower) | (znew > upper)]  # type: ignore[operator]
 
-    return zphot
+    return zphot  # type: ignore[return-value]