From 9b577deda3dd76c179c5e2c9128a972038b29422 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 30 Jun 2023 11:22:39 -0700 Subject: [PATCH 1/6] 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. --- harmonica/_forward/_tesseroid_utils.py | 114 ++++++++++++++++++------- harmonica/_forward/tesseroid.py | 4 +- harmonica/tests/test_tesseroid.py | 11 +-- 3 files changed, 91 insertions(+), 38 deletions(-) diff --git a/harmonica/_forward/_tesseroid_utils.py b/harmonica/_forward/_tesseroid_utils.py index eaa0b5210..6af26299f 100644 --- a/harmonica/_forward/_tesseroid_utils.py +++ b/harmonica/_forward/_tesseroid_utils.py @@ -8,7 +8,7 @@ Utils functions for tesseroid forward modelling """ import numpy as np -from numba import jit +from numba import jit, prange from numpy.polynomial.legendre import leggauss from .utils import distance_spherical @@ -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,42 +423,94 @@ 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. """ + conflicting = _check_points_outside_tesseroids(coordinates, tesseroids) + if conflicting: + err_msg = ( + "Found computation point(s) inside tesseroid(s). " + "Computation points must be outside of tesseroids.\n" + ) + for i, j in conflicting: + longitude, latitude, radius = coordinates[:, i] + west, east, south, north, bottom, top = tesseroids[j, :] + err_msg += ( + f"Computation point '({longitude}, {latitude}, {radius})' inside " + f"tesseroid '({west}, {east}, {south}, {north}, {bottom}, {top})'. " + ) + raise ValueError(err_msg) + + +@jit(nopython=True) +def _check_points_outside_tesseroids(coordinates, tesseroids): longitude, latitude, radius = coordinates[:] - west, east, south, north, bottom, top = tuple(tesseroids[:, i] for i in range(6)) + conflicting = [] + for i in range(longitude.size): + for j in range(tesseroids.shape[0]): + if _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], + ): + conflicting.append((i, j)) + return conflicting + + +@jit(nopython=True) +def _check_point_inside_tesseroid( + longitude, latitude, radius, west, east, south, north, bottom, top +): + """ + Check if observation point falls inside tesseroid. + + Parameters + ---------- + longitude : float + Longitude coordinate of the observation point in degrees. + latitude : float + Spherical latitude coordinate of the observation point in degrees. + radius : float + Geoscentric spherical radius of the observation point in meters. + west : float + West longitude bound of the tesseroid in degrees. + east : float + East longitude bound of the tesseroid in degrees. + south : float + South latitude bound of the tesseroid in degrees. + north : float + North latitude bound of the tesseroid in degrees. + bottom : float + Bottom radius bound of the tesseroid in meters. + top : float + Top radius bound of the tesseroid in meters. + + Returns + ------- + inside : bool + True if the observation point falls inside the tesseroid. + """ # 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(): - err_msg = ( - "Found computation point inside tesseroid. " - + "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, :] - ) - raise ValueError(err_msg) + if ( + (west < longitude_180 < east or west < longitude_360 < east) + and south < latitude < north + and bottom < radius < top + ): + return True + return False def _longitude_continuity(tesseroids): diff --git a/harmonica/_forward/tesseroid.py b/harmonica/_forward/tesseroid.py index 2714b4def..1dbcbe2c5 100644 --- a/harmonica/_forward/tesseroid.py +++ b/harmonica/_forward/tesseroid.py @@ -13,7 +13,7 @@ from ..constants import GRAVITATIONAL_CONST from ._tesseroid_utils import ( _adaptive_discretization, - _check_points_outside_tesseroids, + check_points_outside_tesseroids, _check_tesseroids, _discard_null_tesseroids, gauss_legendre_quadrature, @@ -155,7 +155,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 diff --git a/harmonica/tests/test_tesseroid.py b/harmonica/tests/test_tesseroid.py index 91b064f16..55cada4fa 100644 --- a/harmonica/tests/test_tesseroid.py +++ b/harmonica/tests/test_tesseroid.py @@ -24,7 +24,7 @@ MAX_DISCRETIZATIONS, STACK_SIZE, _adaptive_discretization, - _check_points_outside_tesseroids, + check_points_outside_tesseroids, _check_tesseroids, tesseroid_gravity, ) @@ -207,6 +207,7 @@ def test_disable_checks(): npt.assert_allclose(invalid_result, -valid_result) +@pytest.mark.use_numba 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]) @@ -223,19 +224,19 @@ def test_point_inside_tesseroid(): np.atleast_2d([0, 10, 150]).T, # point on northern surface ] for coordinates in points: - _check_points_outside_tesseroids(coordinates, tesseroids) + 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) + 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) + 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) + check_points_outside_tesseroids(coordinates, tesseroids) @pytest.mark.use_numba From e58f84ea6d558c331967f4158d5a59ee4eb02b45 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 30 Jun 2023 13:07:58 -0700 Subject: [PATCH 2/6] Run isort --- harmonica/_forward/tesseroid.py | 2 +- harmonica/tests/test_tesseroid.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/harmonica/_forward/tesseroid.py b/harmonica/_forward/tesseroid.py index 1dbcbe2c5..5a6d0ef70 100644 --- a/harmonica/_forward/tesseroid.py +++ b/harmonica/_forward/tesseroid.py @@ -13,9 +13,9 @@ 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, ) diff --git a/harmonica/tests/test_tesseroid.py b/harmonica/tests/test_tesseroid.py index 55cada4fa..29efe7a83 100644 --- a/harmonica/tests/test_tesseroid.py +++ b/harmonica/tests/test_tesseroid.py @@ -24,8 +24,8 @@ MAX_DISCRETIZATIONS, STACK_SIZE, _adaptive_discretization, - check_points_outside_tesseroids, _check_tesseroids, + check_points_outside_tesseroids, tesseroid_gravity, ) from ..constants import GRAVITATIONAL_CONST From 4ef171b063465b8ef89c6a555ed092f6614f3c72 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 30 Jun 2023 15:28:46 -0700 Subject: [PATCH 3/6] Fix issue with multiple tesseroids and points Also, extend the tests for those cases of multiple observation points and tesseroids. --- harmonica/_forward/_tesseroid_utils.py | 9 +-- harmonica/tests/test_tesseroid.py | 97 ++++++++++++++++++-------- 2 files changed, 71 insertions(+), 35 deletions(-) diff --git a/harmonica/_forward/_tesseroid_utils.py b/harmonica/_forward/_tesseroid_utils.py index 6af26299f..5ef132653 100644 --- a/harmonica/_forward/_tesseroid_utils.py +++ b/harmonica/_forward/_tesseroid_utils.py @@ -8,7 +8,7 @@ Utils functions for tesseroid forward modelling """ import numpy as np -from numba import jit, prange +from numba import jit from numpy.polynomial.legendre import leggauss from .utils import distance_spherical @@ -429,6 +429,7 @@ def check_points_outside_tesseroids(coordinates, tesseroids): ValueError If any computation point falls inside any tesseroid. """ + longitude, latitude, radius = coordinates conflicting = _check_points_outside_tesseroids(coordinates, tesseroids) if conflicting: err_msg = ( @@ -436,11 +437,11 @@ def check_points_outside_tesseroids(coordinates, tesseroids): "Computation points must be outside of tesseroids.\n" ) for i, j in conflicting: - longitude, latitude, radius = coordinates[:, i] west, east, south, north, bottom, top = tesseroids[j, :] err_msg += ( - f"Computation point '({longitude}, {latitude}, {radius})' inside " - f"tesseroid '({west}, {east}, {south}, {north}, {bottom}, {top})'. " + f" - Computation point '({longitude[i]}, {latitude[i]}, {radius[i]})' " + "inside tesseroid " + f"'({west}, {east}, {south}, {north}, {bottom}, {top})'.\n" ) raise ValueError(err_msg) diff --git a/harmonica/tests/test_tesseroid.py b/harmonica/tests/test_tesseroid.py index 29efe7a83..760be945d 100644 --- a/harmonica/tests/test_tesseroid.py +++ b/harmonica/tests/test_tesseroid.py @@ -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,36 +207,71 @@ def test_disable_checks(): npt.assert_allclose(invalid_result, -valid_result) -@pytest.mark.use_numba -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).T + 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 observation 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 From 696b6e3483666093c383ab9d537c38341dd2e34c Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 30 Jun 2023 15:34:58 -0700 Subject: [PATCH 4/6] Remove unwanted transpose operation --- harmonica/tests/test_tesseroid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/harmonica/tests/test_tesseroid.py b/harmonica/tests/test_tesseroid.py index 760be945d..e93c405b1 100644 --- a/harmonica/tests/test_tesseroid.py +++ b/harmonica/tests/test_tesseroid.py @@ -249,7 +249,7 @@ def test_point_inside_tesseroid(self, points): 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).T + tesseroid = np.atleast_2d(tesseroid) with pytest.raises(ValueError): check_points_outside_tesseroids(point, tesseroid) From 1755a6ce7cda413143859b93d519d27794e073e6 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 30 Jun 2023 16:23:37 -0700 Subject: [PATCH 5/6] Shorten long line --- harmonica/tests/test_tesseroid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/harmonica/tests/test_tesseroid.py b/harmonica/tests/test_tesseroid.py index e93c405b1..e4ab0711e 100644 --- a/harmonica/tests/test_tesseroid.py +++ b/harmonica/tests/test_tesseroid.py @@ -255,7 +255,7 @@ def test_point_phased_longitude(self, point, tesseroid): @pytest.mark.use_numba def test_multiple_points_and_tesseroids(self): - """Check if error is raised with multiple observation points and tesseroids""" + """Check if error is raised with multiple points and tesseroids.""" tesseroids = np.atleast_2d( [ [-10, 10, -10, 10, 100, 200], From d015add0fb17304109acaddc758e911b969916f7 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Wed, 5 Jul 2023 13:24:58 -0700 Subject: [PATCH 6/6] Merge private check functions into a single one Include the if statement in the function that runs the for loop over observation points and tesseroids. --- harmonica/_forward/_tesseroid_utils.py | 70 +++++--------------------- 1 file changed, 13 insertions(+), 57 deletions(-) diff --git a/harmonica/_forward/_tesseroid_utils.py b/harmonica/_forward/_tesseroid_utils.py index 5ef132653..2d73bdb07 100644 --- a/harmonica/_forward/_tesseroid_utils.py +++ b/harmonica/_forward/_tesseroid_utils.py @@ -448,72 +448,28 @@ def check_points_outside_tesseroids(coordinates, tesseroids): @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]): - if _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], + # 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 -@jit(nopython=True) -def _check_point_inside_tesseroid( - longitude, latitude, radius, west, east, south, north, bottom, top -): - """ - Check if observation point falls inside tesseroid. - - Parameters - ---------- - longitude : float - Longitude coordinate of the observation point in degrees. - latitude : float - Spherical latitude coordinate of the observation point in degrees. - radius : float - Geoscentric spherical radius of the observation point in meters. - west : float - West longitude bound of the tesseroid in degrees. - east : float - East longitude bound of the tesseroid in degrees. - south : float - South latitude bound of the tesseroid in degrees. - north : float - North latitude bound of the tesseroid in degrees. - bottom : float - Bottom radius bound of the tesseroid in meters. - top : float - Top radius bound of the tesseroid in meters. - - Returns - ------- - inside : bool - True if the observation point falls inside the tesseroid. - """ - # 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 - if ( - (west < longitude_180 < east or west < longitude_360 < east) - and south < latitude < north - and bottom < radius < top - ): - return True - return False - - def _longitude_continuity(tesseroids): """ Modify longitudinal boundaries of tesseroids to ensure longitude continuity