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

Code before:
        ndim = 2
        mean = np.zeros(ndim)
        cov = np.eye(ndim)

        def func(x):
            # standard multivariate normal PDF in two dimensions
            return xp.exp(-0.5 * xp.sum(x*x, axis=0)) / (2 * xp.pi)

        rng = np.random.default_rng(2879434385674690281)
        qrng = stats.qmc.Sobol(ndim, seed=rng)
        a = np.zeros(ndim)
        b = np.ones(ndim) * signs
        res = qmc_quad(func, xp.asarray(a, dtype=dtype), xp.asarray(b, dtype=dtype),
                       n_points=n_points, n_estimates=n_estimates, qrng=qrng)
        ref = stats.multivariate_normal.cdf(b, mean, cov, lower_limit=a)
        atol = special.stdtrit(n_estimates-1, 0.995) * res.standard_error  # 99% CI
        xp_assert_close(res.integral, xp.asarray(ref, dtype=dtype), atol=atol)
        assert np.prod(signs)*res.integral > 0

        rng = np.random.default_rng(2879434385674690281)
        qrng = stats.qmc.Sobol(ndim, seed=rng)
        logres = qmc_quad(lambda *args: xp.log(func(*args)),
                          xp.asarray(a, dtype=dtype), xp.asarray(b, dtype=dtype),
                          n_points=n_points, n_estimates=n_estimates,
                          log=True, qrng=qrng)
        rtol = 1e-14 if res.integral.dtype == xp.float64 else 1e-6
        xp_assert_close(xp.real(xp.exp(logres.integral)), res.integral, rtol=rtol)
        assert xp.imag(logres.integral + 0j) == (xp.pi if np.prod(signs) < 0 else 0)
        xp_assert_close(xp.exp(logres.standard_error),
                        res.standard_error, rtol=rtol, atol=rtol/100)

    @pytest.mark.parametrize("n_points", [2**8, 2**12])
    @pytest.mark.parametrize("n_estimates", [8, 16])
    def test_basic(self, n_points, n_estimates, xp):
        self.basic_test(n_points, n_estimates, xp=xp)

    @pytest.mark.parametrize("signs", [[1., 1.], [-1., -1.], [-1., 1.], [1., -1.]])
    def test_sign(self, signs, xp):
        self.basic_test(signs=signs, xp=xp)

    @pytest.mark.parametrize("log", [False, True])
    def test_zero(self, log, xp):
        message = "A lower limit was equal to an upper limit, so"
        with pytest.warns(UserWarning, match=message):
            res = qmc_quad(lambda x: 1, xp.asarray([0, 0]), xp.asarray([0, 1]), log=log)
        assert res.integral == (-xp.inf if log else 0)
        assert res.standard_error == 0

    def test_flexible_input(self):
        # check that qrng is not required
        # also checks that for 1d problems, a and b can be scalars

Review comment (suggestion):
Looks like a small tolerance bump is needed
```suggestion
        rtol = 1e-14 if res.integral.dtype == xp.float64 else 1.1e-6
```

Code after:
        ndim = 2
        mean = np.zeros(ndim)
        cov = np.eye(ndim)

        def func(x):
            # standard multivariate normal PDF in two dimensions
            return xp.exp(-0.5 * xp.sum(x*x, axis=0)) / (2 * xp.pi)

        rng = np.random.default_rng(2879434385674690281)
        qrng = stats.qmc.Sobol(ndim, seed=rng)
        a = np.zeros(ndim)
        b = np.ones(ndim) * signs
        res = qmc_quad(func, xp.asarray(a, dtype=dtype), xp.asarray(b, dtype=dtype),
                       n_points=n_points, n_estimates=n_estimates, qrng=qrng)
        ref = stats.multivariate_normal.cdf(b, mean, cov, lower_limit=a)
        atol = special.stdtrit(n_estimates-1, 0.995) * res.standard_error  # 99% CI
        xp_assert_close(res.integral, xp.asarray(ref, dtype=dtype), atol=atol)
        assert np.prod(signs)*res.integral > 0

        rng = np.random.default_rng(2879434385674690281)
        qrng = stats.qmc.Sobol(ndim, seed=rng)
        logres = qmc_quad(lambda *args: xp.log(func(*args)),
                          xp.asarray(a, dtype=dtype), xp.asarray(b, dtype=dtype),
                          n_points=n_points, n_estimates=n_estimates,
                          log=True, qrng=qrng)
        rtol = 1e-14 if res.integral.dtype == xp.float64 else 1.1e-6
        xp_assert_close(xp.real(xp.exp(logres.integral)), res.integral, rtol=rtol)
        assert xp.imag(logres.integral + 0j) == (xp.pi if np.prod(signs) < 0 else 0)
        xp_assert_close(xp.exp(logres.standard_error),
                        res.standard_error, rtol=rtol, atol=rtol/100)

    @pytest.mark.parametrize("n_points", [2**8, 2**12])
    @pytest.mark.parametrize("n_estimates", [8, 16])
    def test_basic(self, n_points, n_estimates, xp):
        self.basic_test(n_points, n_estimates, xp=xp)

    @pytest.mark.parametrize("signs", [[1., 1.], [-1., -1.], [-1., 1.], [1., -1.]])
    def test_sign(self, signs, xp):
        self.basic_test(signs=signs, xp=xp)

    @pytest.mark.parametrize("log", [False, True])
    def test_zero(self, log, xp):
        message = "A lower limit was equal to an upper limit, so"
        with pytest.warns(UserWarning, match=message):
            res = qmc_quad(lambda x: 1, xp.asarray([0, 0]), xp.asarray([0, 1]), log=log)
        assert res.integral == (-xp.inf if log else 0)
        assert res.standard_error == 0

    def test_flexible_input(self):
        # check that qrng is not required
        # also checks that for 1d problems, a and b can be scalars