From 826cf4b583a287ab42d3f9540f8b0ea4c74f4afc Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Thu, 14 Nov 2024 11:27:05 -0800 Subject: [PATCH] Make 3rd-order prism kernels return np.nan on vertices Make the 3rd order kernel functions for prisms to return `np.nan` when the observation points matches any vertex of the prisms. Fixes division-by-zero error when evaluating `kernel_enu` when `radius` is zero, and makes the behaviour consistent with the other third-order kernels. --- choclo/prism/_kernels.py | 40 ++++++++++++++++-------- choclo/tests/test_prism_kernels.py | 49 ++++++++++++++++++++++++++++++ 2 files changed, 77 insertions(+), 12 deletions(-) create mode 100644 choclo/tests/test_prism_kernels.py diff --git a/choclo/prism/_kernels.py b/choclo/prism/_kernels.py index b3f8f845..bff6ab65 100644 --- a/choclo/prism/_kernels.py +++ b/choclo/prism/_kernels.py @@ -932,7 +932,8 @@ def kernel_eee(easting, northing, upward, radius): kernel : float Value of the kernel function for the easting-easting-easting component of the grad tensor due to a rectangular prism evaluated on a single - vertex. + vertex. Return ``np.nan`` if ``radius`` is zero (i.e. observation point + is on the vertex). Notes ----- @@ -990,7 +991,8 @@ def kernel_nnn(easting, northing, upward, radius): kernel : float Value of the kernel function for the northing-northing-northing component of the grad tensor due to a rectangular prism evaluated on - a single vertex. + a single vertex. Return ``np.nan`` if ``radius`` is zero (i.e. + observation point is on the vertex). Notes ----- @@ -1048,7 +1050,8 @@ def kernel_uuu(easting, northing, upward, radius): kernel : float Value of the kernel function for the upward-upward-upward component of the grad tensor due to a rectangular prism evaluated on a single - vertex. + vertex. Return ``np.nan`` if ``radius`` is zero (i.e. observation point + is on the vertex). Notes ----- @@ -1106,7 +1109,8 @@ def kernel_een(easting, northing, upward, radius): kernel : float Value of the kernel function for the easting-easting-northing component of the grad tensor due to a rectangular prism evaluated on a single - vertex. + vertex. Return ``np.nan`` if ``radius`` is zero (i.e. observation point + is on the vertex). Notes ----- @@ -1164,7 +1168,8 @@ def kernel_eeu(easting, northing, upward, radius): kernel : float Value of the kernel function for the easting-easting-upward component of the grad tensor due to a rectangular prism evaluated on a single - vertex. + vertex. Return ``np.nan`` if ``radius`` is zero (i.e. observation point + is on the vertex). Notes ----- @@ -1222,7 +1227,8 @@ def kernel_enn(easting, northing, upward, radius): kernel : float Value of the kernel function for the easting-northing-northing component of the grad tensor due to a rectangular prism evaluated on - a single vertex. + a single vertex. Return ``np.nan`` if ``radius`` is zero (i.e. + observation point is on the vertex). Notes ----- @@ -1280,7 +1286,8 @@ def kernel_nnu(easting, northing, upward, radius): kernel : float Value of the kernel function for the northing-northing-upward component of the grad tensor due to a rectangular prism evaluated on - a single vertex. + a single vertex. Return ``np.nan`` if ``radius`` is zero (i.e. + observation point is on the vertex). Notes ----- @@ -1338,7 +1345,8 @@ def kernel_euu(easting, northing, upward, radius): kernel : float Value of the kernel function for the easting-upward-upward component of the grad tensor due to a rectangular prism evaluated on - a single vertex. + a single vertex. Return ``np.nan`` if ``radius`` is zero (i.e. + observation point is on the vertex). Notes ----- @@ -1396,7 +1404,8 @@ def kernel_nuu(easting, northing, upward, radius): kernel : float Value of the kernel function for the northing-upward-upward component of the grad tensor due to a rectangular prism evaluated on - a single vertex. + a single vertex. Return ``np.nan`` if ``radius`` is zero (i.e. + observation point is on the vertex). Notes ----- @@ -1454,7 +1463,8 @@ def kernel_enu(easting, northing, upward, radius): kernel : float Value of the kernel function for the easting-northing-upward component of the grad tensor due to a rectangular prism evaluated on a single - vertex. + vertex. Return ``np.nan`` if ``radius`` is zero (i.e. observation point + is on the vertex). Notes ----- @@ -1472,6 +1482,8 @@ def kernel_enu(easting, northing, upward, radius): >>> float(kernel_enu(x, y, z, r)) # doctest: +NUMBER -0.1480061 """ + if radius == 0.0: + return np.nan return -1 / radius @@ -1507,7 +1519,7 @@ def _kernel_iii(x_i, x_j, x_k, radius): Returns ------- kernel : float - Value of the kernel function. + Value of the kernel function. Return ``np.nan`` if ``radius`` is zero. Examples -------- @@ -1523,6 +1535,8 @@ def _kernel_iii(x_i, x_j, x_k, radius): >>> k_nnn = _kernel_iii(northing, upward, easting, radius) >>> k_uuu = _kernel_iii(upward, easting, northing, radius) """ + if radius == 0.0: + return np.nan if (x_i == 0 and x_j == 0) or (x_i == 0 and x_k == 0): return 0.0 x_i_sq, x_j_sq, x_k_sq = x_i**2, x_j**2, x_k**2 @@ -1567,7 +1581,7 @@ def _kernel_iij(x_i, x_j, x_k, radius): Returns ------- kernel : float - Value of the kernel function. + Value of the kernel function. Return ``np.nan`` if ``radius`` is zero. Examples -------- @@ -1586,6 +1600,8 @@ def _kernel_iij(x_i, x_j, x_k, radius): >>> k_euu = _kernel_iij(upward, easting, northing, radius) >>> k_nuu = _kernel_iij(upward, northing, upward, radius) """ + if radius == 0.0: + return np.nan if x_i == 0 and x_j == 0: return 0.0 return -x_i / ((x_k + radius) * radius) diff --git a/choclo/tests/test_prism_kernels.py b/choclo/tests/test_prism_kernels.py new file mode 100644 index 00000000..b09d02ee --- /dev/null +++ b/choclo/tests/test_prism_kernels.py @@ -0,0 +1,49 @@ +""" +Additional tests to prism kernels. +""" + +import pytest +import numpy as np + +from choclo.prism import ( + kernel_eee, + kernel_nnn, + kernel_uuu, + kernel_een, + kernel_eeu, + kernel_enn, + kernel_nnu, + kernel_euu, + kernel_nuu, + kernel_enu, +) + + +class TestThirdOrderKernelsOnVertex: + """ + Test if third-order kernels evaluated on a vertex return np.nan. + + This functionality is meant to be supported for the public usage of this + kernels, although the forward modelling functions within Choclo won't use + it since they filter out singular points before they call kernels. + """ + + KERNELS = ( + kernel_eee, + kernel_nnn, + kernel_uuu, + kernel_een, + kernel_eeu, + kernel_enn, + kernel_nnu, + kernel_euu, + kernel_nuu, + kernel_enu, + ) + + @pytest.mark.parametrize("kernel", KERNELS) + def test_third_order_kernels_on_vertex(self, kernel): + easting, northing, upward = 0.0, 0.0, 0.0 + radius = 0.0 + result = kernel(easting, northing, upward, radius) + assert np.isnan(result)