diff --git a/pyproject.toml b/pyproject.toml
index 1cd97274..8a4a8b2e 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -14,6 +14,7 @@ dev = [
     "coverage[toml]",
     "hatch-vcs",
     "hatchling",
+    "jupyter",
     "nox",
     "nox-uv",
     {include-group = "build"},
diff --git a/tests/benchmarks/test_arraytools.py b/tests/benchmarks/test_arraytools.py
index da07fbe8..63409f37 100644
--- a/tests/benchmarks/test_arraytools.py
+++ b/tests/benchmarks/test_arraytools.py
@@ -4,6 +4,8 @@
 
 import pytest
 
+import array_api_extra as xpx
+
 import glass.arraytools
 
 if TYPE_CHECKING:
@@ -11,8 +13,6 @@
 
     from pytest_benchmark.fixture import BenchmarkFixture
 
-    from tests.fixtures.helper_classes import Compare
-
 
 @pytest.mark.unstable
 def test_broadcast_leading_axes(
@@ -44,7 +44,6 @@ def test_broadcast_leading_axes(
 @pytest.mark.unstable
 def test_cumulative_trapezoid_1d(
     benchmark: BenchmarkFixture,
-    compare: type[Compare],
     xpb: ModuleType,
 ) -> None:
     """Benchmark test for glass.arraytools.cumulative_trapezoid."""
@@ -56,7 +55,7 @@ def test_cumulative_trapezoid_1d(
     ct = benchmark(glass.arraytools.cumulative_trapezoid, f, x)
 
     # Compare to int64 as old versions of glass round to int64 if `dtype` is not passed.
-    compare.assert_allclose(
+    xpx.testing.assert_close(
         xpb.asarray(ct[:4], dtype=xpb.int64),
         xpb.asarray([0, 1, 4, 7]),
     )
@@ -66,7 +65,6 @@ def test_cumulative_trapezoid_1d(
 @pytest.mark.unstable
 def test_cumulative_trapezoid_2d(
     benchmark: BenchmarkFixture,
-    compare: type[Compare],
     xpb: ModuleType,
 ) -> None:
     """Benchmark test for glass.arraytools.cumulative_trapezoid."""
@@ -80,16 +78,16 @@ def test_cumulative_trapezoid_2d(
     )
     x = xpb.arange(scaled_length)  # [0, 1, 2, 3,...]
 
-    expected_first_4_out = [0, 1, 4, 7]
-
     ct = benchmark(glass.arraytools.cumulative_trapezoid, f, x)
 
+    expected_first_4_out = xpb.asarray([0, 1, 4, 7])
+
     # Compare to int64 as old versions of glass round to int64 if `dtype` is not passed.
-    compare.assert_allclose(
+    xpx.testing.assert_close(
         xpb.asarray(ct[0, :4], dtype=xpb.int64),
         expected_first_4_out,
     )
-    compare.assert_allclose(
+    xpx.testing.assert_close(
         xpb.asarray(ct[1, :4], dtype=xpb.int64),
         expected_first_4_out,
     )
diff --git a/tests/benchmarks/test_fields.py b/tests/benchmarks/test_fields.py
index 72a8a092..18d429ea 100644
--- a/tests/benchmarks/test_fields.py
+++ b/tests/benchmarks/test_fields.py
@@ -5,6 +5,8 @@
 
 import pytest
 
+import array_api_extra as xpx
+
 import glass
 import glass.fields
 import glass.healpix as hp
@@ -16,7 +18,7 @@
     from pytest_benchmark.fixture import BenchmarkFixture
 
     from glass._types import AngularPowerSpectra, UnifiedGenerator
-    from tests.fixtures.helper_classes import Compare, GeneratorConsumer
+    from tests.fixtures.helper_classes import GeneratorConsumer
 
 
 @pytest.mark.stable
@@ -102,7 +104,6 @@ def function_to_benchmark() -> list[Any]:
 @pytest.mark.stable
 def test_iternorm_k_0(
     benchmark: BenchmarkFixture,
-    compare: type[Compare],
     generator_consumer: type[GeneratorConsumer],
     xpb: ModuleType,
 ) -> None:
@@ -119,13 +120,12 @@ def function_to_benchmark() -> list[Any]:
     j, a, s = results[0]
     assert j is None
     assert a.shape == (0,)
-    compare.assert_allclose(xpb.asarray(s), 1.0)
+    xpx.testing.assert_close(s, xpb.asarray(1.0))
 
 
 @pytest.mark.stable
 def test_cls2cov(
     benchmark: BenchmarkFixture,
-    compare: type[Compare],
     generator_consumer: type[GeneratorConsumer],
     xpb: ModuleType,
 ) -> None:
@@ -148,9 +148,9 @@ def function_to_benchmark() -> list[Any]:
     assert cov.shape == (nl, nc + 1)
     assert cov.dtype == xpb.float64
 
-    compare.assert_allclose(cov[:, 0], xpb.asarray([1.0, 1.5, 2.0]))
-    compare.assert_allclose(cov[:, 1], xpb.asarray([1.5, 2.0, 2.5]))
-    compare.assert_allclose(cov[:, 2], 0)
+    xpx.testing.assert_close(cov[:, 0], xpb.asarray([1.0, 1.5, 2.0]))
+    xpx.testing.assert_close(cov[:, 1], xpb.asarray([1.5, 2.0, 2.5]))
+    xpx.testing.assert_close(cov[:, 2], xpb.asarray(0.0), check_shape=False)
 
 
 @pytest.mark.stable
@@ -185,7 +185,6 @@ def function_to_benchmark() -> list[Any]:
 @pytest.mark.parametrize("ncorr", [None, 1])
 def test_generate(
     benchmark: BenchmarkFixture,
-    compare: type[Compare],
     generator_consumer: type[GeneratorConsumer],
     xpb: ModuleType,
     ncorr: int | None,
@@ -214,13 +213,12 @@ def function_to_benchmark() -> list[Any]:
 
     for field in result:
         assert field.shape == (hp.nside2npix(nside),)
-    compare.assert_allclose(result[1], result[0] ** 2, atol=1e-05)
+    xpx.testing.assert_close(result[1], result[0] ** 2, atol=1e-05)
 
 
 @pytest.mark.unstable
 def test_getcl_lmax_0(
     benchmark: BenchmarkFixture,
-    compare: type[Compare],
     xpb: ModuleType,
 ) -> None:
     """Benchmarks for glass.getcl with lmax of 0."""
@@ -245,13 +243,12 @@ def test_getcl_lmax_0(
     )
     expected = xpb.asarray([max(random_i, random_j)], dtype=xpb.float64)
     assert result.shape[0] == 1
-    compare.assert_allclose(result, expected)
+    xpx.testing.assert_close(result, expected)
 
 
 @pytest.mark.unstable
 def test_getcl_lmax_larger_than_cls(
     benchmark: BenchmarkFixture,
-    compare: type[Compare],
     xpb: ModuleType,
 ) -> None:
     """Benchmarks for glass.getcl with lmax larger than the length of cl."""
@@ -277,4 +274,4 @@ def test_getcl_lmax_larger_than_cls(
     )
     expected = xpb.zeros((lmax - 1,), dtype=xpb.float64)
     assert result.shape[0] == lmax + 1
-    compare.assert_allclose(result[2:], expected)
+    xpx.testing.assert_close(result[2:], expected)
diff --git a/tests/benchmarks/test_harmonics.py b/tests/benchmarks/test_harmonics.py
index c004c36c..d0953bd6 100644
--- a/tests/benchmarks/test_harmonics.py
+++ b/tests/benchmarks/test_harmonics.py
@@ -4,6 +4,8 @@
 
 import pytest
 
+import array_api_extra as xpx
+
 glass_harmonics = pytest.importorskip(
     "glass.harmonics",
     reason="tests require glass.harmonics",
@@ -15,13 +17,10 @@
 
     from pytest_benchmark.fixture import BenchmarkFixture
 
-    from tests.fixtures.helper_classes import Compare
-
 
 @pytest.mark.unstable
 def test_multalm(
     benchmark: BenchmarkFixture,
-    compare: type[Compare],
     xpb: ModuleType,
 ) -> None:
     """Benchmarks for glass.harmonics.multalm."""
@@ -32,7 +31,7 @@ def test_multalm(
 
     result = benchmark(glass_harmonics.multalm, alm, bl)
 
-    compare.assert_allclose(
+    xpx.testing.assert_close(
         result[:5],
         xpb.asarray([scale_factor * x for x in [0.0, 3.0, 6.0, 9.0, 12.0]]),
     )
diff --git a/tests/benchmarks/test_lensing.py b/tests/benchmarks/test_lensing.py
index 94e242f5..1684ae4d 100644
--- a/tests/benchmarks/test_lensing.py
+++ b/tests/benchmarks/test_lensing.py
@@ -4,6 +4,8 @@
 
 import pytest
 
+import array_api_extra as xpx
+
 import glass
 
 if TYPE_CHECKING:
@@ -14,13 +16,11 @@
 
     from glass._types import FloatArray, UnifiedGenerator
     from glass.cosmology import Cosmology
-    from tests.fixtures.helper_classes import Compare
 
 
 @pytest.mark.stable
 def test_multi_plane_matrix(
     benchmark: BenchmarkFixture,
-    compare: type[Compare],
     cosmo: Cosmology,
     urngb: UnifiedGenerator,
     xpb: ModuleType,
@@ -38,8 +38,8 @@ def test_multi_plane_matrix(
     mat = glass.multi_plane_matrix(shells, cosmo)
     deltas = urngb.random((len(shells), 10))
 
-    compare.assert_array_equal(mat, xpb.tril(mat))
-    compare.assert_array_equal(xpb.triu(mat, k=1), 0)
+    xpx.testing.assert_equal(mat, xpb.tril(mat))
+    xpx.testing.assert_equal(xpb.triu(mat, k=1), xpb.asarray(0.0), check_shape=False)
 
     def setup_shells_and_deltas() -> tuple[
         tuple[
@@ -79,7 +79,6 @@ def multi_plane_matrix_add_window(
 @pytest.mark.stable
 def test_multi_plane_weights(
     benchmark: BenchmarkFixture,
-    compare: type[Compare],
     cosmo: Cosmology,
     urngb: UnifiedGenerator,
     xpb: ModuleType,
@@ -100,8 +99,8 @@ def test_multi_plane_weights(
 
     w_out = glass.multi_plane_weights(w_in, shells, cosmo)
 
-    compare.assert_array_equal(w_out, xpb.triu(w_out, 1))
-    compare.assert_array_equal(xpb.tril(w_out), 0)
+    xpx.testing.assert_equal(w_out, xpb.triu(w_out, 1))
+    xpx.testing.assert_equal(xpb.tril(w_out), xpb.asarray(0.0), check_shape=False)
 
     def setup_shells_deltas_and_weights() -> tuple[
         tuple[
diff --git a/tests/benchmarks/test_shells.py b/tests/benchmarks/test_shells.py
index 10d7938e..49ec8982 100644
--- a/tests/benchmarks/test_shells.py
+++ b/tests/benchmarks/test_shells.py
@@ -4,6 +4,8 @@
 
 import pytest
 
+import array_api_extra as xpx
+
 import glass
 
 if TYPE_CHECKING:
@@ -11,8 +13,6 @@
 
     from pytest_benchmark.fixture import BenchmarkFixture
 
-    from tests.fixtures.helper_classes import Compare
-
 
 @pytest.mark.unstable
 def test_radialwindow(
@@ -23,11 +23,11 @@ def test_radialwindow(
     """Benchmark for shells.RadialWindow."""
     # check zeff is computed when not provided
     arr_length = 100_000
-    expected_zeff = 66_666.0
+    expected_zeff = xpb.asarray(66_666.0)
 
     wa = xpb.arange(arr_length)
     za = xpb.arange(arr_length)
 
     w = benchmark(glass.RadialWindow, za, wa)
 
-    compare.assert_allclose(w.zeff, expected_zeff)
+    xpx.testing.assert_close(w.zeff, expected_zeff)
diff --git a/tests/core/grf/test_solver.py b/tests/core/grf/test_solver.py
index d5283b60..d3f3e8ce 100644
--- a/tests/core/grf/test_solver.py
+++ b/tests/core/grf/test_solver.py
@@ -5,11 +5,12 @@
 import numpy as np
 import pytest
 
+import array_api_extra as xpx
+
 import glass.grf
 
 if TYPE_CHECKING:
     from glass._types import FloatArray
-    from tests.fixtures.helper_classes import Compare
 
 
 @pytest.fixture(scope="session")
@@ -21,7 +22,6 @@ def cl() -> FloatArray:
 
 def test_one_transformation(
     cl: FloatArray,
-    compare: type[Compare],
     rng: np.random.Generator,
 ) -> None:
     lam = rng.random()
@@ -30,7 +30,7 @@ def test_one_transformation(
     gl1, _, _ = glass.grf.solve(cl, t)
     gl2, _, _ = glass.grf.solve(cl, t, t)
 
-    compare.assert_array_equal(gl1, gl2)
+    xpx.testing.assert_equal(gl1, gl2)
 
 
 def test_pad(
@@ -51,7 +51,6 @@ def test_pad(
 
 def test_initial(
     cl: FloatArray,
-    compare: type[Compare],
     rng: np.random.Generator,
 ) -> None:
     lam = rng.random()
@@ -62,24 +61,20 @@ def test_initial(
     gl1, _, _ = glass.grf.solve(cl, t)
     gl2, _, _ = glass.grf.solve(cl, t, initial=gl)
 
-    compare.assert_array_equal(gl1, gl2)
+    xpx.testing.assert_equal(gl1, gl2)
 
 
-def test_no_iterations(
-    cl: FloatArray,
-    compare: type[Compare],
-) -> None:
+def test_no_iterations(cl: FloatArray) -> None:
     t = glass.grf.Lognormal()
 
     gl1 = glass.grf.compute(cl, t)
     gl2, _, _ = glass.grf.solve(cl, t, maxiter=0)
 
-    compare.assert_array_equal(gl1, gl2)
+    xpx.testing.assert_equal(gl1, gl2)
 
 
 def test_lognormal(
     cl: FloatArray,
-    compare: type[Compare],
     rng: np.random.Generator,
 ) -> None:
     t1 = glass.grf.Lognormal()
@@ -93,17 +88,16 @@ def test_lognormal(
 
     assert info > 0
 
-    compare.assert_allclose(cl_[1 : cl.shape[0]], cl[1:], atol=0.0, rtol=cltol)
+    xpx.testing.assert_close(cl_[1 : cl.shape[0]], cl[1:], atol=0.0, rtol=cltol)
 
     gl_ = glass.grf.compute(cl_, t1, t2)
 
     assert gl[0] == gl0
-    compare.assert_allclose(gl_[1 : gl.shape[0]], gl[1:])
+    xpx.testing.assert_close(gl_[1 : gl.shape[0]], gl[1:])
 
 
 def test_monopole(
     cl: FloatArray,
-    compare: type[Compare],
     rng: np.random.Generator,
 ) -> None:
     t = glass.grf.Lognormal()
@@ -114,10 +108,10 @@ def test_monopole(
     gl, cl_out, _ = glass.grf.solve(cl, t, monopole=None, gltol=1e-8)
 
     assert gl[0] != 0.0
-    compare.assert_allclose(cl_out[0], cl[0])
+    xpx.testing.assert_close(cl_out[0], cl[0])
 
     gl, cl_out, _ = glass.grf.solve(cl, t, monopole=gl0, gltol=1e-8)
 
     assert gl[0] == gl0
     with pytest.raises(AssertionError, match="Not equal to tolerance"):
-        compare.assert_allclose(cl_out[0], cl[0])
+        xpx.testing.assert_close(cl_out[0], cl[0])
diff --git a/tests/core/grf/test_transformations.py b/tests/core/grf/test_transformations.py
index 1c87f6ba..a24ad492 100644
--- a/tests/core/grf/test_transformations.py
+++ b/tests/core/grf/test_transformations.py
@@ -2,26 +2,23 @@
 
 from typing import TYPE_CHECKING
 
+import array_api_extra as xpx
+
 import glass.grf
 
 if TYPE_CHECKING:
     from types import ModuleType
 
     from glass._types import UnifiedGenerator
-    from tests.fixtures.helper_classes import Compare
 
 
-def test_normal(
-    compare: type[Compare],
-    urng: UnifiedGenerator,
-) -> None:
+def test_normal(urng: UnifiedGenerator) -> None:
     t = glass.grf.Normal()
     x = urng.standard_normal(10)
-    compare.assert_array_equal(t(x, 1.0), x)
+    xpx.testing.assert_equal(t(x, 1.0), x)
 
 
 def test_lognormal(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -30,11 +27,10 @@ def test_lognormal(
         t = glass.grf.Lognormal(lam)
         x = urng.standard_normal(10)
         y = lam * xp.expm1(x - var / 2)
-        compare.assert_array_equal(t(x, var), y)
+        xpx.testing.assert_equal(t(x, var), y)
 
 
 def test_sqnormal(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -44,24 +40,22 @@ def test_sqnormal(
         t = glass.grf.SquaredNormal(a, lam)
         x = urng.standard_normal(10)
         y = lam * ((x - a) ** 2 - 1)
-        compare.assert_array_equal(t(x, var), y)
+        xpx.testing.assert_equal(t(x, var), y)
 
 
 def test_normal_normal(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
     t1 = glass.grf.Normal()
     t2 = glass.grf.Normal()
     x = urng.random(10)
-    compare.assert_array_equal(glass.grf.corr(t1, t2, x), x)
-    compare.assert_array_equal(glass.grf.icorr(t1, t2, x), x)
-    compare.assert_array_equal(glass.grf.dcorr(t1, t2, x), xp.ones_like(x))
+    xpx.testing.assert_equal(glass.grf.corr(t1, t2, x), x)
+    xpx.testing.assert_equal(glass.grf.icorr(t1, t2, x), x)
+    xpx.testing.assert_equal(glass.grf.dcorr(t1, t2, x), xp.ones_like(x))
 
 
 def test_lognormal_lognormal(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -75,13 +69,12 @@ def test_lognormal_lognormal(
     y = lam1 * lam2 * xp.expm1(x)
     dy = lam1 * lam2 * xp.exp(x)
 
-    compare.assert_array_equal(glass.grf.corr(t1, t2, x), y)
-    compare.assert_array_almost_equal_nulp(glass.grf.icorr(t1, t2, y), x)
-    compare.assert_array_equal(glass.grf.dcorr(t1, t2, x), dy)
+    xpx.testing.assert_equal(glass.grf.corr(t1, t2, x), y)
+    xpx.testing.assert_close_nulp(glass.grf.icorr(t1, t2, y), x)
+    xpx.testing.assert_equal(glass.grf.dcorr(t1, t2, x), dy)
 
 
 def test_lognormal_normal(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -94,13 +87,12 @@ def test_lognormal_normal(
     y = lam1 * x
     dy = lam1 * xp.ones_like(x)
 
-    compare.assert_array_equal(glass.grf.corr(t1, t2, x), y)
-    compare.assert_array_almost_equal_nulp(glass.grf.icorr(t1, t2, y), x)
-    compare.assert_array_equal(glass.grf.dcorr(t1, t2, x), dy)
+    xpx.testing.assert_equal(glass.grf.corr(t1, t2, x), y)
+    xpx.testing.assert_close_nulp(glass.grf.icorr(t1, t2, y), x)
+    xpx.testing.assert_equal(glass.grf.dcorr(t1, t2, x), dy)
 
 
 def test_sqnormal_sqnormal(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -117,6 +109,6 @@ def test_sqnormal_sqnormal(
     y = 2 * lam1 * lam2 * x * (x + 2 * a1 * a2)
     dy = 4 * lam1 * lam2 * (x + a1 * a2)
 
-    compare.assert_array_equal(glass.grf.corr(t1, t2, x), y)
-    compare.assert_array_almost_equal_nulp(glass.grf.icorr(t1, t2, y), x, nulp=8)
-    compare.assert_array_equal(glass.grf.dcorr(t1, t2, x), dy)
+    xpx.testing.assert_equal(glass.grf.corr(t1, t2, x), y)
+    xpx.testing.assert_close_nulp(glass.grf.icorr(t1, t2, y), x, nulp=8)
+    xpx.testing.assert_equal(glass.grf.dcorr(t1, t2, x), dy)
diff --git a/tests/core/test_algorithm.py b/tests/core/test_algorithm.py
index c6ec0639..83c1e433 100644
--- a/tests/core/test_algorithm.py
+++ b/tests/core/test_algorithm.py
@@ -14,11 +14,9 @@
     from pytest_mock import MockerFixture
 
     from glass._types import UnifiedGenerator
-    from tests.fixtures.helper_classes import Compare
 
 
 def test_nnls(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -40,7 +38,7 @@ def test_nnls(
         x,
         tol=500 * xp.linalg.matrix_norm(a, ord=1) * xp.finfo(xp.float64).eps,
     )
-    compare.assert_allclose(res, b, rtol=0.0, atol=1e-10)
+    xpx.testing.assert_close(res, b, rtol=0.0, atol=1e-10)
 
     # check matrix and vector's shape
 
@@ -56,7 +54,6 @@ def test_nnls(
 
 
 def test_cov_clip(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -77,13 +74,10 @@ def test_cov_clip(
 
     # make sure all eigenvalues are positive
     h = xp.max(xp.linalg.eigvalsh(a))
-    compare.assert_allclose(xp.linalg.eigvalsh(cov), h)
+    xpx.testing.assert_close(xp.linalg.eigvalsh(cov), h, check_shape=False)
 
 
-def test_nearcorr(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_nearcorr(xp: ModuleType) -> None:
     # from Higham (2002)
     a = xp.asarray(
         [
@@ -101,11 +95,11 @@ def test_nearcorr(
     )
 
     x = glass.algorithm.nearcorr(a)
-    compare.assert_allclose(x, b, atol=1e-4)
+    xpx.testing.assert_close(x, b, atol=1e-4)
 
     # explicit tolerance
     x = glass.algorithm.nearcorr(a, tol=1e-10)
-    compare.assert_allclose(x, b, atol=1e-4)
+    xpx.testing.assert_close(x, b, atol=1e-4)
 
     # no iterations
     with pytest.warns(
@@ -113,7 +107,7 @@ def test_nearcorr(
         match="Nearest correlation matrix not found in 0 iterations",
     ):
         x = glass.algorithm.nearcorr(a, niter=0)
-    compare.assert_allclose(x, a)
+    xpx.testing.assert_close(x, a)
 
     # non-square matrix should raise
     with pytest.raises(ValueError, match="non-square matrix"):
@@ -128,7 +122,6 @@ def test_nearcorr(
 
 
 def test_cov_nearest(
-    compare: type[Compare],
     mocker: MockerFixture,
     urng: UnifiedGenerator,
     xp: ModuleType,
@@ -154,7 +147,7 @@ def test_cov_nearest(
 
     # make sure nearcorr was called with correct input
     nearcorr.assert_called_once()
-    compare.assert_array_almost_equal_nulp(
+    xpx.testing.assert_close_nulp(
         nearcorr.call_args_list[0].args[0],
         xp.divide(a, norm),
     )
diff --git a/tests/core/test_arraytools.py b/tests/core/test_arraytools.py
index 49b8b8dc..841a29b8 100644
--- a/tests/core/test_arraytools.py
+++ b/tests/core/test_arraytools.py
@@ -4,13 +4,13 @@
 
 import pytest
 
+import array_api_extra as xpx
+
 import glass.arraytools
 
 if TYPE_CHECKING:
     from types import ModuleType
 
-    from tests.fixtures.helper_classes import Compare
-
 
 def test_broadcast_first(xp: ModuleType) -> None:
     a = xp.ones((2, 3, 4))
@@ -59,10 +59,7 @@ def test_broadcast_leading_axes(xp: ModuleType) -> None:
     assert c_out.shape == (3, 4, 5, 6)
 
 
-def test_ndinterp(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_ndinterp(xp: ModuleType) -> None:
     # test 1d interpolation
 
     xq = xp.asarray([0, 1, 2, 3, 4])
@@ -71,17 +68,17 @@ def test_ndinterp(
     x = 0.5
     y = glass.arraytools.ndinterp(x, xq, yq)
     assert y.shape == ()
-    compare.assert_allclose(y, 1.15, atol=1e-15)
+    xpx.testing.assert_close(y, xp.asarray(1.15), atol=1e-15)
 
     x = xp.asarray([0.5, 1.5, 2.5])
     y = glass.arraytools.ndinterp(x, xq, yq)
     assert y.shape == (3,)
-    compare.assert_allclose(y, [1.15, 1.25, 1.35], atol=1e-15)
+    xpx.testing.assert_close(y, xp.asarray([1.15, 1.25, 1.35]), atol=1e-15)
 
     x = xp.asarray([[0.5, 1.5], [2.5, 3.5]])
     y = glass.arraytools.ndinterp(x, xq, yq)
     assert y.shape == (2, 2)
-    compare.assert_allclose(y, [[1.15, 1.25], [1.35, 1.45]], atol=1e-15)
+    xpx.testing.assert_close(y, xp.asarray([[1.15, 1.25], [1.35, 1.45]]), atol=1e-15)
 
     # test n-dimensional interpolation in final axis
 
@@ -90,19 +87,23 @@ def test_ndinterp(
     x = 0.5
     y = glass.arraytools.ndinterp(x, xq, yq)
     assert y.shape == (2,)
-    compare.assert_allclose(y, [1.15, 2.15], atol=1e-15)
+    xpx.testing.assert_close(y, xp.asarray([1.15, 2.15]), atol=1e-15)
 
     x = xp.asarray([0.5, 1.5, 2.5])
     y = glass.arraytools.ndinterp(x, xq, yq)
     assert y.shape == (2, 3)
-    compare.assert_allclose(y, [[1.15, 1.25, 1.35], [2.15, 2.25, 2.35]], atol=1e-15)
+    xpx.testing.assert_close(
+        y,
+        xp.asarray([[1.15, 1.25, 1.35], [2.15, 2.25, 2.35]]),
+        atol=1e-15,
+    )
 
     x = xp.asarray([[0.5, 1.5], [2.5, 3.5]])
     y = glass.arraytools.ndinterp(x, xq, yq)
     assert y.shape == (2, 2, 2)
-    compare.assert_allclose(
+    xpx.testing.assert_close(
         y,
-        [[[1.15, 1.25], [1.35, 1.45]], [[2.15, 2.25], [2.35, 2.45]]],
+        xp.asarray([[[1.15, 1.25], [1.35, 1.45]], [[2.15, 2.25], [2.35, 2.45]]]),
         atol=1e-15,
     )
 
@@ -115,23 +116,23 @@ def test_ndinterp(
     x = 0.5
     y = glass.arraytools.ndinterp(x, xq, yq, axis=1)
     assert y.shape == (2, 1)
-    compare.assert_allclose(y, [[1.15], [2.15]], atol=1e-15)
+    xpx.testing.assert_close(y, xp.asarray([[1.15], [2.15]]), atol=1e-15)
 
     x = xp.asarray([0.5, 1.5, 2.5])
     y = glass.arraytools.ndinterp(x, xq, yq, axis=1)
     assert y.shape == (2, 3, 1)
-    compare.assert_allclose(
+    xpx.testing.assert_close(
         y,
-        [[[1.15], [1.25], [1.35]], [[2.15], [2.25], [2.35]]],
+        xp.asarray([[[1.15], [1.25], [1.35]], [[2.15], [2.25], [2.35]]]),
         atol=1e-15,
     )
 
     x = xp.asarray([[0.5, 1.5, 2.5, 3.5], [3.5, 2.5, 1.5, 0.5], [0.5, 3.5, 1.5, 2.5]])
     y = glass.arraytools.ndinterp(x, xq, yq, axis=1)
     assert y.shape == (2, 3, 4, 1)
-    compare.assert_allclose(
+    xpx.testing.assert_close(
         y,
-        [
+        xp.asarray([
             [
                 [[1.15], [1.25], [1.35], [1.45]],
                 [[1.45], [1.35], [1.25], [1.15]],
@@ -142,15 +143,12 @@ def test_ndinterp(
                 [[2.45], [2.35], [2.25], [2.15]],
                 [[2.15], [2.45], [2.25], [2.35]],
             ],
-        ],
+        ]),
         atol=1e-15,
     )
 
 
-def test_trapezoid_product(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_trapezoid_product(xp: ModuleType) -> None:
     x1 = xp.linspace(0, 2, 100)
     f1 = xp.full_like(x1, 2.0)
 
@@ -159,27 +157,24 @@ def test_trapezoid_product(
 
     s = glass.arraytools.trapezoid_product((x1, f1), (x2, f2))
 
-    compare.assert_allclose(s, 1.0)
+    xpx.testing.assert_close(s, xp.asarray(1.0))
 
 
-def test_cumulative_trapezoid(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_cumulative_trapezoid(xp: ModuleType) -> None:
     # 1D f and x
 
     f = xp.asarray([1, 2, 3, 4])
     x = xp.asarray([0, 1, 2, 3])
 
     ct = glass.arraytools.cumulative_trapezoid(f, x)
-    compare.assert_allclose(ct, xp.asarray([0.0, 1.5, 4.0, 7.5]))
+    xpx.testing.assert_close(ct, xp.asarray([0.0, 1.5, 4.0, 7.5]))
 
     # 2D f and 1D x
     f = xp.asarray([[1, 4, 9, 16], [2, 3, 5, 7]])
     x = xp.asarray([0, 1, 2.5, 4])
 
     ct = glass.arraytools.cumulative_trapezoid(f, x)
-    compare.assert_allclose(
+    xpx.testing.assert_close(
         ct,
         xp.asarray([[0.0, 2.5, 12.25, 31.0], [0.0, 2.5, 8.5, 17.5]]),
     )
diff --git a/tests/core/test_fields.py b/tests/core/test_fields.py
index 98addc79..2c07ebcd 100644
--- a/tests/core/test_fields.py
+++ b/tests/core/test_fields.py
@@ -3,8 +3,11 @@
 import importlib.util
 from typing import TYPE_CHECKING
 
+import numpy as np
 import pytest
 
+import array_api_extra as xpx
+
 import glass
 import glass.fields
 import glass.healpix as hp
@@ -16,7 +19,6 @@
     from pytest_mock import MockerFixture
 
     from glass._types import AngularPowerSpectra, UnifiedGenerator
-    from tests.fixtures.helper_classes import Compare
 
 HAVE_JAX = importlib.util.find_spec("jax") is not None
 
@@ -140,10 +142,7 @@ def test_iternorm(xp: ModuleType) -> None:
 
 
 @pytest.mark.skipif(not HAVE_JAX, reason="test requires jax")
-def test_cls2cov_jax(
-    compare: type[Compare],
-    jnp: ModuleType,
-) -> None:
+def test_cls2cov_jax(jnp: ModuleType) -> None:
     nl, nf, nc = 3, 3, 2
 
     generator = glass.cls2cov(
@@ -176,20 +175,17 @@ def test_cls2cov_jax(
     assert cov3.dtype == jnp.float64
 
     # cov1 has the expected value for the first iteration (different to cov1_copy)
-    compare.assert_allclose(cov1[:, 0], jnp.asarray([0.5, 0.25, 0.15]))
+    xpx.testing.assert_close(cov1[:, 0], jnp.asarray([0.5, 0.25, 0.15]))
 
     # The copies should not be equal
     with pytest.raises(AssertionError, match="Not equal to tolerance"):
-        compare.assert_allclose(cov1, cov2)
+        xpx.testing.assert_close(cov1, cov2)
 
     with pytest.raises(AssertionError, match="Not equal to tolerance"):
-        compare.assert_allclose(cov2, cov3)
+        xpx.testing.assert_close(cov2, cov3)
 
 
-def test_cls2cov_no_jax(
-    compare: type[Compare],
-    xpb: ModuleType,
-) -> None:
+def test_cls2cov_no_jax(xpb: ModuleType) -> None:
     # check output values and shape
 
     nl, nf, nc = 3, 2, 2
@@ -205,9 +201,9 @@ def test_cls2cov_no_jax(
     assert cov.shape == (nl, nc + 1)
     assert cov.dtype == xpb.float64
 
-    compare.assert_allclose(cov[:, 0], xpb.asarray([0.5, 0.25, 0.15]))
-    compare.assert_allclose(cov[:, 1], 0)
-    compare.assert_allclose(cov[:, 2], 0)
+    xpx.testing.assert_close(cov[:, 0], xpb.asarray([0.5, 0.25, 0.15]))
+    xpx.testing.assert_close(cov[:, 1], xpb.asarray(0.0), check_shape=False)
+    xpx.testing.assert_close(cov[:, 2], xpb.asarray(0.0), check_shape=False)
 
     # test negative value error
 
@@ -263,19 +259,19 @@ def test_cls2cov_no_jax(
     assert cov3.dtype == xpb.float64
 
     # cov1|2|3 reuse the same data, so should all equal the third result
-    compare.assert_allclose(cov1[:, 0], xpb.asarray([0.45, 0.25, 0.15]))
-    compare.assert_allclose(cov1, cov2)
-    compare.assert_allclose(cov2, cov3)
+    xpx.testing.assert_close(cov1[:, 0], xpb.asarray([0.45, 0.25, 0.15]))
+    xpx.testing.assert_close(cov1, cov2)
+    xpx.testing.assert_close(cov2, cov3)
 
     # cov1 has the expected value for the first iteration (different to cov1_copy)
-    compare.assert_allclose(cov1_copy[:, 0], xpb.asarray([0.5, 0.25, 0.15]))
+    xpx.testing.assert_close(cov1_copy[:, 0], xpb.asarray([0.5, 0.25, 0.15]))
 
     # The copies should not be equal
     with pytest.raises(AssertionError, match="Not equal to tolerance"):
-        compare.assert_allclose(cov1_copy, cov2_copy)
+        xpx.testing.assert_close(cov1_copy, cov2_copy)
 
     with pytest.raises(AssertionError, match="Not equal to tolerance"):
-        compare.assert_allclose(cov2_copy, cov3)
+        xpx.testing.assert_close(cov2_copy, cov3)
 
 
 def test_lognormal_gls(xp: ModuleType) -> None:
@@ -300,10 +296,7 @@ def test_lognormal_gls(xp: ModuleType) -> None:
     assert out[2].shape[0] == 3
 
 
-def test_discretized_cls(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_discretized_cls(xp: ModuleType) -> None:
     # empty cls
 
     result = glass.discretized_cls([])
@@ -348,13 +341,10 @@ def test_discretized_cls(
     for cl in result:
         n = min(cl.shape[0], pw.shape[0])  # ty: ignore[unresolved-attribute]
         expected = xp.ones(n) * pw[:n] ** 2
-        compare.assert_allclose(cl[:n], expected)
+        xpx.testing.assert_close(cl[:n], expected)
 
 
-def test_effective_cls(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_effective_cls(xp: ModuleType) -> None:
     # empty cls
 
     result = glass.effective_cls([], xp.asarray([]))
@@ -383,7 +373,7 @@ def test_effective_cls(
     result = glass.effective_cls(cls, weights1, lmax=5)
 
     assert result.shape == (1, 1, 6)
-    compare.assert_allclose(result[..., 6:], 0)
+    xpx.testing.assert_close(result[..., 6:], xp.asarray(0.0), check_shape=False)
 
     # check with weights1 and weights2 and weights1 is weights2
 
@@ -391,10 +381,7 @@ def test_effective_cls(
     assert result.shape == (1, 1, 15)
 
 
-def test_generate_grf(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_generate_grf(xp: ModuleType) -> None:
     gls: AngularPowerSpectra = [xp.asarray([1.0, 0.5, 0.1])]
     nside = 4
     ncorr = 1
@@ -417,7 +404,7 @@ def test_generate_grf(
 
     assert new_gaussian_fields[0].shape == (hp.nside2npix(nside),)
 
-    compare.assert_allclose(new_gaussian_fields[0], gaussian_fields[0])
+    xpx.testing.assert_close(new_gaussian_fields[0], gaussian_fields[0])
 
     with pytest.raises(ValueError, match="all gls are empty"):
         list(glass.fields._generate_grf([xp.asarray([])], nside))
@@ -435,10 +422,7 @@ def test_generate_lognormal(xp: ModuleType) -> None:
     next(result)
 
 
-def test_generate(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_generate(xp: ModuleType) -> None:
     # shape mismatch error
 
     fields = [lambda x, var: x, lambda x, var: x]  # noqa: ARG005
@@ -473,13 +457,10 @@ def test_generate(
 
     result = list(glass.generate(fields, gls, nside=nside))
 
-    compare.assert_allclose(result[1], result[0] ** 2, atol=1e-05)
+    xpx.testing.assert_close(result[1], result[0] ** 2, atol=1e-05)
 
 
-def test_getcl(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_getcl(xp: ModuleType) -> None:
     # make a mock Cls array with the index pairs as entries
     cls: AngularPowerSpectra = [
         xp.asarray([i, j], dtype=xp.float64)
@@ -491,19 +472,19 @@ def test_getcl(
         for j in range(10):
             result = glass.getcl(cls, i, j)
             expected = xp.asarray([min(i, j), max(i, j)], dtype=xp.float64)
-            compare.assert_allclose(xp.sort(result), expected)
+            xpx.testing.assert_close(xp.sort(result), expected)
 
             # check slicing
             result = glass.getcl(cls, i, j, lmax=0)
             expected = xp.asarray([max(i, j)], dtype=xp.float64)
             assert result.shape[0] == 1
-            compare.assert_allclose(result, expected)
+            xpx.testing.assert_close(result, expected)
 
             # check padding
             result = glass.getcl(cls, i, j, lmax=50)
             expected = xp.zeros((49,), dtype=xp.float64)
             assert result.shape[0] == 51
-            compare.assert_allclose(result[2:], expected)
+            xpx.testing.assert_close(result[2:], expected)
 
 
 def test_is_inv_triangle_number(not_triangle_numbers: list[int]) -> None:
@@ -524,10 +505,7 @@ def test_nfields_from_nspectra(not_triangle_numbers: list[int]) -> None:
             glass.nfields_from_nspectra(t)
 
 
-def test_enumerate_spectra(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_enumerate_spectra(xp: ModuleType) -> None:
     n = 100
     tn = n * (n + 1) // 2
 
@@ -542,27 +520,27 @@ def test_enumerate_spectra(
 
     # go through expected indices and values and compare
     for k, (i, j) in enumerate(indices):
-        compare.assert_allclose(next(it), (i, j, k))
+        assert next(it) == (i, j, k)
 
     # make sure iterator is exhausted
     with pytest.raises(StopIteration):
         next(it)
 
 
-def test_spectra_indices(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
-    compare.assert_array_equal(glass.spectra_indices(0), xp.zeros((0, 2)))
-    compare.assert_array_equal(glass.spectra_indices(0, xp=xp), xp.zeros((0, 2)))
-    compare.assert_array_equal(glass.spectra_indices(1, xp=xp), [[0, 0]])
-    compare.assert_array_equal(
+def test_spectra_indices(xp: ModuleType) -> None:
+    xpx.testing.assert_equal(glass.spectra_indices(0), np.zeros((0, 2), dtype=np.int64))
+    xpx.testing.assert_equal(
+        glass.spectra_indices(0, xp=xp),
+        xp.zeros((0, 2), dtype=xp.int64),
+    )
+    xpx.testing.assert_equal(glass.spectra_indices(1, xp=xp), xp.asarray([[0, 0]]))
+    xpx.testing.assert_equal(
         glass.spectra_indices(2, xp=xp),
-        [[0, 0], [1, 1], [1, 0]],
+        xp.asarray([[0, 0], [1, 1], [1, 0]]),
     )
-    compare.assert_array_equal(
+    xpx.testing.assert_equal(
         glass.spectra_indices(3, xp=xp),
-        [[0, 0], [1, 1], [1, 0], [2, 2], [2, 1], [2, 0]],
+        xp.asarray([[0, 0], [1, 1], [1, 0], [2, 2], [2, 1], [2, 0]]),
     )
 
 
@@ -677,44 +655,47 @@ def test_solve_gaussian_spectra(
         glass.solve_gaussian_spectra(fields, spectra[:2])
 
 
-def test_glass_to_healpix_spectra(compare: type[Compare]) -> None:
+def test_glass_to_healpix_spectra() -> None:
     inp = [11, 22, 21, 33, 32, 31, 44, 43, 42, 41]
     out = glass.glass_to_healpix_spectra(inp)
-    compare.assert_array_equal(out, [11, 22, 33, 44, 21, 32, 43, 31, 42, 41])
+    assert out == [11, 22, 33, 44, 21, 32, 43, 31, 42, 41]
 
 
-def test_healpix_to_glass_spectra(compare: type[Compare]) -> None:
+def test_healpix_to_glass_spectra() -> None:
     inp = [11, 22, 33, 44, 21, 32, 43, 31, 42, 41]
     out = glass.healpix_to_glass_spectra(inp)
-    compare.assert_array_equal(out, [11, 22, 21, 33, 32, 31, 44, 43, 42, 41])
+    assert out == [11, 22, 21, 33, 32, 31, 44, 43, 42, 41]
 
 
-def test_glass_to_healpix_alm(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_glass_to_healpix_alm(xp: ModuleType) -> None:
     inp = xp.asarray([00, 10, 11, 20, 21, 22, 30, 31, 32, 33], dtype=xp.complex128)
     out = glass.fields._glass_to_healpix_alm(inp)
-    compare.assert_array_equal(
+    xpx.testing.assert_equal(
         out,
         xp.asarray([00, 10, 20, 30, 11, 21, 31, 22, 32, 33], dtype=xp.complex128),
     )
 
 
-def test_lognormal_shift_hilbert2011(compare: type[Compare]) -> None:
+def test_lognormal_shift_hilbert2011(xp: ModuleType) -> None:
     zs = [0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0]
-    shifts = [glass.lognormal_shift_hilbert2011(z) for z in zs]
+    shifts = xp.asarray([glass.lognormal_shift_hilbert2011(z) for z in zs])
 
     # computed by hand
-    check = [0.0103031, 0.02975, 0.0538781, 0.0792, 0.103203, 0.12435, 0.142078, 0.1568]
+    check = xp.asarray([
+        0.0103031,
+        0.02975,
+        0.0538781,
+        0.0792,
+        0.103203,
+        0.12435,
+        0.142078,
+        0.1568,
+    ])
 
-    compare.assert_allclose(shifts, check, atol=1e-4, rtol=1e-4)
+    xpx.testing.assert_close(shifts, check, atol=1e-4, rtol=1e-4)
 
 
-def test_cov_from_spectra(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_cov_from_spectra(xp: ModuleType) -> None:
     spectra: AngularPowerSpectra = [
         xp.asarray(x)
         for x in [
@@ -727,9 +708,9 @@ def test_cov_from_spectra(
         ]
     ]
 
-    compare.assert_array_equal(
+    xpx.testing.assert_equal(
         glass.cov_from_spectra(spectra),
-        [
+        xp.asarray([
             [
                 [110, 210, 310],
                 [210, 220, 320],
@@ -750,12 +731,12 @@ def test_cov_from_spectra(
                 [213, 223, 323],
                 [313, 323, 333],
             ],
-        ],
+        ]),
     )
 
-    compare.assert_array_equal(
+    xpx.testing.assert_equal(
         glass.cov_from_spectra(spectra, lmax=1),
-        [
+        xp.asarray([
             [
                 [110, 210, 310],
                 [210, 220, 320],
@@ -766,12 +747,12 @@ def test_cov_from_spectra(
                 [211, 221, 321],
                 [311, 321, 331],
             ],
-        ],
+        ]),
     )
 
-    compare.assert_array_equal(
+    xpx.testing.assert_equal(
         glass.cov_from_spectra(spectra, lmax=4),
-        [
+        xp.asarray([
             [
                 [110, 210, 310],
                 [210, 220, 320],
@@ -797,7 +778,7 @@ def test_cov_from_spectra(
                 [0, 0, 0],
                 [0, 0, 0],
             ],
-        ],
+        ]),
     )
 
 
diff --git a/tests/core/test_galaxies.py b/tests/core/test_galaxies.py
index b9172cab..9e108864 100644
--- a/tests/core/test_galaxies.py
+++ b/tests/core/test_galaxies.py
@@ -4,6 +4,8 @@
 
 import pytest
 
+import array_api_extra as xpx
+
 import glass
 
 if TYPE_CHECKING:
@@ -12,7 +14,6 @@
     from pytest_mock import MockerFixture
 
     from glass._types import FloatArray, UnifiedGenerator
-    from tests.fixtures.helper_classes import Compare
 
 
 def test_redshifts(
@@ -205,7 +206,6 @@ def test_galaxy_shear(
 
 
 def test_gaussian_phz(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -218,13 +218,13 @@ def test_gaussian_phz(
 
     phz = glass.gaussian_phz(z, sigma_0)
 
-    compare.assert_array_equal(z, phz)
+    xpx.testing.assert_equal(z, phz)
 
     # test with rng
 
     phz = glass.gaussian_phz(z, sigma_0, rng=urng)
 
-    compare.assert_array_equal(z, phz)
+    xpx.testing.assert_equal(z, phz)
 
     # case: truncated normal
 
@@ -270,7 +270,7 @@ def test_gaussian_phz(
 
     assert phz.__array_namespace__() == xp
     assert phz.shape == (10,)
-    compare.assert_array_equal(z, phz)
+    xpx.testing.assert_equal(z, phz)
 
     # case: scalar redshift, array sigma_0
 
@@ -281,7 +281,7 @@ def test_gaussian_phz(
 
     assert phz.__array_namespace__() == xp
     assert phz.shape == (10,)
-    compare.assert_array_equal(z, phz)
+    xpx.testing.assert_equal(xp.broadcast_to(xp.asarray(1.0), phz.shape), phz)
 
     # case: array redshift, array sigma_0
 
@@ -292,7 +292,7 @@ def test_gaussian_phz(
 
     assert phz.__array_namespace__() == xp
     assert phz.shape == (11, 10)
-    compare.assert_array_equal(xp.broadcast_to(z, (11, 10)), phz)
+    xpx.testing.assert_equal(xp.broadcast_to(z, (11, 10)), phz)
 
     # shape mismatch
 
diff --git a/tests/core/test_harmonics.py b/tests/core/test_harmonics.py
index f74684b8..e00d2b91 100644
--- a/tests/core/test_harmonics.py
+++ b/tests/core/test_harmonics.py
@@ -4,18 +4,15 @@
 
 import pytest
 
+import array_api_extra as xpx
+
 import glass.harmonics
 
 if TYPE_CHECKING:
     from types import ModuleType
 
-    from tests.fixtures.helper_classes import Compare
-
 
-def test_multalm(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_multalm(xp: ModuleType) -> None:
     # check output values and shapes
 
     alm = xp.asarray([1.0, 2.0, 3.0, 4.0, 5.0, 6.0])
@@ -25,9 +22,9 @@ def test_multalm(
     result = glass.harmonics.multalm(alm, bl)
 
     expected_result = xp.asarray([2.0, 1.0, 1.5, 4.0, 5.0, 6.0])
-    compare.assert_allclose(result, expected_result)
+    xpx.testing.assert_close(result, expected_result)
     with pytest.raises(AssertionError, match="Not equal to tolerance"):
-        compare.assert_allclose(alm_copy, result)
+        xpx.testing.assert_close(alm_copy, result)
 
     # multiple with 1s
 
@@ -35,7 +32,7 @@ def test_multalm(
     bl = xp.ones(3)
 
     result = glass.harmonics.multalm(alm, bl)
-    compare.assert_allclose(result, alm)
+    xpx.testing.assert_close(result, alm)
 
     # multiple with 0s
 
@@ -44,7 +41,7 @@ def test_multalm(
     result = glass.harmonics.multalm(alm, bl)
 
     expected_result = xp.asarray([0.0, 2.0, 3.0, 0.0, 0.0, 0.0])
-    compare.assert_allclose(result, expected_result)
+    xpx.testing.assert_close(result, expected_result)
 
     # empty arrays
 
@@ -52,4 +49,4 @@ def test_multalm(
     bl = xp.asarray([])
 
     result = glass.harmonics.multalm(alm, bl)
-    compare.assert_allclose(result, alm)
+    xpx.testing.assert_close(result, alm)
diff --git a/tests/core/test_healpix.py b/tests/core/test_healpix.py
index 230a8a0c..116e7f4f 100644
--- a/tests/core/test_healpix.py
+++ b/tests/core/test_healpix.py
@@ -9,6 +9,8 @@
 import numpy as np
 import pytest
 
+import array_api_extra as xpx
+
 import glass.healpix as hp
 from glass import _rng
 
@@ -16,7 +18,7 @@
     from types import ModuleType
 
     from glass._types import UnifiedGenerator
-    from tests.fixtures.helper_classes import Compare, HealpixInputs
+    from tests.fixtures.helper_classes import HealpixInputs
 
 HAVE_ARRAY_API_STRICT = importlib.util.find_spec("array_api_strict") is not None
 
@@ -31,30 +33,29 @@
     ],
 )
 def test_alm2map_individual(
-    compare: type[Compare],
     healpix_inputs: type[HealpixInputs],
     pixwin: bool,  # noqa: FBT001
     pol: bool,  # noqa: FBT001
     urng: UnifiedGenerator,
+    xp: ModuleType,
 ) -> None:
     """Compare ``glass.healpix.alm2map`` against ``healpy.alm2map``."""
     alm = healpix_inputs.alm(rng=urng)
-    compare.assert_array_equal(
-        healpy.alm2map(
-            np.asarray(alm),
-            healpix_inputs.nside,
-            lmax=healpix_inputs.lmax,
-            pixwin=pixwin,
-            pol=pol,
-        ),
-        hp.alm2map(
-            alm,
-            healpix_inputs.nside,
-            lmax=healpix_inputs.lmax,
-            pixwin=pixwin,
-            pol=pol,
-        ),
+    old = healpy.alm2map(
+        np.asarray(alm),
+        healpix_inputs.nside,
+        lmax=healpix_inputs.lmax,
+        pixwin=pixwin,
+        pol=pol,
     )
+    new = hp.alm2map(
+        alm,
+        healpix_inputs.nside,
+        lmax=healpix_inputs.lmax,
+        pixwin=pixwin,
+        pol=pol,
+    )
+    xpx.testing.assert_equal(xp.asarray(old), new)
 
 
 @pytest.mark.parametrize(
@@ -67,11 +68,11 @@ def test_alm2map_individual(
     ],
 )
 def test_alm2map_sequence(
-    compare: type[Compare],
     healpix_inputs: type[HealpixInputs],
     pixwin: bool,  # noqa: FBT001
     pol: bool,  # noqa: FBT001
     urng: UnifiedGenerator,
+    xp: ModuleType,
 ) -> None:
     """Compare ``glass.healpix.alm2map`` against ``healpy.alm2map``."""
     alm = healpix_inputs.alm(rng=urng)
@@ -91,15 +92,15 @@ def test_alm2map_sequence(
         pixwin=pixwin,
         pol=pol,
     )
-    compare.assert_array_equal(old, new)
+    xpx.testing.assert_equal(xp.asarray(old), new)
 
 
 @pytest.mark.parametrize("spin", [1, 2])
 def test_alm2map_spin(
-    compare: type[Compare],
     healpix_inputs: type[HealpixInputs],
     spin: int,
     urng: UnifiedGenerator,
+    xp: ModuleType,
 ) -> None:
     """Compare ``glass.healpix.alm2map_spin`` against ``healpy.alm2map_spin``."""
     alm = healpix_inputs.alm(rng=urng)
@@ -114,21 +115,20 @@ def test_alm2map_spin(
     assert type(old) is type(new)
     assert len(old) == len(new)
     for i in range(len(old)):
-        compare.assert_array_equal(old[i], new[i])
+        xpx.testing.assert_equal(xp.asarray(old[i]), new[i])
 
 
 def test_almxfl(
-    compare: type[Compare],
     healpix_inputs: type[HealpixInputs],
     urng: UnifiedGenerator,
+    xp: ModuleType,
 ) -> None:
     """Compare ``glass.healpix.almxfl`` against ``healpy.almxfl``."""
     alm = healpix_inputs.alm(rng=urng)
     fl = healpix_inputs.fl(rng=urng)
-    compare.assert_array_equal(
-        healpy.almxfl(alm, fl),
-        hp.almxfl(alm, fl),
-    )
+    old = healpy.almxfl(alm, fl)
+    new = hp.almxfl(alm, fl)
+    xpx.testing.assert_equal(xp.asarray(old), new)
 
 
 @pytest.mark.parametrize(
@@ -139,7 +139,6 @@ def test_almxfl(
     ],
 )
 def test_ang2pix(  # noqa: PLR0913
-    compare: type[Compare],
     healpix_inputs: type[HealpixInputs],
     lonlat: bool,  # noqa: FBT001
     max_phi: float,
@@ -152,7 +151,7 @@ def test_ang2pix(  # noqa: PLR0913
     phis = healpix_inputs.latitudes(max_phi, rng=urng)
     old = healpix.ang2pix(healpix_inputs.nside, thetas, phis, lonlat=lonlat)
     new = hp.ang2pix(healpix_inputs.nside, thetas, phis, lonlat=lonlat, xp=xp)
-    compare.assert_array_equal(old, new)
+    xpx.testing.assert_equal(xp.asarray(old), new)
 
 
 @pytest.mark.parametrize(
@@ -163,7 +162,6 @@ def test_ang2pix(  # noqa: PLR0913
     ],
 )
 def test_ang2vec(  # noqa: PLR0913
-    compare: type[Compare],
     healpix_inputs: type[HealpixInputs],
     lonlat: bool,  # noqa: FBT001
     max_phi: float,
@@ -179,7 +177,7 @@ def test_ang2vec(  # noqa: PLR0913
     assert type(old) is type(new)
     assert len(old) == len(new)
     for i in range(len(old)):
-        compare.assert_array_equal(old[i], new[i])
+        xpx.testing.assert_equal(xp.asarray(old[i]), new[i])
 
 
 def test_get_nside(
@@ -201,28 +199,27 @@ def test_get_nside(
     ],
 )
 def test_map2alm_individual(
-    compare: type[Compare],
     healpix_inputs: type[HealpixInputs],
     pol: bool,  # noqa: FBT001
     urng: UnifiedGenerator,
     use_pixel_weights: bool,  # noqa: FBT001
+    xp: ModuleType,
 ) -> None:
     """Compare ``glass.healpix.map2alm`` against ``healpy.map2alm``."""
     kappa = healpix_inputs.kappa(rng=urng)
-    compare.assert_array_equal(
-        healpy.map2alm(
-            np.asarray(kappa),
-            lmax=healpix_inputs.lmax,
-            pol=pol,
-            use_pixel_weights=use_pixel_weights,
-        ),
-        hp.map2alm(
-            kappa,
-            lmax=healpix_inputs.lmax,
-            pol=pol,
-            use_pixel_weights=use_pixel_weights,
-        ),
+    old = healpy.map2alm(
+        np.asarray(kappa),
+        lmax=healpix_inputs.lmax,
+        pol=pol,
+        use_pixel_weights=use_pixel_weights,
+    )
+    new = hp.map2alm(
+        kappa,
+        lmax=healpix_inputs.lmax,
+        pol=pol,
+        use_pixel_weights=use_pixel_weights,
     )
+    xpx.testing.assert_equal(xp.asarray(old), new)
 
 
 @pytest.mark.parametrize(
@@ -235,42 +232,37 @@ def test_map2alm_individual(
     ],
 )
 def test_map2alm_sequence(
-    compare: type[Compare],
     healpix_inputs: type[HealpixInputs],
     pol: bool,  # noqa: FBT001
     urng: UnifiedGenerator,
     use_pixel_weights: bool,  # noqa: FBT001
+    xp: ModuleType,
 ) -> None:
     """Compare ``glass.healpix.map2alm`` against ``healpy.map2alm``."""
     kappa1 = healpix_inputs.kappa(rng=urng)
     kappa2 = healpix_inputs.kappa(rng=urng)
     kappa3 = healpix_inputs.kappa(rng=urng)
-    compare.assert_array_equal(
-        healpy.map2alm(
-            [np.asarray(kappa1), np.asarray(kappa2), np.asarray(kappa3)],
-            lmax=healpix_inputs.lmax,
-            pol=pol,
-            use_pixel_weights=use_pixel_weights,
-        ),
-        hp.map2alm(
-            [kappa1, kappa2, kappa3],
-            lmax=healpix_inputs.lmax,
-            pol=pol,
-            use_pixel_weights=use_pixel_weights,
-        ),
+    old = healpy.map2alm(
+        [np.asarray(kappa1), np.asarray(kappa2), np.asarray(kappa3)],
+        lmax=healpix_inputs.lmax,
+        pol=pol,
+        use_pixel_weights=use_pixel_weights,
     )
+    new = hp.map2alm(
+        [kappa1, kappa2, kappa3],
+        lmax=healpix_inputs.lmax,
+        pol=pol,
+        use_pixel_weights=use_pixel_weights,
+    )
+    xpx.testing.assert_equal(xp.asarray(old), new)
 
 
-def test_npix2nside(
-    healpix_inputs: type[HealpixInputs],
-) -> None:
+def test_npix2nside(healpix_inputs: type[HealpixInputs]) -> None:
     """Compare ``glass.healpix.npix2nside`` against ``healpix.npix2nside``."""
     assert healpix.npix2nside(healpix_inputs.npix) == hp.npix2nside(healpix_inputs.npix)
 
 
-def test_nside2npix(
-    healpix_inputs: type[HealpixInputs],
-) -> None:
+def test_nside2npix(healpix_inputs: type[HealpixInputs]) -> None:
     """Compare ``glass.healpix.nside2npix`` against ``healpix.nside2npix``."""
     assert healpix.nside2npix(healpix_inputs.nside) == hp.nside2npix(
         healpix_inputs.nside,
@@ -279,7 +271,6 @@ def test_nside2npix(
 
 @pytest.mark.parametrize("pol", [False, True])
 def test_pixwin(
-    compare: type[Compare],
     healpix_inputs: type[HealpixInputs],
     pol: bool,  # noqa: FBT001
     xp: ModuleType,
@@ -294,12 +285,11 @@ def test_pixwin(
 
     assert len(old) == len(new)
     for i in range(len(old)):
-        compare.assert_array_equal(old[i], new[i])
+        xpx.testing.assert_equal(xp.asarray(old[i], dtype=xp.float64), new[i])
 
 
 @pytest.mark.parametrize("thetas", [((20, 80)), ((30, 90))])
 def test_query_strip_float64(
-    compare: type[Compare],
     healpix_inputs: type[HealpixInputs],
     thetas: tuple[int, int],
     xp: ModuleType,
@@ -316,7 +306,7 @@ def test_query_strip_float64(
     old[healpy.query_strip(healpix_inputs.nside, *thetas)] = 0
     new = xp.ones(healpix_inputs.npix)
     new *= 1 - hp.query_strip(healpix_inputs.nside, thetas, dtype=xp.float64, xp=xp)
-    compare.assert_array_equal(old, new)
+    xpx.testing.assert_equal(xp.asarray(old), new)
 
 
 @pytest.mark.skipif(not HAVE_ARRAY_API_STRICT, reason="test requires array_api_strict")
@@ -340,7 +330,6 @@ def test_query_strip_none(
 
 @pytest.mark.parametrize("lonlat", [False, True])
 def test_randang(
-    compare: type[Compare],
     healpix_inputs: type[HealpixInputs],
     lonlat: bool,  # noqa: FBT001
     xp: ModuleType,
@@ -368,15 +357,15 @@ def test_randang(
     assert type(old) is type(new)
     assert len(old) == len(new)
     for i in range(len(old)):
-        compare.assert_array_equal(old[i], new[i])
+        xpx.testing.assert_equal(xp.asarray(old[i]), new[i])
 
 
 @pytest.mark.parametrize("coord", ["CE", "GC"])
 def test_rotate_map_pixel(
-    compare: type[Compare],
     coord: str,
     healpix_inputs: type[HealpixInputs],
     urng: UnifiedGenerator,
+    xp: ModuleType,
 ) -> None:
     """
     Compare ``glass.healpix.Rotator.rotate_map_pixel`` against
@@ -384,7 +373,6 @@ def test_rotate_map_pixel(
 
     """  # noqa: D205
     kappa = healpix_inputs.kappa(rng=urng)
-    compare.assert_array_equal(
-        healpy.Rotator(coord=coord).rotate_map_pixel(np.asarray(kappa)),
-        hp.Rotator(coord=coord).rotate_map_pixel(kappa),
-    )
+    old = healpy.Rotator(coord=coord).rotate_map_pixel(np.asarray(kappa))
+    new = hp.Rotator(coord=coord).rotate_map_pixel(kappa)
+    xpx.testing.assert_equal(xp.asarray(old), new)
diff --git a/tests/core/test_lensing.py b/tests/core/test_lensing.py
index 8f97ab51..17d114c4 100644
--- a/tests/core/test_lensing.py
+++ b/tests/core/test_lensing.py
@@ -5,6 +5,8 @@
 
 import pytest
 
+import array_api_extra as xpx
+
 import glass
 import glass.healpix as hp
 from glass._array_api_utils import xp_additions as uxpx
@@ -14,13 +16,9 @@
 
     from glass._types import FloatArray, UnifiedGenerator
     from glass.cosmology import Cosmology
-    from tests.fixtures.helper_classes import Compare
 
 
-def test_from_convergence(
-    compare: type[Compare],
-    urng: UnifiedGenerator,
-) -> None:
+def test_from_convergence(urng: UnifiedGenerator) -> None:
     """Add unit tests for :func:`glass.from_convergence`."""
     # l_max = 100  # noqa: ERA001
     n_side = 32
@@ -32,30 +30,30 @@ def test_from_convergence(
     # check with all False
 
     results = glass.from_convergence(kappa)
-    compare.assert_array_equal(results, ())
+    assert results == ()
 
     # check all combinations of potential, deflection, shear being True
 
     results = glass.from_convergence(kappa, potential=True)
-    compare.assert_array_equal(len(results), 1)
+    assert len(results) == 1
 
     results = glass.from_convergence(kappa, deflection=True)
-    compare.assert_array_equal(len(results), 1)
+    assert len(results) == 1
 
     results = glass.from_convergence(kappa, shear=True)
-    compare.assert_array_equal(len(results), 1)
+    assert len(results) == 1
 
     results = glass.from_convergence(kappa, potential=True, deflection=True)
-    compare.assert_array_equal(len(results), 2)
+    assert len(results) == 2
 
     results = glass.from_convergence(kappa, potential=True, shear=True)
-    compare.assert_array_equal(len(results), 2)
+    assert len(results) == 2
 
     results = glass.from_convergence(kappa, deflection=True, shear=True)
-    compare.assert_array_equal(len(results), 2)
+    assert len(results) == 2
 
     results = glass.from_convergence(kappa, potential=True, deflection=True, shear=True)
-    compare.assert_array_equal(len(results), 3)
+    assert len(results) == 3
 
 
 def test_shear_from_convergence() -> None:
@@ -64,7 +62,6 @@ def test_shear_from_convergence() -> None:
 
 
 def test_multi_plane_matrix(
-    compare: type[Compare],
     cosmo: Cosmology,
     shells: list[glass.RadialWindow],
     urng: UnifiedGenerator,
@@ -72,8 +69,8 @@ def test_multi_plane_matrix(
 ) -> None:
     mat = glass.multi_plane_matrix(shells, cosmo)
 
-    compare.assert_array_equal(mat, xp.tril(mat))
-    compare.assert_array_equal(xp.triu(mat, k=1), 0)
+    xpx.testing.assert_equal(mat, xp.tril(mat))
+    xpx.testing.assert_equal(xp.triu(mat, k=1), xp.asarray(0.0), check_shape=False)
 
     convergence = glass.MultiPlaneConvergence(cosmo)
 
@@ -86,11 +83,10 @@ def test_multi_plane_matrix(
         if convergence.kappa is not None:
             kappas.append(xp.asarray(convergence.kappa, copy=True))
 
-    compare.assert_allclose(mat @ deltas, kappas)
+    xpx.testing.assert_close(mat @ deltas, xp.stack(kappas))
 
 
 def test_multi_plane_weights(
-    compare: type[Compare],
     cosmo: Cosmology,
     urng: UnifiedGenerator,
     shells: list[glass.RadialWindow],
@@ -99,8 +95,8 @@ def test_multi_plane_weights(
     w_in = xp.eye(len(shells))
     w_out = glass.multi_plane_weights(w_in, shells, cosmo)
 
-    compare.assert_array_equal(w_out, xp.triu(w_out, k=1))
-    compare.assert_array_equal(xp.tril(w_out), 0)
+    xpx.testing.assert_equal(w_out, xp.triu(w_out, k=1))
+    xpx.testing.assert_equal(xp.tril(w_out), xp.asarray(0.0), check_shape=False)
 
     convergence = glass.MultiPlaneConvergence(cosmo)
 
@@ -118,12 +114,11 @@ def test_multi_plane_weights(
 
     wmat = glass.multi_plane_weights(weights, shells, cosmo)
 
-    compare.assert_allclose(uxpx.einsum("ij,ik", wmat, deltas), kappa)
+    xpx.testing.assert_close(uxpx.einsum("ij,ik", wmat, deltas), kappa)
 
 
 @pytest.mark.parametrize("usecomplex", [True, False])
 def test_deflect_nsew(
-    compare: type[Compare],
     usecomplex: bool,  # noqa: FBT001
     xp: ModuleType,
 ) -> None:
@@ -140,19 +135,19 @@ def alpha(
 
     # north
     lon, lat = glass.deflect(0.0, 0.0, alpha(r, 0, usecomplex=usecomplex), xp=xp)
-    compare.assert_allclose([lon, lat], [0.0, d], atol=1e-15)
+    xpx.testing.assert_close(xp.stack([lon, lat]), xp.asarray([0.0, d]), atol=1e-15)
 
     # south
     lon, lat = glass.deflect(0.0, 0.0, alpha(-r, 0, usecomplex=usecomplex), xp=xp)
-    compare.assert_allclose([lon, lat], [0.0, -d], atol=1e-15)
+    xpx.testing.assert_close(xp.stack([lon, lat]), xp.asarray([0.0, -d]), atol=1e-15)
 
     # east
     lon, lat = glass.deflect(0.0, 0.0, alpha(0, r, usecomplex=usecomplex), xp=xp)
-    compare.assert_allclose([lon, lat], [-d, 0.0], atol=1e-15)
+    xpx.testing.assert_close(xp.stack([lon, lat]), xp.asarray([-d, 0.0]), atol=1e-15)
 
     # west
     lon, lat = glass.deflect(0.0, 0.0, alpha(0, -r, usecomplex=usecomplex), xp=xp)
-    compare.assert_allclose([lon, lat], [d, 0.0], atol=1e-15)
+    xpx.testing.assert_close(xp.stack([lon, lat]), xp.asarray([d, 0.0]), atol=1e-15)
 
     # At least one input is an array
     lon, lat = glass.deflect(
@@ -160,15 +155,15 @@ def alpha(
         xp.asarray(0.0),
         alpha(0, -r, usecomplex=usecomplex),
     )
-    compare.assert_allclose([lon, lat], [d, 0.0], atol=1e-15)
+    xpx.testing.assert_close(xp.stack([lon, lat]), xp.asarray([d, 0.0]), atol=1e-15)
 
     lon, lat = glass.deflect(
         xp.asarray([0.0, 0.0]),
         xp.asarray([0.0, 0.0]),
         alpha(0, -r, usecomplex=usecomplex),
     )
-    compare.assert_allclose(lon, xp.asarray([d, d]), atol=1e-15)
-    compare.assert_allclose(lat, 0.0, atol=1e-15)
+    xpx.testing.assert_close(lon, xp.asarray([d, d]), atol=1e-15)
+    xpx.testing.assert_close(lat, xp.asarray(0.0), atol=1e-15, check_shape=False)
 
     # No inputs are arrays and xp not provided
     with pytest.raises(
@@ -179,7 +174,6 @@ def alpha(
 
 
 def test_deflect_many(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -197,4 +191,4 @@ def test_deflect_many(
 
     dotp = x * x_ + y * y_ + z * z_
 
-    compare.assert_allclose(dotp, xp.cos(abs_alpha))
+    xpx.testing.assert_close(dotp, xp.cos(abs_alpha))
diff --git a/tests/core/test_observations.py b/tests/core/test_observations.py
index 3b9b802f..cac05ccf 100644
--- a/tests/core/test_observations.py
+++ b/tests/core/test_observations.py
@@ -6,6 +6,8 @@
 import numpy as np
 import pytest
 
+import array_api_extra as xpx
+
 import glass
 import glass.healpix as hp
 
@@ -13,31 +15,27 @@
     from types import ModuleType
 
     from glass._types import UnifiedGenerator
-    from tests.fixtures.helper_classes import Compare
 
 
-def test_vmap_galactic_ecliptic(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_vmap_galactic_ecliptic(xp: ModuleType) -> None:
     """Add unit tests for :func:`glass.vmap_galactic_ecliptic`."""
     n_side = 4
 
     # check shape
 
     vmap = glass.vmap_galactic_ecliptic(n_side, xp=xp)
-    compare.assert_array_equal(vmap.shape[0], hp.nside2npix(n_side))
+    assert vmap.shape[0] == hp.nside2npix(n_side)
 
     # Test without xp
 
     vmap = glass.vmap_galactic_ecliptic(n_side)
     assert vmap.__array_namespace__().__name__ == "numpy"
-    compare.assert_array_equal(vmap.shape[0], hp.nside2npix(n_side))
+    assert vmap.shape[0] == hp.nside2npix(n_side)
 
     # no rotation
 
     vmap = glass.vmap_galactic_ecliptic(n_side, galactic=(0, 0), ecliptic=(0, 0), xp=xp)
-    compare.assert_array_equal(vmap, xp.zeros_like(vmap))
+    xpx.testing.assert_equal(vmap, xp.zeros_like(vmap))
 
     # check errors raised
 
@@ -55,7 +53,6 @@ def test_vmap_galactic_ecliptic(
 
 
 def test_gaussian_nz(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -67,13 +64,17 @@ def test_gaussian_nz(
     # check passing in the norm
 
     nz = glass.gaussian_nz(z, mean, sigma, norm=0)
-    compare.assert_array_equal(nz, xp.zeros_like(nz))
+    xpx.testing.assert_equal(nz, xp.zeros_like(nz))
 
     # check the value of each entry is close to the norm
 
     norm = 1
     nz = glass.gaussian_nz(z, mean, sigma, norm=norm)
-    compare.assert_allclose(xp.sum(nz) / nz.shape[0], norm, rtol=1e-2)
+    xpx.testing.assert_close(
+        xp.sum(nz) / nz.shape[0],
+        xp.asarray(norm, dtype=xp.float64),
+        rtol=1e-2,
+    )
 
     # check multidimensionality size
 
@@ -86,10 +87,7 @@ def test_gaussian_nz(
     assert nz.shape == (z.shape[0], z.shape[0])
 
 
-def test_smail_nz(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_smail_nz(xp: ModuleType) -> None:
     """Add unit tests for :func:`glass.smail_nz`."""
     alpha = 1
     beta = 1
@@ -99,10 +97,10 @@ def test_smail_nz(
     # check passing in the norm
 
     pz = glass.smail_nz(z, mode, alpha, beta, norm=0)
-    compare.assert_array_equal(pz, xp.zeros_like(pz))
+    xpx.testing.assert_equal(pz, xp.zeros_like(pz))
 
 
-def test_fixed_zbins_default_xp(compare: type[Compare]) -> None:
+def test_fixed_zbins_default_xp() -> None:
     """Add unit tests for :func:`glass.fixed_zbins` with default xp."""
     zmin = 0.0
     zmax = 1.0
@@ -118,14 +116,14 @@ def test_fixed_zbins_default_xp(compare: type[Compare]) -> None:
     )
     zbins = glass.fixed_zbins(zmin, zmax, nbins=nbins)
     assert len(zbins) == nbins
-    compare.assert_allclose(zbins, expected_zbins, rtol=1e-15)
+    xpx.testing.assert_close(np.asarray(zbins), expected_zbins, rtol=1e-15)
 
     # check dz input
 
     dz = 0.2
     zbins = glass.fixed_zbins(zmin, zmax, dz=dz)
     assert len(zbins) == math.ceil((zmax - zmin) / dz)
-    compare.assert_allclose(zbins, expected_zbins, rtol=1e-15)
+    xpx.testing.assert_close(np.asarray(zbins), expected_zbins, rtol=1e-15)
 
     # check dz for spacing which results in a max value above zmax
 
@@ -140,10 +138,7 @@ def test_fixed_zbins_default_xp(compare: type[Compare]) -> None:
         glass.fixed_zbins(zmin, zmax, nbins=nbins, dz=dz)
 
 
-def test_fixed_zbins_xp_provided(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_fixed_zbins_xp_provided(xp: ModuleType) -> None:
     """Add unit tests for :func:`glass.fixed_zbins`."""
     zmin = 0.0
     zmax = 1.0
@@ -159,14 +154,14 @@ def test_fixed_zbins_xp_provided(
     )
     zbins = glass.fixed_zbins(zmin, zmax, nbins=nbins, xp=xp)
     assert len(zbins) == nbins
-    compare.assert_allclose(zbins, expected_zbins, rtol=1e-15)
+    xpx.testing.assert_close(xp.asarray(zbins), expected_zbins, rtol=1e-15)
 
     # check dz input
 
     dz = 0.2
     zbins = glass.fixed_zbins(zmin, zmax, dz=dz, xp=xp)
     assert len(zbins) == math.ceil((zmax - zmin) / dz)
-    compare.assert_allclose(zbins, expected_zbins, rtol=1e-15)
+    xpx.testing.assert_close(xp.asarray(zbins), expected_zbins, rtol=1e-15)
 
     # check dz for spacing which results in a max value above zmax
 
@@ -179,10 +174,7 @@ def test_fixed_zbins_xp_provided(
         glass.fixed_zbins(zmin, zmax, nbins=nbins, dz=dz, xp=xp)
 
 
-def test_equal_dens_zbins(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_equal_dens_zbins(xp: ModuleType) -> None:
     """Add unit tests for :func:`glass.equal_dens_zbins`."""
     z = xp.linspace(0, 1, 11)
     nbins = 5
@@ -196,17 +188,14 @@ def test_equal_dens_zbins(
         ],
     )
     zbins = glass.equal_dens_zbins(z, xp.ones_like(z), nbins)
-    compare.assert_allclose(zbins, expected_zbins, rtol=1e-15)
+    xpx.testing.assert_close(xp.asarray(zbins), expected_zbins, rtol=1e-15)
 
     # check output shape
 
     assert len(zbins) == nbins
 
 
-def test_tomo_nz_gausserr(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_tomo_nz_gausserr(xp: ModuleType) -> None:
     """Add unit tests for :func:`glass.tomo_nz_gausserr`."""
     sigma_0 = 0.1
     z = xp.linspace(0, 1, 11)
@@ -215,7 +204,7 @@ def test_tomo_nz_gausserr(
     # check zeros returned
 
     binned_nz = glass.tomo_nz_gausserr(z, xp.zeros_like(z), sigma_0, zbins)
-    compare.assert_array_equal(binned_nz, xp.zeros_like(binned_nz))
+    xpx.testing.assert_equal(binned_nz, xp.zeros_like(binned_nz))
 
     # check the shape of the output
 
diff --git a/tests/core/test_points.py b/tests/core/test_points.py
index 6bb9af64..424feec4 100644
--- a/tests/core/test_points.py
+++ b/tests/core/test_points.py
@@ -9,6 +9,7 @@
 
 import glass
 import glass.healpix as hp
+import glass.points
 from glass._array_api_utils import xp_additions as uxpx
 
 if TYPE_CHECKING:
@@ -19,11 +20,10 @@
     from pytest_mock import MockerFixture
 
     from glass._types import UnifiedGenerator
-    from tests.fixtures.helper_classes import Compare, DataTransformer
+    from tests.fixtures.helper_classes import DataTransformer
 
 
 def test_effective_bias(
-    compare: type[Compare],
     mocker: MockerFixture,
     xp: ModuleType,
 ) -> None:
@@ -34,21 +34,20 @@ def test_effective_bias(
 
     z = xp.linspace(0, 1, 10)
     bz = xp.zeros((10,))
-    compare.assert_allclose(glass.effective_bias(z, bz, w), 0.0)
+    xpx.testing.assert_close(glass.effective_bias(z, bz, w), xp.asarray(0.0))
 
     z = xp.zeros((10,))
     bz = xp.full_like(z, 0.5)
 
-    compare.assert_allclose(glass.effective_bias(z, bz, w), 0.0)
+    xpx.testing.assert_close(glass.effective_bias(z, bz, w), xp.asarray(0.0))
 
     z = xp.linspace(0, 1, 10)
     bz = xp.full_like(z, 0.5)
 
-    compare.assert_allclose(glass.effective_bias(z, bz, w), 0.25)
+    xpx.testing.assert_close(glass.effective_bias(z, bz, w), xp.asarray(0.25))
 
 
 def test_linear_bias(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -57,25 +56,24 @@ def test_linear_bias(
     delta = xp.zeros((2, 2))
     b = 2.0
 
-    compare.assert_allclose(glass.linear_bias(delta, b), xp.zeros((2, 2)))
+    xpx.testing.assert_close(glass.linear_bias(delta, b), xp.zeros((2, 2)))
 
     # test with 0 b
 
     delta = urng.normal(5, 1, size=(2, 2))
     b = 0.0
 
-    compare.assert_allclose(glass.linear_bias(delta, b), xp.zeros((2, 2)))
+    xpx.testing.assert_close(glass.linear_bias(delta, b), xp.zeros((2, 2)))
 
     # compare with original implementation
 
     delta = urng.normal(5, 1, size=(2, 2))
     b = 2.0
 
-    compare.assert_allclose(glass.linear_bias(delta, b), b * delta)
+    xpx.testing.assert_close(glass.linear_bias(delta, b), b * delta)
 
 
 def test_loglinear_bias(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -84,30 +82,27 @@ def test_loglinear_bias(
     delta = xp.zeros((2, 2))
     b = 2.0
 
-    compare.assert_allclose(glass.loglinear_bias(delta, b), xp.zeros((2, 2)))
+    xpx.testing.assert_close(glass.loglinear_bias(delta, b), xp.zeros((2, 2)))
 
     # test with 0 b
 
     delta = urng.normal(5, 1, size=(2, 2))
     b = 0.0
 
-    compare.assert_allclose(glass.loglinear_bias(delta, b), xp.zeros((2, 2)))
+    xpx.testing.assert_close(glass.loglinear_bias(delta, b), xp.zeros((2, 2)))
 
     # compare with numpy implementation
 
     delta = urng.normal(5, 1, size=(2, 2))
     b = 2.0
 
-    compare.assert_allclose(
+    xpx.testing.assert_close(
         glass.loglinear_bias(delta, b),
         xp.expm1(b * xp.log1p(delta)),
     )
 
 
-def test_broadcast_inputs(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_broadcast_inputs(xp: ModuleType) -> None:
     bias_in = 0.8
     delta_in = xp.zeros((3, 1, 12))
     ngal_in = xp.asarray([1e-3, 2e-3])
@@ -123,11 +118,11 @@ def test_broadcast_inputs(
     assert dims == (3, 2)
     assert bias.shape == dims  # ty: ignore[unresolved-attribute]
     assert xp.all(bias == bias_in)
-    compare.assert_array_equal(delta, xp.zeros_like(delta))
+    xpx.testing.assert_equal(delta, xp.zeros_like(delta))
     assert ngal.shape == dims  # ty: ignore[unresolved-attribute]
-    compare.assert_array_equal(ngal[0, :], ngal_in)  # ty: ignore[not-subscriptable]
+    xpx.testing.assert_equal(ngal[0, :], ngal_in)  # ty: ignore[not-subscriptable]
     assert vis.shape == delta.shape  # ty: ignore[unresolved-attribute]
-    compare.assert_array_equal(vis[0, 0, :], vis_in)  # ty: ignore[not-subscriptable]
+    xpx.testing.assert_equal(vis[0, 0, :], vis_in)  # ty: ignore[not-subscriptable]
 
 
 @pytest.mark.parametrize(
@@ -139,7 +134,6 @@ def test_broadcast_inputs(
 )
 def test_compute_density_contrast(
     bias_model: Callable[..., Any],
-    compare: type[Compare],
     xp: ModuleType,
 ) -> None:
     bias = 0.8 * xp.ones((3, 2))
@@ -155,7 +149,7 @@ def test_compute_density_contrast(
     )
 
     assert n.shape[0] == delta.shape[-1]
-    compare.assert_array_equal(n, xp.zeros_like(n))
+    xpx.testing.assert_equal(n, xp.zeros_like(n))
 
 
 @pytest.mark.parametrize("remove_monopole", [False, True])
@@ -179,10 +173,7 @@ def test_compute_expected_count(
     assert xp.all(n == n[0])
 
 
-def test_apply_visibility(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_apply_visibility(xp: ModuleType) -> None:
     k = (1, 1)
     n_in = 24751.77674965 * xp.ones(12)
     vis = xp.tile(xp.repeat(xp.asarray([0.0, 1.0]), 6), (3, 2, 1))
@@ -193,24 +184,20 @@ def test_apply_visibility(
         vis,
     )
 
-    compare.assert_array_equal(n[:6], xp.zeros_like(n[:6]))
-    compare.assert_array_equal(n[6:], n_in[6:])
+    xpx.testing.assert_equal(n[:6], xp.zeros_like(n[:6]))
+    xpx.testing.assert_equal(n[6:], n_in[6:])
 
 
-def test_sample_number_galaxies(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_sample_number_galaxies(xp: ModuleType) -> None:
     n_in = xp.repeat(xp.asarray([0.0, 24751.77674965]), 6)
 
     n = glass.points._sample_number_galaxies(n_in)
 
-    compare.assert_array_equal(n[:6], xp.zeros_like(n[:6]))
-    compare.assert_allclose(n[6:], n_in[6:], atol=250)
+    xpx.testing.assert_equal(n[:6], xp.zeros_like(n[:6], dtype=xp.int64))
+    xpx.testing.assert_close(xp.asarray(n[6:], dtype=xp.float64), n_in[6:], atol=250)
 
 
 def test_sample_number_galaxies_rng(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -218,8 +205,8 @@ def test_sample_number_galaxies_rng(
 
     n = glass.points._sample_number_galaxies(n_in, rng=urng)
 
-    compare.assert_array_equal(n[:6], xp.zeros_like(n[:6]))
-    compare.assert_allclose(n[6:], n_in[6:], atol=250)
+    xpx.testing.assert_equal(n[:6], xp.zeros_like(n[:6]))
+    xpx.testing.assert_close(xp.asarray(n[6:], dtype=xp.float64), n_in[6:], atol=250)
 
 
 def test_sample_galaxies_per_pixel(
@@ -247,7 +234,6 @@ def test_sample_galaxies_per_pixel(
 
 
 def test_positions_from_delta(  # noqa: PLR0915
-    compare: type[Compare],
     data_transformer: type[DataTransformer],
     urng: UnifiedGenerator,
     xp: ModuleType,
@@ -319,8 +305,8 @@ def test_positions_from_delta(  # noqa: PLR0915
     )
 
     assert int(count) == count
-    compare.assert_allclose(lon, [])
-    compare.assert_allclose(lat, [])
+    xpx.testing.assert_equal(lon, xp.asarray([]))
+    xpx.testing.assert_equal(lat, xp.asarray([]))
 
     # case: large delta
 
@@ -453,7 +439,6 @@ def test_uniform_positions(
 
 
 def test_position_weights(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -479,13 +464,10 @@ def test_position_weights(
                     )
                 expected = bias * expected
 
-            compare.assert_allclose(weights, expected)
+            xpx.testing.assert_close(weights, expected)
 
 
-def test_displace_arg_complex(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_displace_arg_complex(xp: ModuleType) -> None:
     """Test displace function with complex-valued displacement."""
     d = 5.0  # deg
     r = d / 180 * math.pi
@@ -496,25 +478,22 @@ def test_displace_arg_complex(
 
     # north
     lon, lat = glass.displace(lon0, lat0, xp.asarray(r + 0j))
-    compare.assert_allclose([lon, lat], [0.0, d])
+    xpx.testing.assert_close(xp.stack([lon, lat]), xp.asarray([0.0, d]))
 
     # south
     lon, lat = glass.displace(lon0, lat0, xp.asarray(-r + 0j))
-    compare.assert_allclose([lon, lat], [0.0, -d], atol=1e-15)
+    xpx.testing.assert_close(xp.stack([lon, lat]), xp.asarray([0.0, -d]), atol=1e-15)
 
     # east
     lon, lat = glass.displace(lon0, lat0, xp.asarray(1j * r))
-    compare.assert_allclose([lon, lat], [d, 0.0], atol=1e-15)
+    xpx.testing.assert_close(xp.stack([lon, lat]), xp.asarray([d, 0.0]), atol=1e-15)
 
     # west
     lon, lat = glass.displace(lon0, lat0, xp.asarray(-1j * r))
-    compare.assert_allclose([lon, lat], [-d, 0.0], atol=1e-15)
+    xpx.testing.assert_close(xp.stack([lon, lat]), xp.asarray([-d, 0.0]), atol=1e-15)
 
 
-def test_displace_arg_real(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_displace_arg_real(xp: ModuleType) -> None:
     """Test displace function with real-valued argument."""
     d = 5.0  # deg
     r = d / 180 * math.pi
@@ -525,23 +504,22 @@ def test_displace_arg_real(
 
     # north
     lon, lat = glass.displace(lon0, lat0, xp.asarray([r, 0]))
-    compare.assert_allclose([lon, lat], [0.0, d])
+    xpx.testing.assert_close(xp.stack([lon, lat]), xp.asarray([0.0, d]))
 
     # south
     lon, lat = glass.displace(lon0, lat0, xp.asarray([-r, 0]))
-    compare.assert_allclose([lon, lat], [0.0, -d], atol=1e-15)
+    xpx.testing.assert_close(xp.stack([lon, lat]), xp.asarray([0.0, -d]), atol=1e-15)
 
     # east
     lon, lat = glass.displace(lon0, lat0, xp.asarray([0, r]))
-    compare.assert_allclose([lon, lat], [d, 0.0], atol=1e-15)
+    xpx.testing.assert_close(xp.stack([lon, lat]), xp.asarray([d, 0.0]), atol=1e-15)
 
     # west
     lon, lat = glass.displace(lon0, lat0, xp.asarray([0, -r]))
-    compare.assert_allclose([lon, lat], [-d, 0.0], atol=1e-15)
+    xpx.testing.assert_close(xp.stack([lon, lat]), xp.asarray([-d, 0.0]), atol=1e-15)
 
 
 def test_displace_abs(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -561,11 +539,10 @@ def test_displace_abs(
 
     cos_a = xp.cos(th) * xp.cos(th_) + xp.cos(delta) * xp.sin(th) * xp.sin(th_)
 
-    compare.assert_allclose(cos_a, xp.cos(abs_alpha))
+    xpx.testing.assert_close(cos_a, xp.cos(abs_alpha))
 
 
 def test_displacement(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -598,7 +575,7 @@ def test_displacement(
     # test each displacement individually
     for from_lon, from_lat, to_lon, to_lat, alpha in data:
         alpha_ = glass.displacement(from_lon, from_lat, to_lon, to_lat)
-        compare.assert_allclose(alpha_, alpha)
+        xpx.testing.assert_close(alpha_, alpha)
 
     # test on an array
     alpha = glass.displacement(
@@ -611,7 +588,6 @@ def test_displacement(
 
 
 def test_displacement_zerodist(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -619,14 +595,13 @@ def test_displacement_zerodist(
     lon = urng.uniform(-180.0, 180.0, size=100)
     lat = urng.uniform(-90.0, 90.0, size=100)
 
-    compare.assert_allclose(
+    xpx.testing.assert_close(
         glass.displacement(lon, lat, lon, lat),
-        xp.zeros(100),
+        xp.zeros(100, dtype=xp.complex128),
     )
 
 
 def test_displacement_consistent(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -650,11 +625,10 @@ def test_displacement_consistent(
     # measure displacement
     alpha_out = glass.displacement(from_lon, from_lat, to_lon, to_lat)
 
-    compare.assert_allclose(alpha_out, alpha_in, atol=0.0, rtol=1e-10)
+    xpx.testing.assert_close(alpha_out, alpha_in, atol=0.0, rtol=1e-10)
 
 
 def test_displacement_random(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -699,11 +673,11 @@ def test_displacement_random(
         ],
         axis=1,
     )
-    compare.assert_allclose(rot @ xp.asarray([0.0, 0.0, 1.0]), u)
+    xpx.testing.assert_close(rot @ xp.asarray([0.0, 0.0, 1.0]), u)
 
     # meta-check that recovering theta and phi from vector works
-    compare.assert_allclose(xp.atan2(xp.hypot(u[:, 0], u[:, 1]), u[:, 2]), theta)
-    compare.assert_allclose(xp.atan2(u[:, 1], u[:, 0]), phi)
+    xpx.testing.assert_close(xp.atan2(xp.hypot(u[:, 0], u[:, 1]), u[:, 2]), theta)
+    xpx.testing.assert_close(xp.atan2(u[:, 1], u[:, 0]), phi)
 
     # build the displaced points near (0, 0, 1) and rotate near theta and phi
     v = xp.stack(
@@ -730,4 +704,4 @@ def test_displacement_random(
     # compute displacement and compare to input
     alpha_in = r * xp.exp(1j * x)
     alpha_out = glass.displacement(from_lon, from_lat, to_lon, to_lat)
-    compare.assert_allclose(alpha_out, alpha_in, atol=0.0, rtol=1e-10)
+    xpx.testing.assert_close(alpha_out, alpha_in, atol=0.0, rtol=1e-10)
diff --git a/tests/core/test_shapes.py b/tests/core/test_shapes.py
index a4d09522..554cba2e 100644
--- a/tests/core/test_shapes.py
+++ b/tests/core/test_shapes.py
@@ -4,13 +4,14 @@
 
 import pytest
 
+import array_api_extra as xpx
+
 import glass
 
 if TYPE_CHECKING:
     from types import ModuleType
 
     from glass._types import UnifiedGenerator
-    from tests.fixtures.helper_classes import Compare
 
 
 def test_triaxial_axis_ratio(
@@ -127,7 +128,6 @@ def test_ellipticity_ryden04(
 
 @pytest.mark.flaky(rerun=5, only_rerun=["AssertionError"])
 def test_ellipticity_gaussian(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -151,25 +151,44 @@ def test_ellipticity_gaussian(
 
     assert eps.shape == (n,)
 
-    compare.assert_array_less(xp.abs(eps), 1)
+    xpx.testing.assert_less(xp.abs(eps), xp.asarray(1.0), check_shape=False)
 
-    compare.assert_allclose(xp.std(xp.real(eps)), 0.256, atol=1e-3, rtol=0)
-    compare.assert_allclose(xp.std(xp.imag(eps)), 0.256, atol=1e-3, rtol=0)
+    xpx.testing.assert_close(xp.std(xp.real(eps)), xp.asarray(0.256), atol=1e-3, rtol=0)
+    xpx.testing.assert_close(xp.std(xp.imag(eps)), xp.asarray(0.256), atol=1e-3, rtol=0)
 
     eps = glass.ellipticity_gaussian(xp.asarray([n, n]), xp.asarray([0.128, 0.256]))
 
     assert eps.shape == (2 * n,)
 
-    compare.assert_array_less(xp.abs(eps), 1)
+    xpx.testing.assert_less(xp.abs(eps), xp.asarray(1.0), check_shape=False)
 
-    compare.assert_allclose(xp.std(xp.real(eps)[:n]), 0.128, atol=1e-3, rtol=0)
-    compare.assert_allclose(xp.std(xp.imag(eps)[:n]), 0.128, atol=1e-3, rtol=0)
-    compare.assert_allclose(xp.std(xp.real(eps)[n:]), 0.256, atol=1e-3, rtol=0)
-    compare.assert_allclose(xp.std(xp.imag(eps)[n:]), 0.256, atol=1e-3, rtol=0)
+    xpx.testing.assert_close(
+        xp.std(xp.real(eps)[:n]),
+        xp.asarray(0.128),
+        atol=1e-3,
+        rtol=0,
+    )
+    xpx.testing.assert_close(
+        xp.std(xp.imag(eps)[:n]),
+        xp.asarray(0.128),
+        atol=1e-3,
+        rtol=0,
+    )
+    xpx.testing.assert_close(
+        xp.std(xp.real(eps)[n:]),
+        xp.asarray(0.256),
+        atol=1e-3,
+        rtol=0,
+    )
+    xpx.testing.assert_close(
+        xp.std(xp.imag(eps)[n:]),
+        xp.asarray(0.256),
+        atol=1e-3,
+        rtol=0,
+    )
 
 
 def test_ellipticity_intnorm(
-    compare: type[Compare],
     urng: UnifiedGenerator,
     xp: ModuleType,
 ) -> None:
@@ -193,21 +212,41 @@ def test_ellipticity_intnorm(
 
     assert eps.shape == (n,)
 
-    compare.assert_array_less(xp.abs(eps), 1)
+    xpx.testing.assert_less(xp.abs(eps), xp.asarray(1.0), check_shape=False)
 
-    compare.assert_allclose(xp.std(xp.real(eps)), 0.256, atol=1e-3, rtol=0)
-    compare.assert_allclose(xp.std(xp.imag(eps)), 0.256, atol=1e-3, rtol=0)
+    xpx.testing.assert_close(xp.std(xp.real(eps)), xp.asarray(0.256), atol=1e-3, rtol=0)
+    xpx.testing.assert_close(xp.std(xp.imag(eps)), xp.asarray(0.256), atol=1e-3, rtol=0)
 
     eps = glass.ellipticity_intnorm(xp.asarray([n, n]), xp.asarray([0.128, 0.256]))
 
     assert eps.shape == (2 * n,)
 
-    compare.assert_array_less(xp.abs(eps), 1)
+    xpx.testing.assert_less(xp.abs(eps), xp.asarray(1.0), check_shape=False)
 
-    compare.assert_allclose(xp.std(xp.real(eps)[:n]), 0.128, atol=1e-3, rtol=0)
-    compare.assert_allclose(xp.std(xp.imag(eps)[:n]), 0.128, atol=1e-3, rtol=0)
-    compare.assert_allclose(xp.std(xp.real(eps)[n:]), 0.256, atol=1e-3, rtol=0)
-    compare.assert_allclose(xp.std(xp.imag(eps)[n:]), 0.256, atol=1e-3, rtol=0)
+    xpx.testing.assert_close(
+        xp.std(xp.real(eps)[:n]),
+        xp.asarray(0.128),
+        atol=1e-3,
+        rtol=0,
+    )
+    xpx.testing.assert_close(
+        xp.std(xp.imag(eps)[:n]),
+        xp.asarray(0.128),
+        atol=1e-3,
+        rtol=0,
+    )
+    xpx.testing.assert_close(
+        xp.std(xp.real(eps)[n:]),
+        xp.asarray(0.256),
+        atol=1e-3,
+        rtol=0,
+    )
+    xpx.testing.assert_close(
+        xp.std(xp.imag(eps)[n:]),
+        xp.asarray(0.256),
+        atol=1e-3,
+        rtol=0,
+    )
 
     with pytest.raises(ValueError, match="sigma must be between"):
         glass.ellipticity_intnorm(1, 0.71, xp=xp)
diff --git a/tests/core/test_shells.py b/tests/core/test_shells.py
index 072283ff..29895915 100644
--- a/tests/core/test_shells.py
+++ b/tests/core/test_shells.py
@@ -7,6 +7,8 @@
 import numpy as np
 import pytest
 
+import array_api_extra as xpx
+
 import glass
 from glass._array_api_utils import xp_additions as uxpx
 
@@ -14,11 +16,9 @@
     from types import ModuleType
 
     from glass.cosmology import Cosmology
-    from tests.fixtures.helper_classes import Compare
 
 
 def test_distance_weight(
-    compare: type[Compare],
     cosmo: Cosmology,
     xp: ModuleType,
 ) -> None:
@@ -28,7 +28,7 @@ def test_distance_weight(
     # check shape
 
     w = glass.DistanceWeight(cosmo)(z)
-    compare.assert_array_equal(w.shape, z.shape)
+    assert w.shape == z.shape
 
     # check first value is 1
 
@@ -36,11 +36,10 @@ def test_distance_weight(
 
     # check values are decreasing
 
-    compare.assert_array_less(w[1:], w[:-1])
+    xpx.testing.assert_less(w[1:], w[:-1])
 
 
 def test_volume_weight(
-    compare: type[Compare],
     cosmo: Cosmology,
     xp: ModuleType,
 ) -> None:
@@ -50,7 +49,7 @@ def test_volume_weight(
     # check shape
 
     w = glass.VolumeWeight(cosmo)(z)
-    compare.assert_array_equal(w.shape, z.shape)
+    assert w.shape == z.shape
 
     # check first value is 0
 
@@ -58,20 +57,17 @@ def test_volume_weight(
 
     # check values are increasing
 
-    compare.assert_array_less(w[:-1], w[1:])
+    xpx.testing.assert_less(w[:-1], w[1:])
 
 
-def test_density_weight(
-    compare: type[Compare],
-    cosmo: Cosmology,
-) -> None:
+def test_density_weight(cosmo: Cosmology) -> None:
     """Add unit tests for :class:`glass.DensityWeight`."""
     z = np.linspace(0, 1, 6)
 
     # check shape
 
     w = glass.DensityWeight(cosmo)(z)
-    compare.assert_array_equal(w.shape, z.shape)
+    assert w.shape == z.shape
 
     # check first value is 0
 
@@ -79,7 +75,7 @@ def test_density_weight(
 
     # check values are increasing
 
-    compare.assert_array_less(w[:-1], w[1:])
+    xpx.testing.assert_less(w[:-1], w[1:])
 
 
 def test_tophat_windows(xp: ModuleType) -> None:
@@ -104,12 +100,9 @@ def test_tophat_windows(xp: ModuleType) -> None:
     assert all(xp.all(w.wa == 1) for w in ws)
 
 
-def test_linear_windows(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_linear_windows(xp: ModuleType) -> None:
     """Add unit tests for :func:`glass.linear_windows`."""
-    dz = 1e-2
+    dz = xp.asarray(1e-2)
     zgrid = xp.asarray(
         [
             0.0,
@@ -123,7 +116,7 @@ def test_linear_windows(
     # check spacing of redshift grid
 
     ws = glass.linear_windows(zgrid)
-    compare.assert_allclose(dz, xp.mean(xp.diff(ws[0].za)), atol=1e-2)
+    xpx.testing.assert_close(dz, xp.mean(xp.diff(ws[0].za)), atol=1e-2)
 
     # check number of windows
 
@@ -131,7 +124,7 @@ def test_linear_windows(
 
     # check values of zeff
 
-    compare.assert_allclose([w.zeff for w in ws], zgrid[1:-1])
+    xpx.testing.assert_close(xp.stack([w.zeff for w in ws]), zgrid[1:-1])
 
     # check weight function input
 
@@ -140,7 +133,7 @@ def test_linear_windows(
         weight=lambda _: 0,
     )
     for w in ws:
-        compare.assert_allclose(w.wa, xp.zeros_like(w.wa))
+        xpx.testing.assert_close(w.wa, xp.zeros_like(w.wa))
 
     # check error raised
 
@@ -156,12 +149,9 @@ def test_linear_windows(
         glass.linear_windows(xp.asarray([0.1, 0.2, 0.3]))
 
 
-def test_cubic_windows(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_cubic_windows(xp: ModuleType) -> None:
     """Add unit tests for :func:`glass.cubic_windows`."""
-    dz = 1e-2
+    dz = xp.asarray(1e-2)
     zgrid = xp.asarray(
         [
             0.0,
@@ -175,7 +165,7 @@ def test_cubic_windows(
     # check spacing of redshift grid
 
     ws = glass.cubic_windows(zgrid)
-    compare.assert_allclose(dz, xp.mean(xp.diff(ws[0].za)), atol=1e-2)
+    xpx.testing.assert_close(dz, xp.mean(xp.diff(ws[0].za)), atol=1e-2)
 
     # check number of windows
 
@@ -183,7 +173,7 @@ def test_cubic_windows(
 
     # check values of zeff
 
-    compare.assert_allclose([w.zeff for w in ws], zgrid[1:-1])
+    xpx.testing.assert_close(xp.stack([w.zeff for w in ws]), zgrid[1:-1])
 
     # check weight function input
 
@@ -192,7 +182,7 @@ def test_cubic_windows(
         weight=lambda _: 0,
     )
     for w in ws:
-        compare.assert_allclose(w.wa, xp.zeros_like(w.wa))
+        xpx.testing.assert_close(w.wa, xp.zeros_like(w.wa))
 
     # check error raised
 
@@ -243,7 +233,6 @@ def test_restrict(xp: ModuleType) -> None:
 
 @pytest.mark.parametrize("method", ["lstsq", "nnls", "restrict"])
 def test_partition(
-    compare: type[Compare],
     method: str,
     xp: ModuleType,
 ) -> None:
@@ -286,7 +275,7 @@ def test_partition(
 
     assert part.shape == (len(shells), 3, 2)
 
-    compare.assert_allclose(xp.sum(part, axis=0), uxpx.trapezoid(fz, z))
+    xpx.testing.assert_close(xp.sum(part, axis=0), uxpx.trapezoid(fz, z))
 
 
 def test_redshift_grid_default_xp() -> None:
@@ -349,10 +338,7 @@ def test_redshift_grid(xp: ModuleType) -> None:
         glass.redshift_grid(zmin, zmax, dz=dz, num=num, xp=xp)
 
 
-def test_distance_grid(
-    compare: type[Compare],
-    cosmo: Cosmology,
-) -> None:
+def test_distance_grid(cosmo: Cosmology) -> None:
     """Add unit tests for :func:`glass.distance_grid`."""
     zmin = 0.0
     zmax = 1.0
@@ -372,7 +358,7 @@ def test_distance_grid(
     # check decrease in distance
 
     x = glass.distance_grid(cosmo, zmin, zmax, dx=0.3)
-    compare.assert_array_less(x[1:], x[:-1])
+    xpx.testing.assert_less(x[1:], x[:-1])
 
     # check error raised
 
@@ -389,10 +375,7 @@ def test_distance_grid(
         glass.distance_grid(cosmo, zmin, zmax, dx=dx, num=num)
 
 
-def test_combine(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_combine(xp: ModuleType) -> None:
     """Add unit tests for :func:`glass.combine`."""
     z = xp.linspace(0.0, 5.0, 1_000)
     weights = xp.asarray(
@@ -428,10 +411,10 @@ def test_combine(
     assert result.shape == z.shape
 
     # Check sum of result
-    compare.assert_allclose(sum(result), 929.267284)
+    xpx.testing.assert_close(xp.sum(result), xp.asarray(929.267284))
 
     # Check integral w.r.t z has not changed
-    compare.assert_allclose(uxpx.trapezoid(result, z), 4.643139, rtol=1e-6)
+    xpx.testing.assert_close(uxpx.trapezoid(result, z), xp.asarray(4.643139), rtol=1e-6)
 
 
 def test_radial_window_immutable(xp: ModuleType) -> None:
@@ -461,10 +444,7 @@ def test_radial_window_immutable(xp: ModuleType) -> None:
         w.zeff = zeff  # ty: ignore[invalid-assignment]
 
 
-def test_radial_window_zeff_none(
-    compare: type[Compare],
-    xp: ModuleType,
-) -> None:
+def test_radial_window_zeff_none(xp: ModuleType) -> None:
     """Checks ``zeff`` is computed when not provided to :class:`RadialWindow`."""
     # check zeff is computed when not provided
 
@@ -473,7 +453,7 @@ def test_radial_window_zeff_none(
 
     w = glass.RadialWindow(za, wa)
 
-    compare.assert_allclose(w.zeff, 1.0)
+    xpx.testing.assert_close(w.zeff, xp.asarray(1.0))
 
     # check zeff is NaN when redshift array is empty
 
diff --git a/tests/core/test_user.py b/tests/core/test_user.py
index b22c7245..ecf6956a 100644
--- a/tests/core/test_user.py
+++ b/tests/core/test_user.py
@@ -11,7 +11,6 @@
 
 if TYPE_CHECKING:
     from glass._types import AngularPowerSpectra, FloatArray
-    from tests.fixtures.helper_classes import Compare
 
 # check if available for testing
 HAVE_FITSIO = importlib.util.find_spec("fitsio") is not None
@@ -24,7 +23,6 @@
 
 
 def test_read_write_cls(
-    compare: type[Compare],
     rng: np.random.Generator,
     tmp_path: pathlib.Path,
 ) -> None:
@@ -37,15 +35,15 @@ def test_read_write_cls(
         values = npz["values"]
         split = npz["split"]
 
-    compare.assert_array_equal(values, np.concat(cls))
-    compare.assert_array_equal(
+    np.testing.assert_array_equal(values, np.concat(cls))
+    np.testing.assert_array_equal(
         split,
         np.cumulative_sum([cl.shape[0] for cl in cls[:-1]]),
     )
-    compare.assert_array_equal(cls, np.split(values, split))
+    np.testing.assert_array_equal(cls, np.split(values, split))
 
     npz = glass.load_cls(tmp_path / cls_file)
-    compare.assert_array_equal(npz, cls)
+    np.testing.assert_array_equal(npz, cls)
 
 
 @pytest.mark.skipif(not HAVE_FITSIO, reason="test requires fitsio")
diff --git a/tests/fixtures/helper_classes.py b/tests/fixtures/helper_classes.py
index 05f08f28..bf83bbb3 100644
--- a/tests/fixtures/helper_classes.py
+++ b/tests/fixtures/helper_classes.py
@@ -4,7 +4,6 @@
 
 from typing import TYPE_CHECKING
 
-import numpy as np
 import pytest
 
 if TYPE_CHECKING:
@@ -13,7 +12,6 @@
     from typing import Any
 
     from glass._types import (
-        AnyArray,
         ComplexArray,
         FloatArray,
         IntArray,
@@ -21,54 +19,6 @@
     )
 
 
-class Compare:
-    """
-    Helper class for array comparisons in tests.
-
-    This class wraps numpy testing functions to provide a consistent interface
-    for comparing arrays in tests. Ultimately, it would be great if we can
-    make the array testing backend-agnostic.
-
-    """
-
-    @staticmethod
-    def assert_allclose(
-        actual: AnyArray,
-        desired: AnyArray,
-        *,
-        rtol: float = 1e-7,
-        atol: float = 0,
-    ) -> None:
-        """Check if two objects are not equal up to desired tolerance."""
-        np.testing.assert_allclose(actual, desired, rtol=rtol, atol=atol)
-
-    @staticmethod
-    def assert_array_almost_equal_nulp(
-        actual: AnyArray,
-        desired: AnyArray,
-        *,
-        nulp: int = 1,
-    ) -> None:
-        """Compare two arrays relatively to their spacing."""
-        np.testing.assert_array_almost_equal_nulp(actual, desired, nulp=nulp)
-
-    @staticmethod
-    def assert_array_equal(actual: AnyArray, desired: AnyArray) -> None:
-        """Check if two array objects are not equal."""
-        np.testing.assert_array_equal(actual, desired)
-
-    @staticmethod
-    def assert_array_less(actual: AnyArray, desired: AnyArray) -> None:
-        """Check if two array objects are not ordered by less than."""
-        np.testing.assert_array_less(actual, desired)
-
-
-@pytest.fixture(scope="session")
-def compare() -> type[Compare]:
-    """Fixture for array comparison utility."""
-    return Compare
-
-
 class DataTransformer:
     """Helper class for transforming various data structures into others."""