Repository: scipy/scipy
File: scipy/interpolate/tests/test_interpnd.py (Python)

Code before:

        worker_thread_1 = threading.Thread(
            target=do_interp,
            args=(interpolator, slice(0, 2100), slice(0, 2100)))
        worker_thread_2 = threading.Thread(
            target=do_interp,
            args=(interpolator, slice(2100, 4200), slice(0, 2100)))
        worker_thread_3 = threading.Thread(
            target=do_interp,
            args=(interpolator, slice(0, 2100), slice(2100, 4200)))
        worker_thread_4 = threading.Thread(
            target=do_interp,
            args=(interpolator, slice(2100, 4200), slice(2100, 4200)))

        worker_thread_1.start()
        worker_thread_2.start()
        worker_thread_3.start()
        worker_thread_4.start()

        worker_thread_1.join()
        worker_thread_2.join()
        worker_thread_3.join()
        worker_thread_4.join()

    @pytest.mark.parametrize("simplex_tolerance", [1, 10])
    def test_interp_from_boundary(self, simplex_tolerance):
        coords = np.block(
            [[0, np.arange(400, 500), np.zeros(100), 800],
             [0, np.zeros(100), np.arange(400, 500), 800]]
        )
        values = np.ones(202)
        interpolator = interpnd.LinearNDInterpolator(coords.T, values)
        results = interpolator(
            np.mgrid[:500, :500].T, simplex_tolerance=simplex_tolerance
        )
        num_nans = np.count_nonzero(np.isnan(results))
        if simplex_tolerance == 1:
            # The default tolerance doesn't work when filling from a
            # dense boundary
            assert num_nans > 0
        else:
            # The larger tolerance allows QHull to assign a simplex
            # for each point
            assert num_nans == 0


class TestEstimateGradients2DGlobal:
    def test_smoketest(self):
        x = np.array([(0, 0), (0, 2),
                      (1, 0), (1, 2), (0.25, 0.75), (0.6, 0.8)], dtype=float)
        tri = qhull.Delaunay(x)

Review comment (suggestion):
```suggestion
    def test_interp_from_boundary(self, simplex_tolerance):
         # regression test for https://github.com/scipy/scipy/issues/21910
```

Code after:

        worker_thread_1 = threading.Thread(
            target=do_interp,
            args=(interpolator, slice(0, 2100), slice(0, 2100)))
        worker_thread_2 = threading.Thread(
            target=do_interp,
            args=(interpolator, slice(2100, 4200), slice(0, 2100)))
        worker_thread_3 = threading.Thread(
            target=do_interp,
            args=(interpolator, slice(0, 2100), slice(2100, 4200)))
        worker_thread_4 = threading.Thread(
            target=do_interp,
            args=(interpolator, slice(2100, 4200), slice(2100, 4200)))

        worker_thread_1.start()
        worker_thread_2.start()
        worker_thread_3.start()
        worker_thread_4.start()

        worker_thread_1.join()
        worker_thread_2.join()
        worker_thread_3.join()
        worker_thread_4.join()


class TestEstimateGradients2DGlobal:
    def test_smoketest(self):
        x = np.array([(0, 0), (0, 2),
                      (1, 0), (1, 2), (0.25, 0.75), (0.6, 0.8)], dtype=float)
        tri = qhull.Delaunay(x)

        # Should be exact for linear functions, independent of triangulation

        funcs = [
            (lambda x, y: 0*x + 1, (0, 0)),
            (lambda x, y: 0 + x, (1, 0)),
            (lambda x, y: -2 + y, (0, 1)),
            (lambda x, y: 3 + 3*x + 14.15*y, (3, 14.15))
        ]

        for j, (func, grad) in enumerate(funcs):
            z = func(x[:,0], x[:,1])
            dz = interpnd.estimate_gradients_2d_global(tri, z, tol=1e-6)

            assert dz.shape == (6, 2)
            xp_assert_close(
                dz, np.array(grad)[None, :] + 0*dz, rtol=1e-5, atol=1e-5, 
                err_msg=f"item {j}"
            )

    def test_regression_2359(self):