Repository: fatiando/harmonica (python)

Issue:
Check for observartion point inside tesseroid eats too much memory
**Description of the problem:**

The tesseroid forward modelling function checks if no observation point falls inside any tesseroid through the  `_check_points_outside_tesseroids` function:
https://github.com/fatiando/harmonica/blob/b8f037e6d72b04cb0b4ff62bb12b63dd13e38808/harmonica/_forward/_tesseroid_utils.py#L405-L461

This functions makes use of Numpy arrays to run these checks, but this is super inefficient, since Numpy will try to allocate massive arrays of shapes (`n_coords`, `n_tesseroids`).

**Solution**
We should rewrite this function using Numba instead, making sure that we aren't allocating any additional array and the function just checks if observation points don't fall inside any tesseroids. We could also benefit from Numba parallelizations and run the first for loop in parallel (this should follow the user's choice set through the `parallel` argument in `tesseroid_gravity`).

As a quick design, we could do something like:
```python

@numba.jit(nopython=True, parallel=True)
def check_points_outside_tesseroids(coordinates, tesseroids):
    """..."""
    longitude, latitude, radius = coordinates[:]
    for i in numba.prange(longitude.size):
        for j in range(tesseroids.shape[0]):
            _check_point_inside_tesseroid(
                longitude[i],
                latitude[i],
                radius[i],
                tesseroids[j, 0],
                tesseroids[j, 1],
                tesseroids[j, 2],
                tesseroids[j, 3],
                tesseroids[j, 4],
                tesseroids[j, 5],
            )
            
@numba.jit(nopython=True)
def _check_point_inside_tesseroid(
    longitude, latitude, radius, west, east, south, north, bottom, top
):
    """Raise error if observation point falls inside tesseroid"""
    ...
    if ... :
        raise ValueError(...)
```

Check out https://github.com/fatiando/harmonica/blob/main/harmonica/_forward/prism.py for inspiration on how to code this.

**Temporary solution**

For the moment, the check can be bypassed by passing `disable_checks=True`, but it also disables other checks.

Pull request:
Rewrite check for point inside tesseroid with Numba
Rewrite the functions that checks if a computation point falls inside a tesseroid using Numba instad of Numpy. This greatly reduce the memory need of the Numpy functions, that required to build large matrices even for medium-size problems.


**Relevant issues/PRs:**

Fixes #417


Patch:
diff --git a/harmonica/_forward/_tesseroid_utils.py b/harmonica/_forward/_tesseroid_utils.py
index eaa0b52..2d73bdb 100644
--- a/harmonica/_forward/_tesseroid_utils.py
+++ b/harmonica/_forward/_tesseroid_utils.py
@@ -402,7 +402,7 @@ def _check_tesseroids(tesseroids):
     return tesseroids
 
 
-def _check_points_outside_tesseroids(coordinates, tesseroids):
+def check_points_outside_tesseroids(coordinates, tesseroids):
     """
     Check if computation points are not inside the tesseroids
 
@@ -423,44 +423,53 @@ def _check_points_outside_tesseroids(coordinates, tesseroids):
         This array of tesseroids must have longitude continuity and valid
         boundaries.
         Run ``_check_tesseroids`` before.
+
+    Raises
+    ------
+    ValueError
+        If any computation point falls inside any tesseroid.
     """
-    longitude, latitude, radius = coordinates[:]
-    west, east, south, north, bottom, top = tuple(tesseroids[:, i] for i in range(6))
-    # Longitudinal boundaries of the tesseroid must be compared with
-    # longitudinal coordinates of computation points when moved to
-    # [0, 360) and [-180, 180).
-    longitude_360 = longitude % 360
-    longitude_180 = ((longitude + 180) % 360) - 180
-    inside_longitude = np.logical_or(
-        np.logical_and(
-            west < longitude_360[:, np.newaxis], longitude_360[:, np.newaxis] < east
-        ),
-        np.logical_and(
-            west < longitude_180[:, np.newaxis], longitude_180[:, np.newaxis] < east
-        ),
-    )
-    inside_latitude = np.logical_and(
-        south < latitude[:, np.newaxis], latitude[:, np.newaxis] < north
-    )
-    inside_radius = np.logical_and(
-        bottom < radius[:, np.newaxis], radius[:, np.newaxis] < top
-    )
-    # Build array of booleans.
-    # The (i, j) element is True if the computation point i is inside the
-    # tesseroid j.
-    inside = inside_longitude * inside_latitude * inside_radius
-    if inside.any():
+    longitude, latitude, radius = coordinates
+    conflicting = _check_points_outside_tesseroids(coordinates, tesseroids)
+    if conflicting:
         err_msg = (
-            "Found computation point inside tesseroid. "
-            + "Computation points must be outside of tesseroids.\n"
+            "Found computation point(s) inside tesseroid(s). "
+            "Computation points must be outside of tesseroids.\n"
         )
-        for point_i, tess_i in np.argwhere(inside):
-            err_msg += "\tComputation point '{}' found inside tesseroid '{}'\n".format(
-                coordinates[:, point_i], tesseroids[tess_i, :]
+        for i, j in conflicting:
+            west, east, south, north, bottom, top = tesseroids[j, :]
+            err_msg += (
+                f" - Computation point '({longitude[i]}, {latitude[i]}, {radius[i]})' "
+                "inside tesseroid "
+                f"'({west}, {east}, {south}, {north}, {bottom}, {top})'.\n"
             )
         raise ValueError(err_msg)
 
 
+@jit(nopython=True)
+def _check_points_outside_tesseroids(coordinates, tesseroids):
+    """
+    Check if observation points fall inside tesseroids.
+    """
+    longitude, latitude, radius = coordinates[:]
+    conflicting = []
+    for i in range(longitude.size):
+        for j in range(tesseroids.shape[0]):
+            # Longitudinal boundaries of the tesseroid must be compared with
+            # longitudinal coordinates of computation points when moved to
+            # [0, 360) and [-180, 180).
+            longitude_360 = longitude[i] % 360
+            longitude_180 = ((longitude[i] + 180) % 360) - 180
+            west, east, south, north, bottom, top = tesseroids[j, :]
+            if (
+                (west < longitude_180 < east or west < longitude_360 < east)
+                and south < latitude[i] < north
+                and bottom < radius[i] < top
+            ):
+                conflicting.append((i, j))
+    return conflicting
+
+
 def _longitude_continuity(tesseroids):
     """
     Modify longitudinal boundaries of tesseroids to ensure longitude continuity
diff --git a/harmonica/_forward/tesseroid.py b/harmonica/_forward/tesseroid.py
index 4154f2e..3472695 100644
--- a/harmonica/_forward/tesseroid.py
+++ b/harmonica/_forward/tesseroid.py
@@ -13,9 +13,9 @@ from numba import jit, prange
 from ..constants import GRAVITATIONAL_CONST
 from ._tesseroid_utils import (
     _adaptive_discretization,
-    _check_points_outside_tesseroids,
     _check_tesseroids,
     _discard_null_tesseroids,
+    check_points_outside_tesseroids,
     gauss_legendre_quadrature,
     glq_nodes_weights,
 )
@@ -163,7 +163,7 @@ def tesseroid_gravity(
     # Sanity checks for tesseroids and computation points
     if not disable_checks:
         tesseroids = _check_tesseroids(tesseroids)
-        _check_points_outside_tesseroids(coordinates, tesseroids)
+        check_points_outside_tesseroids(coordinates, tesseroids)
     # Check if density is a function or constant values
     if callable(density):
         # Run density-based discretization on each tesseroid


Test patch:
diff --git a/harmonica/tests/test_tesseroid.py b/harmonica/tests/test_tesseroid.py
index 91b064f..e4ab071 100644
--- a/harmonica/tests/test_tesseroid.py
+++ b/harmonica/tests/test_tesseroid.py
@@ -24,8 +24,8 @@ from .._forward.tesseroid import (
     MAX_DISCRETIZATIONS,
     STACK_SIZE,
     _adaptive_discretization,
-    _check_points_outside_tesseroids,
     _check_tesseroids,
+    check_points_outside_tesseroids,
     tesseroid_gravity,
 )
 from ..constants import GRAVITATIONAL_CONST
@@ -129,7 +129,7 @@ def test_valid_tesseroid():
     _check_tesseroids(np.atleast_2d([w, e, s, n, bottom, bottom]))
     # Check if valid tesseroid with west > east is not caught
     _check_tesseroids(np.atleast_2d([350, 10, s, n, bottom, top]))
-    # Check if tesseroids on [-180, 180) are not caught
+    # Check if tesseroids on -180, 180 are not caught
     _check_tesseroids(np.atleast_2d([-70, -60, s, n, bottom, top]))
     _check_tesseroids(np.atleast_2d([-10, 10, s, n, bottom, top]))
     _check_tesseroids(np.atleast_2d([-150, 150, s, n, bottom, top]))
@@ -207,35 +207,71 @@ def test_disable_checks():
     npt.assert_allclose(invalid_result, -valid_result)
 
 
-def test_point_inside_tesseroid():
-    "Check if a computation point inside the tesseroid is caught"
-    tesseroids = np.atleast_2d([-10, 10, -10, 10, 100, 200])
-    # Test if outside point is not caught
-    points = [
-        np.atleast_2d([0, 0, 250]).T,  # outside point on radius
-        np.atleast_2d([20, 0, 150]).T,  # outside point on longitude
-        np.atleast_2d([0, 20, 150]).T,  # outside point on latitude
-        np.atleast_2d([0, 0, 200]).T,  # point on top surface
-        np.atleast_2d([0, 0, 100]).T,  # point on bottom surface
-        np.atleast_2d([-10, 0, 150]).T,  # point on western surface
-        np.atleast_2d([10, 0, 150]).T,  # point on eastern surface
-        np.atleast_2d([0, -10, 150]).T,  # point on southern surface
-        np.atleast_2d([0, 10, 150]).T,  # point on northern surface
-    ]
-    for coordinates in points:
-        _check_points_outside_tesseroids(coordinates, tesseroids)
-    # Test if computation point is inside the tesseroid
-    coordinates = np.atleast_2d([0, 0, 150]).T
-    with pytest.raises(ValueError):
-        _check_points_outside_tesseroids(coordinates, tesseroids)
-    # Test if computation point with phased longitude is inside the tesseroid
-    coordinates = np.atleast_2d([360, 0, 150]).T
-    with pytest.raises(ValueError):
-        _check_points_outside_tesseroids(coordinates, tesseroids)
-    tesseroids = np.atleast_2d([260, 280, -10, 10, 100, 200])
-    coordinates = np.atleast_2d([-90, 0, 150]).T
-    with pytest.raises(ValueError):
-        _check_points_outside_tesseroids(coordinates, tesseroids)
+class TestPointInsideTesseroid:
+    """Test checks for points inside tesseroids."""
+
+    @pytest.mark.use_numba
+    def test_points_outside_tesseroid(self):
+        """Check if error is not raised when no point is inside a tesseroid."""
+        points = np.array(
+            [
+                [0, 0, 250],  # outside point on radius
+                [20, 0, 150],  # outside point on longitude
+                [0, 20, 150],  # outside point on latitude
+                [0, 0, 200],  # point on top surface
+                [0, 0, 100],  # point on bottom surface
+                [-10, 0, 150],  # point on western surface
+                [10, 0, 150],  # point on eastern surface
+                [0, -10, 150],  # point on southern surface
+                [0, 10, 150],  # point on northern surface
+            ]
+        ).T
+        tesseroid = np.atleast_2d([-10, 10, -10, 10, 100, 200])
+        check_points_outside_tesseroids(points, tesseroid)
+
+    @pytest.mark.use_numba
+    @pytest.mark.parametrize("points", [[0, 0, 150], [360, 0, 150]])
+    def test_point_inside_tesseroid(self, points):
+        """Check if error is raised when point fall inside the tesseroid."""
+        points = np.atleast_2d(points).T
+        tesseroid = np.atleast_2d([-10, 10, -10, 10, 100, 200])
+        with pytest.raises(ValueError):
+            check_points_outside_tesseroids(points, tesseroid)
+
+    @pytest.mark.use_numba
+    @pytest.mark.parametrize(
+        "point, tesseroid",
+        [
+            ([360, 0, 150], [-10, 10, -10, 10, 100, 200]),
+            ([-90, 0, 150], [260, 280, -10, 10, 100, 200]),
+        ],
+    )
+    def test_point_phased_longitude(self, point, tesseroid):
+        """Test if error is raised when the longitude coords are phased."""
+        point = np.atleast_2d(point).T
+        tesseroid = np.atleast_2d(tesseroid)
+        with pytest.raises(ValueError):
+            check_points_outside_tesseroids(point, tesseroid)
+
+    @pytest.mark.use_numba
+    def test_multiple_points_and_tesseroids(self):
+        """Check if error is raised with multiple points and tesseroids."""
+        tesseroids = np.atleast_2d(
+            [
+                [-10, 10, -10, 10, 100, 200],
+                [20, 30, 20, 30, 400, 500],
+                [-50, -40, -30, -20, 100, 500],
+            ]
+        )
+        points = np.array(
+            [
+                [0, 0, 150],
+                [80, 82, 4000],
+                [10, 10, 450],
+            ]
+        ).T
+        with pytest.raises(ValueError):
+            check_points_outside_tesseroids(points, tesseroids)
 
 
 @pytest.mark.use_numba