Merge branch 'main' into saransh/array-api-test

.github/test-constraints.txt

@@ -1,9 +1,6 @@
 numpy
 --only-binary numpy
 
-scipy
---only-binary scipy
-
 healpy
 --only-binary healpy
 

.pre-commit-config.yaml

@@ -33,7 +33,7 @@ repos:
         args:
           - --strict
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: v0.8.0
+    rev: v0.8.6
     hooks:
       - id: ruff
       - id: ruff-format
@@ -48,7 +48,7 @@ repos:
     hooks:
       - id: toml-sort-fix
   - repo: https://github.com/rbubley/mirrors-prettier
-    rev: v3.3.3
+    rev: v3.4.2
     hooks:
       - id: prettier
   - repo: https://github.com/igorshubovych/markdownlint-cli
@@ -62,11 +62,11 @@ repos:
     hooks:
       - id: forbid-tabs
   - repo: https://github.com/rhysd/actionlint
-    rev: v1.7.4
+    rev: v1.7.6
     hooks:
       - id: actionlint
   - repo: https://github.com/pre-commit/mirrors-mypy
-    rev: v1.13.0
+    rev: v1.14.1
     hooks:
       - id: mypy
         files: ^glass/

glass/core/array.py

@@ -166,7 +166,7 @@ def trapezoid_product(
     y = np.interp(x, *f)
     for f_ in ff:
         y *= np.interp(x, *f_)
-    return np.trapezoid(y, x, axis=axis)  # type: ignore[no-any-return]
+    return np.atleast_1d(np.trapezoid(y, x, axis=axis))
 
 
 def cumulative_trapezoid(

glass/points.py

@@ -85,7 +85,7 @@ def effective_bias(
 
     """
     norm = np.trapezoid(w.wa, w.za)
-    return trapezoid_product((z, bz), (w.za, w.wa)) / norm
+    return (trapezoid_product((z, bz), (w.za, w.wa)) / norm).astype(np.float64)
 
 
 def linear_bias(
@@ -107,7 +107,7 @@ def linear_bias(
         The density contrast after biasing.
 
     """
-    return b * delta
+    return (b * delta).astype(np.float64)
 
 
 def loglinear_bias(
@@ -413,6 +413,6 @@ def position_weights(
     densities = densities / np.sum(densities, axis=0)
     # apply bias after normalisation
     if bias is not None:
-        densities = densities * bias
+        densities = (densities * bias).astype(np.float64)
     # densities now contains the relative contribution with bias applied
     return densities

glass/shells.py

@@ -235,8 +235,8 @@ def tophat_windows(
     ws = []
     for zmin, zmax in itertools.pairwise(zbins):
         n = max(round((zmax - zmin) / dz), 2)
-        z = np.linspace(zmin, zmax, n)
-        w = wht(z)
+        z = np.linspace(zmin, zmax, n, dtype=np.float64)
+        w = np.asarray(wht(z), dtype=np.float64)
         zeff = np.trapezoid(w * z, z) / np.trapezoid(w, z)
         ws.append(RadialWindow(z, w, zeff))
     return ws
@@ -410,7 +410,7 @@ def restrict(
     z_ = np.compress(np.greater(z, w.za[0]) & np.less(z, w.za[-1]), z)
     zr = np.union1d(w.za, z_)
     fr = ndinterp(zr, z, f, left=0.0, right=0.0) * ndinterp(zr, w.za, w.wa)
-    return zr, fr
+    return zr.astype(np.float64), fr.astype(np.float64)
 
 
 def partition(
@@ -565,7 +565,7 @@ def partition_lstsq(
     # compute the union of all given redshift grids
     zp = z
     for w in shells:
-        zp = np.union1d(zp, w.za)
+        zp = np.union1d(zp, w.za).astype(np.float64)
 
     # get extra leading axes of fz
     *dims, _ = np.shape(fz)
@@ -629,7 +629,7 @@ def partition_nnls(
     # compute the union of all given redshift grids
     zp = z
     for w in shells:
-        zp = np.union1d(zp, w.za)
+        zp = np.union1d(zp, w.za).astype(np.float64)
 
     # get extra leading axes of fz
     *dims, _ = np.shape(fz)
@@ -673,7 +673,8 @@ def partition_nnls(
     # for each dim, find non-negative weights x such that y == r @ x
     x = np.empty([len(shells), *dims])
     for i in np.ndindex(*dims):
-        x[(..., *i)] = nnls(r, y[i])
+        index = (slice(None),) + i  # noqa: RUF005
+        x[index] = nnls(r, y[i])
 
     # all done
     return x
@@ -740,9 +741,9 @@ def _uniform_grid(
 
     """
     if step is not None and num is None:
-        return np.arange(start, np.nextafter(stop + step, stop), step)
+        return np.arange(start, np.nextafter(stop + step, stop), step, dtype=np.float64)
     if step is None and num is not None:
-        return np.linspace(start, stop, num + 1)
+        return np.linspace(start, stop, num + 1, dtype=np.float64)
     msg = "exactly one of grid step size or number of steps must be given"
     raise ValueError(msg)
 

pyproject.toml

@@ -61,7 +61,6 @@ test = [
     "pytest-doctestplus",
     "pytest-mock",
     "pytest-rerunfailures",
-    "scipy",
 ]
 
 [project.urls]

tests/core/test_algorithm.py

@@ -1,33 +1,33 @@
-import importlib.util
-
 import numpy as np
 import pytest
 
-from glass.core.algorithm import nnls as nnls_glass
-
-# check if scipy is available for testing
-HAVE_SCIPY = importlib.util.find_spec("scipy") is not None
+from glass.core.algorithm import nnls
 
 
-@pytest.mark.skipif(not HAVE_SCIPY, reason="test requires SciPy")
 def test_nnls(rng: np.random.Generator) -> None:
-    from scipy.optimize import nnls as nnls_scipy
-
-    # cross-check output with scipy's nnls
+    # check output
 
-    a = rng.standard_normal((100, 20))
-    b = rng.standard_normal((100,))
+    a = np.arange(25.0).reshape(-1, 5)
+    b = np.arange(5.0)
+    y = a @ b
+    res = nnls(a, y)
+    assert np.linalg.norm((a @ res) - y) < 1e-7
 
-    x_glass = nnls_glass(a, b)
-    x_scipy, _ = nnls_scipy(a, b)
-
-    np.testing.assert_allclose(x_glass, x_scipy)
+    a = rng.uniform(low=-10, high=10, size=[50, 10])
+    b = np.abs(rng.uniform(low=-2, high=2, size=[10]))
+    b[::2] = 0
+    x = a @ b
+    res = nnls(a, x, tol=500 * np.linalg.norm(a, 1) * np.spacing(1.0))
+    np.testing.assert_allclose(res, b, rtol=0.0, atol=1e-10)
 
     # check matrix and vector's shape
 
+    a = rng.standard_normal((100, 20))
+    b = rng.standard_normal((100,))
+
     with pytest.raises(ValueError, match="input `a` is not a matrix"):
-        nnls_glass(b, a)
+        nnls(b, a)
     with pytest.raises(ValueError, match="input `b` is not a vector"):
-        nnls_glass(a, a)
+        nnls(a, a)
     with pytest.raises(ValueError, match="the shapes of `a` and `b` do not match"):
-        nnls_glass(a.T, b)
+        nnls(a.T, b)

tests/core/test_array.py

@@ -1,5 +1,3 @@
-import importlib.util
-
 import numpy as np
 import numpy.typing as npt
 import pytest
@@ -12,9 +10,6 @@
     trapezoid_product,
 )
 
-# check if scipy is available for testing
-HAVE_SCIPY = importlib.util.find_spec("scipy") is not None
-
 
 def test_broadcast_first() -> None:
     a = np.ones((2, 3, 4))
@@ -156,50 +151,47 @@ def test_trapezoid_product() -> None:
     np.testing.assert_allclose(s, 1.0)
 
 
-@pytest.mark.skipif(not HAVE_SCIPY, reason="test requires SciPy")
 def test_cumulative_trapezoid() -> None:
-    import scipy.integrate as spi
-
     # 1D f and x
 
     f = np.array([1, 2, 3, 4])
     x = np.array([0, 1, 2, 3])
 
-    # default dtype (int - not supported by scipy)
+    # default dtype (int)
 
-    glass_ct = cumulative_trapezoid(f, x)
-    np.testing.assert_allclose(glass_ct, np.array([0, 1, 4, 7]))
+    ct = cumulative_trapezoid(f, x)
+    np.testing.assert_allclose(ct, np.array([0, 1, 4, 7]))
 
     # explicit dtype (float)
 
-    glass_ct = cumulative_trapezoid(f, x, dtype=float)
-    scipy_ct = spi.cumulative_trapezoid(f, x, initial=0)
-    np.testing.assert_allclose(glass_ct, scipy_ct)
+    ct = cumulative_trapezoid(f, x, dtype=float)
+    np.testing.assert_allclose(ct, np.array([0.0, 1.5, 4.0, 7.5]))
 
     # explicit return array
 
-    result = cumulative_trapezoid(f, x, dtype=float, out=np.zeros((4,)))
-    scipy_ct = spi.cumulative_trapezoid(f, x, initial=0)
-    np.testing.assert_allclose(result, scipy_ct)
+    out = np.zeros((4,))
+    ct = cumulative_trapezoid(f, x, dtype=float, out=out)
+    np.testing.assert_equal(ct, out)
 
     # 2D f and 1D x
 
     f = np.array([[1, 4, 9, 16], [2, 3, 5, 7]])
     x = np.array([0, 1, 2.5, 4])
 
-    # default dtype (int - not supported by scipy)
+    # default dtype (int)
 
-    glass_ct = cumulative_trapezoid(f, x)
-    np.testing.assert_allclose(glass_ct, np.array([[0, 2, 12, 31], [0, 2, 8, 17]]))
+    ct = cumulative_trapezoid(f, x)
+    np.testing.assert_allclose(ct, np.array([[0, 2, 12, 31], [0, 2, 8, 17]]))
 
     # explicit dtype (float)
 
-    glass_ct = cumulative_trapezoid(f, x, dtype=float)
-    scipy_ct = spi.cumulative_trapezoid(f, x, initial=0)
-    np.testing.assert_allclose(glass_ct, scipy_ct)
+    ct = cumulative_trapezoid(f, x, dtype=float)
+    np.testing.assert_allclose(
+        ct, np.array([[0.0, 2.5, 12.25, 31.0], [0.0, 2.5, 8.5, 17.5]])
+    )
 
     # explicit return array
 
-    glass_ct = cumulative_trapezoid(f, x, dtype=float, out=np.zeros((2, 4)))
-    scipy_ct = spi.cumulative_trapezoid(f, x, initial=0)
-    np.testing.assert_allclose(glass_ct, scipy_ct)
+    out = np.zeros((2, 4))
+    ct = cumulative_trapezoid(f, x, dtype=float, out=out)
+    np.testing.assert_equal(ct, out)