gh-163: full support of NumPy v2

.github/workflows/test.yml

@@ -27,8 +27,6 @@ jobs:
     strategy:
       matrix:
         python-version:
-          - "3.8"
-          - "3.9"
           - "3.10"
           - "3.11"
           - "3.12"

examples/1-basic/photoz.ipynb

@@ -55,7 +55,7 @@
     "dndz = n_arcmin2 * glass.smail_nz(z, 1.0, 2.2, 1.5)\n",
     "\n",
     "# compute the over galaxy number density on the sphere\n",
-    "ngal = np.trapz(dndz, z)"
+    "ngal = np.trapezoid(dndz, z)"
    ]
   },
   {
@@ -212,7 +212,7 @@
     "tomo_nz = glass.tomo_nz_gausserr(z, dndz, phz_sigma_0, zbins)\n",
     "tomo_nz *= ARCMIN2_SPHERE * (z[-1] - z[0]) / 40\n",
     "\n",
-    "for (z1, z2), nz in zip(zbins, tomo_nz):\n",
+    "for (z1, z2), nz in zip(zbins, tomo_nz, strict=False):\n",
     "    plt.hist(\n",
     "        ztrue[(z1 <= zphot) & (zphot < z2)],\n",
     "        bins=40,\n",

examples/2-advanced/cosmic_shear.ipynb

@@ -161,7 +161,7 @@
     "# localised redshift distribution with the given density\n",
     "z = np.arange(0.0, 2.0, 0.01)\n",
     "dndz = np.exp(-((z - 0.5) ** 2) / (0.1) ** 2)\n",
-    "dndz *= n_arcmin2 / np.trapz(dndz, z)\n",
+    "dndz *= n_arcmin2 / np.trapezoid(dndz, z)\n",
     "\n",
     "# distribute dN/dz over the radial window functions\n",
     "ngal = glass.partition(z, dndz, ws)"

glass/core/array.py

@@ -45,7 +45,9 @@ def broadcast_leading_axes(*args):
         shapes.append(s[:i])
         trails.append(s[i:])
     dims = np.broadcast_shapes(*shapes)
-    arrs = (np.broadcast_to(a, dims + t) for (a, _), t in zip(args, trails))
+    arrs = (
+        np.broadcast_to(a, dims + t) for (a, _), t in zip(args, trails, strict=False)
+    )
     return (dims, *arrs)
 
 
@@ -72,7 +74,7 @@ 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.trapezoid(y, x, axis=axis)
 
 
 def cumtrapz(f, x, dtype=None, out=None):

glass/fields.py

@@ -27,9 +27,10 @@
 from __future__ import annotations
 
 import warnings
+from collections.abc import Callable, Generator, Iterable, Sequence
 
 # typing
-from typing import Any, Callable, Generator, Iterable, Optional, Sequence, Tuple, Union
+from typing import Any
 
 import healpy as hp
 import numpy as np
@@ -38,10 +39,10 @@
 
 # types
 Array = np.ndarray
-Size = Union[None, int, Tuple[int, ...]]
-Iternorm = Tuple[Optional[int], Array, Array]
-ClTransform = Union[str, Callable[[Array], Array]]
-Cls = Sequence[Union[Array, Sequence[float]]]
+Size = None | int | tuple[int, ...]
+Iternorm = tuple[int | None, Array, Array]
+ClTransform = str | Callable[[Array], Array]
+Cls = Sequence[Array | Sequence[float]]
 Alms = np.ndarray
 
 

glass/lensing.py

@@ -32,13 +32,15 @@
 from __future__ import annotations
 
 # typing support
-from typing import TYPE_CHECKING, Sequence
+from typing import TYPE_CHECKING
 
 import healpy as hp
 import numpy as np
 
 if TYPE_CHECKING:
     # to prevent circular dependencies, only import these for type checking
+    from collections.abc import Sequence
+
     from numpy.typing import ArrayLike, NDArray
 
     from cosmology import Cosmology
@@ -295,7 +297,7 @@ def add_window(self, delta: np.ndarray, w: RadialWindow) -> None:
 
         """
         zsrc = w.zeff
-        lens_weight = np.trapz(w.wa, w.za) / np.interp(zsrc, w.za, w.wa)
+        lens_weight = np.trapezoid(w.wa, w.za) / np.interp(zsrc, w.za, w.wa)
 
         self.add_plane(delta, zsrc, lens_weight)
 

glass/observations.py

@@ -28,6 +28,7 @@
 
 from __future__ import annotations
 
+import itertools
 import math
 from typing import TYPE_CHECKING
 
@@ -122,7 +123,7 @@ def gaussian_nz(
     sigma = np.reshape(sigma, np.shape(sigma) + (1,) * np.ndim(z))
 
     nz = np.exp(-(((z - mean) / sigma) ** 2) / 2)
-    nz /= np.trapz(nz, z, axis=-1)[..., np.newaxis]
+    nz /= np.trapezoid(nz, z, axis=-1)[..., np.newaxis]
 
     if norm is not None:
         nz *= norm
@@ -184,7 +185,7 @@ def smail_nz(
     beta = np.asanyarray(beta)[..., np.newaxis]
 
     pz = z**alpha * np.exp(-alpha / beta * (z / z_mode) ** beta)
-    pz /= np.trapz(pz, z, axis=-1)[..., np.newaxis]
+    pz /= np.trapezoid(pz, z, axis=-1)[..., np.newaxis]
 
     if norm is not None:
         pz *= norm
@@ -228,7 +229,7 @@ def fixed_zbins(
         msg = "exactly one of nbins and dz must be given"
         raise ValueError(msg)
 
-    return list(zip(zbinedges, zbinedges[1:]))
+    return list(itertools.pairwise(zbinedges))
 
 
 def equal_dens_zbins(
@@ -263,7 +264,7 @@ def equal_dens_zbins(
     cuml_nz /= cuml_nz[[-1]]
     zbinedges = np.interp(np.linspace(0, 1, nbins + 1), cuml_nz, z)
 
-    return list(zip(zbinedges, zbinedges[1:]))
+    return list(itertools.pairwise(zbinedges))
 
 
 def tomo_nz_gausserr(

glass/points.py

@@ -68,7 +68,7 @@ def effective_bias(z, bz, w):
         \\;.
 
     """
-    norm = np.trapz(w.wa, w.za)
+    norm = np.trapezoid(w.wa, w.za)
     return trapz_product((z, bz), (w.za, w.wa)) / norm
 
 

glass/shells.py

@@ -44,10 +44,12 @@
 
 from __future__ import annotations
 
+import itertools
 import warnings
+from collections.abc import Callable, Sequence
 
 # type checking
-from typing import TYPE_CHECKING, Callable, NamedTuple, Sequence, Union
+from typing import TYPE_CHECKING, NamedTuple
 
 import numpy as np
 
@@ -60,7 +62,7 @@
 
 
 # types
-ArrayLike1D = Union[Sequence[float], np.ndarray]
+ArrayLike1D = Sequence[float] | np.ndarray
 WeightFunc = Callable[[ArrayLike1D], np.ndarray]
 
 
@@ -178,11 +180,11 @@ def tophat_windows(
     wht: WeightFunc
     wht = weight if weight is not None else np.ones_like
     ws = []
-    for zmin, zmax in zip(zbins, zbins[1:]):
+    for zmin, zmax in itertools.pairwise(zbins):
         n = max(round((zmax - zmin) / dz), 2)
         z = np.linspace(zmin, zmax, n)
         w = wht(z)
-        zeff = np.trapz(w * z, z) / np.trapz(w, z)
+        zeff = np.trapezoid(w * z, z) / np.trapezoid(w, z)
         ws.append(RadialWindow(z, w, zeff))
     return ws
 
@@ -233,7 +235,7 @@ def linear_windows(
         warnings.warn("first triangular window does not start at z=0", stacklevel=2)
 
     ws = []
-    for zmin, zmid, zmax in zip(zgrid, zgrid[1:], zgrid[2:]):
+    for zmin, zmid, zmax in zip(zgrid, zgrid[1:], zgrid[2:], strict=False):
         n = max(round((zmid - zmin) / dz), 2) - 1
         m = max(round((zmax - zmid) / dz), 2)
         z = np.concatenate(
@@ -295,7 +297,7 @@ def cubic_windows(
         warnings.warn("first cubic spline window does not start at z=0", stacklevel=2)
 
     ws = []
-    for zmin, zmid, zmax in zip(zgrid, zgrid[1:], zgrid[2:]):
+    for zmin, zmid, zmax in zip(zgrid, zgrid[1:], zgrid[2:], strict=False):
         n = max(round((zmid - zmin) / dz), 2) - 1
         m = max(round((zmax - zmid) / dz), 2)
         z = np.concatenate(
@@ -485,7 +487,7 @@ def partition_lstsq(
 
     # create the window function matrix
     a = [np.interp(zp, za, wa, left=0.0, right=0.0) for za, wa, _ in shells]
-    a /= np.trapz(a, zp, axis=-1)[..., None]
+    a /= np.trapezoid(a, zp, axis=-1)[..., None]
     a = a * dz
 
     # create the target vector of distribution values
@@ -495,7 +497,7 @@ def partition_lstsq(
     # append a constraint for the integral
     mult = 1 / sumtol
     a = np.concatenate([a, mult * np.ones((len(shells), 1))], axis=-1)
-    b = np.concatenate([b, mult * np.reshape(np.trapz(fz, z), (*dims, 1))], axis=-1)
+    b = np.concatenate([b, mult * np.reshape(np.trapezoid(fz, z), (*dims, 1))], axis=-1)
 
     # now a is a matrix of shape (len(shells), len(zp) + 1)
     # and b is a matrix of shape (*dims, len(zp) + 1)
@@ -539,7 +541,7 @@ def partition_nnls(
 
     # create the window function matrix
     a = [np.interp(zp, za, wa, left=0.0, right=0.0) for za, wa, _ in shells]
-    a /= np.trapz(a, zp, axis=-1)[..., None]
+    a /= np.trapezoid(a, zp, axis=-1)[..., None]
     a = a * dz
 
     # create the target vector of distribution values
@@ -549,7 +551,7 @@ def partition_nnls(
     # append a constraint for the integral
     mult = 1 / sumtol
     a = np.concatenate([a, mult * np.ones((len(shells), 1))], axis=-1)
-    b = np.concatenate([b, mult * np.reshape(np.trapz(fz, z), (*dims, 1))], axis=-1)
+    b = np.concatenate([b, mult * np.reshape(np.trapezoid(fz, z), (*dims, 1))], axis=-1)
 
     # now a is a matrix of shape (len(shells), len(zp) + 1)
     # and b is a matrix of shape (*dims, len(zp) + 1)
@@ -577,7 +579,7 @@ def partition_restrict(
     part = np.empty((len(shells),) + np.shape(fz)[:-1])
     for i, w in enumerate(shells):
         zr, fr = restrict(z, fz, w)
-        part[i] = np.trapz(fr, zr, axis=-1)
+        part[i] = np.trapezoid(fr, zr, axis=-1)
     return part
 
 
@@ -647,9 +649,9 @@ def combine(
         * np.interp(
             z,
             shell.za,
-            shell.wa / np.trapz(shell.wa, shell.za),
+            shell.wa / np.trapezoid(shell.wa, shell.za),
             left=0.0,
             right=0.0,
         )
-        for shell, weight in zip(shells, weights)
+        for shell, weight in zip(shells, weights, strict=False)
     )

pyproject.toml

@@ -11,8 +11,6 @@ classifiers = [
     "Operating System :: OS Independent",
     "Programming Language :: Python :: 3",
     "Programming Language :: Python :: 3 :: Only",
-    "Programming Language :: Python :: 3.8",
-    "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
@@ -24,7 +22,7 @@ dependencies = [
     "gaussiancl>=2022.10.21",
     "healpix>=2022.11.1",
     "healpy>=1.15.0",
-    "numpy>=1.20.0",
+    "numpy~=2.0",
 ]
 description = "Generator for Large Scale Structure"
 dynamic = [
@@ -35,7 +33,7 @@ maintainers = [
 ]
 name = "glass"
 readme = "README.md"
-requires-python = ">=3.8"
+requires-python = ">=3.10"
 license.file = "LICENSE"
 
 [project.optional-dependencies]

tests/test_fits.py

@@ -32,9 +32,10 @@ def test_basic_write(tmp_path):
     filename_gfits = "gfits.fits"  # what GLASS creates
     filename_tfits = "tfits.fits"  # file created on the fly to test against
 
-    with user.write_catalog(
-        tmp_path / filename_gfits, ext="CATALOG"
-    ) as out, fitsio.FITS(tmp_path / filename_tfits, "rw", clobber=True) as my_fits:
+    with (
+        user.write_catalog(tmp_path / filename_gfits, ext="CATALOG") as out,
+        fitsio.FITS(tmp_path / filename_tfits, "rw", clobber=True) as my_fits,
+    ):
         for i in range(my_max):
             array = np.arange(i, i + 1, delta)  # array of size 1/delta
             array2 = np.arange(i + 1, i + 2, delta)  # array of size 1/delta
@@ -43,9 +44,10 @@ def test_basic_write(tmp_path):
             names = ["RA", "RB"]
             _test_append(my_fits, arrays, names)
 
-    with fitsio.FITS(tmp_path / filename_gfits) as g_fits, fitsio.FITS(
-        tmp_path / filename_tfits
-    ) as t_fits:
+    with (
+        fitsio.FITS(tmp_path / filename_gfits) as g_fits,
+        fitsio.FITS(tmp_path / filename_tfits) as t_fits,
+    ):
         glass_data = g_fits[1].read()
         test_data = t_fits[1].read()
         assert glass_data["RA"].size == test_data["RA"].size

tests/test_lensing.py

@@ -101,7 +101,7 @@ def test_multi_plane_matrix(shells, cosmo, rng):
 
     deltas = rng.random((len(shells), 10))
     kappas = []
-    for shell, delta in zip(shells, deltas):
+    for shell, delta in zip(shells, deltas, strict=False):
         convergence.add_window(delta, shell)
         kappas.append(convergence.kappa.copy())
 
@@ -122,7 +122,7 @@ def test_multi_plane_weights(shells, cosmo, rng):
     deltas = rng.random((len(shells), 10))
     weights = rng.random((len(shells), 3))
     kappa = 0
-    for shell, delta, weight in zip(shells, deltas, weights):
+    for shell, delta, weight in zip(shells, deltas, weights, strict=False):
         convergence.add_window(delta, shell)
         kappa = kappa + weight[..., None] * convergence.kappa
     kappa /= weights.sum(axis=0)[..., None]

tests/test_shells.py

@@ -12,10 +12,14 @@ def test_tophat_windows():
 
     assert len(ws) == len(zb) - 1
 
-    assert all(z0 == w.za[0] and zn == w.za[-1] for w, z0, zn in zip(ws, zb, zb[1:]))
+    assert all(
+        z0 == w.za[0] and zn == w.za[-1]
+        for w, z0, zn in zip(ws, zb, zb[1:], strict=False)
+    )
 
     assert all(
-        zn <= z0 + len(w.za) * dz <= zn + dz for w, z0, zn in zip(ws, zb, zb[1:])
+        zn <= z0 + len(w.za) * dz <= zn + dz
+        for w, z0, zn in zip(ws, zb, zb[1:], strict=False)
     )
 
     assert all(np.all(w.wa == 1) for w in ws)
@@ -38,12 +42,12 @@ def test_restrict():
 
     assert fr[0] == fr[-1] == 0.0
 
-    for zi, wi in zip(w.za, w.wa):
+    for zi, wi in zip(w.za, w.wa, strict=False):
         i = np.searchsorted(zr, zi)
         assert zr[i] == zi
         assert fr[i] == wi * np.interp(zi, z, f)
 
-    for zi, fi in zip(z, f):
+    for zi, fi in zip(z, f, strict=False):
         if w.za[0] <= zi <= w.za[-1]:
             i = np.searchsorted(zr, zi)
             assert zr[i] == zi
@@ -75,4 +79,4 @@ def test_partition(method):
 
     assert part.shape == (len(shells), 3, 2)
 
-    assert np.allclose(part.sum(axis=0), np.trapz(fz, z))
+    assert np.allclose(part.sum(axis=0), np.trapezoid(fz, z))