Repository: scipy/scipy
File: scipy/integrate/tests/test_quadrature.py (Python)

Code before:
        wz = np.ones_like(z) * (z[1] - z[0])
        wz[0] /= 2
        wz[-1] /= 2

        q = x[:, None, None] + y[None,:, None] + z[None, None,:]

        qx = (q * wx[:, None, None]).sum(axis=0)
        qy = (q * wy[None, :, None]).sum(axis=1)
        qz = (q * wz[None, None, :]).sum(axis=2)

        # n-d `x`
        r = trapezoid(q, x=x[:, None, None], axis=0)
        assert_allclose(r, qx)
        r = trapezoid(q, x=y[None,:, None], axis=1)
        assert_allclose(r, qy)
        r = trapezoid(q, x=z[None, None,:], axis=2)
        assert_allclose(r, qz)

        # 1-d `x`
        r = trapezoid(q, x=x, axis=0)
        assert_allclose(r, qx)
        r = trapezoid(q, x=y, axis=1)
        assert_allclose(r, qy)
        r = trapezoid(q, x=z, axis=2)
        assert_allclose(r, qz)

    def test_masked(self):
        # Testing that masked arrays behave as if the function is 0 where
        # masked
        x = np.arange(5)
        y = x * x
        mask = x == 2
        ym = np.ma.array(y, mask=mask)
        r = 13.0  # sum(0.5 * (0 + 1) * 1.0 + 0.5 * (9 + 16))
        assert_allclose(trapezoid(ym, x), r)

        xm = np.ma.array(x, mask=mask)
        assert_allclose(trapezoid(ym, xm), r)

        xm = np.ma.array(x, mask=mask)
        assert_allclose(trapezoid(y, xm), r)


@array_api_compatible
class TestSimpson:
    def test_simple(self, xp):
        x = xp.arange(-10, 10, .1)
        r = simpson(xp.exp(-.5 * x ** 2) / xp.sqrt(2 * xp.asarray(xp.pi)), dx=0.1)
        # check integral of normal equals 1
        xp_assert_close(r, xp.asarray(1.0))


Review comment (suggestion):
This seems to fix the `torch` issue locally (I guess you're going to replace the tests anyway) by allowing the expected value to overcome the default `dtype` used by `torch`.

```diff
--- a/scipy/integrate/tests/test_quadrature.py
+++ b/scipy/integrate/tests/test_quadrature.py
@@ -336,9 +336,9 @@ class TestSimpson:
         xp_assert_close(r, xp.asarray(1.0))
 
     def test_ndim(self, xp):
-        x = xp.linspace(0, 1, 3)
-        y = xp.linspace(0, 2, 8)
-        z = xp.linspace(0, 3, 13)
+        x = xp.linspace(0, 1, 3, dtype=xp.float64)
+        y = xp.linspace(0, 2, 8, dtype=xp.float64)
+        z = xp.linspace(0, 3, 13, dtype=xp.float64)
```

Code after:
    @pytest.mark.parametrize(
        "initial", [1, 0.5]
    )
    def test_initial_error(self, initial, xp):
        """If initial is not None or 0, a ValueError is raised."""
        y = xp.linspace(0, 10, num=10)
        with pytest.raises(ValueError, match="`initial`"):
            cumulative_trapezoid(y, initial=initial)

    def test_zero_len_y(self, xp):
        with pytest.raises(ValueError, match="At least one point is required"):
            cumulative_trapezoid(y=xp.asarray([]))


class CommonTrapezoidSimpsonTests:
    def test_simple(self, xp):
        x = xp.arange(-10, 10, .1)
        r = self.quadrature_func(xp.exp(-.5 * x ** 2) / xp.sqrt(2 * xp.asarray(xp.pi)),
                                 dx=0.1)
        # check integral of normal equals 1
        xp_assert_close(r, xp.asarray(1.0))

    def test_ndim(self, xp):
        x = xp.linspace(0, 1, 3)
        y = xp.linspace(0, 2, 8)
        z = xp.linspace(0, 3, 13)

        wx = xp.ones_like(x) * (x[1] - x[0])
        wx = xpx.at(wx)[0].divide(2)
        wx = xpx.at(wx)[-1].divide(2)
        wy = xp.ones_like(y) * (y[1] - y[0])
        wy = xpx.at(wy)[0].divide(2)
        wy = xpx.at(wy)[-1].divide(2)
        wz = xp.ones_like(z) * (z[1] - z[0])
        wz = xpx.at(wz)[0].divide(2)
        wz = xpx.at(wz)[-1].divide(2)

        q = x[:, None, None] + y[None,:, None] + z[None, None,:]

        qx = xp.sum(q * wx[:, None, None], axis=0)
        qy = xp.sum(q * wy[None, :, None], axis=1)
        qz = xp.sum(q * wz[None, None, :], axis=2)

        # n-d `x`
        r = self.quadrature_func(q, x=x[:, None, None], axis=0)
        xp_assert_close(r, qx)
        r = self.quadrature_func(q, x=y[None,:, None], axis=1)
        xp_assert_close(r, qy)
        r = self.quadrature_func(q, x=z[None, None,:], axis=2)
        xp_assert_close(r, qz)