Repository: cta-observatory/pyirf

Issue:
Nearest Simplex finding in ParametrizedNearestSimplexExtrapolator not working
Apparently, `scipy.spatial.Delaunay.plane_distance` does not work in all cases to find the nearest simplex to some point outside the grid's convex hull. For non-rectangular grids this does somewhat work
![non_rect_scipy](https://github.com/cta-observatory/pyirf/assets/22680504/3202575a-baf4-4c8b-9e98-171e7d4f8cf9)
but it fails completely for rectangular grids
![rect_scipy](https://github.com/cta-observatory/pyirf/assets/22680504/816f6c2e-8ed1-49b0-8833-7bb169c1b338)

With some assumptions the problem can be solved geometrically to behave as expected:
![non_rect_geom](https://github.com/cta-observatory/pyirf/assets/22680504/da4775ee-f627-4178-b993-3d36a234c2b2)
![rect_geom](https://github.com/cta-observatory/pyirf/assets/22680504/61504474-a3fd-445f-972d-22badc59ebd4)


I'll create a PR implementing the solution used to create above plots, some of these functions could even be reused for later extrapolation methods (namely blending over visible edges). I'd like, however, some advise where to put the functions used for this as they are quite general. The first place I have in mind would ofc be just as free functions inside `pyirf.interpolation.nearest_simplex_extrapolator` as they are first needed there, but since they are so general and can partially be reused later, maybe we want something as `pyirf.interpolation.utils` @maxnoe?

As of now, the `ParametrizedNearestSimplexExtrapolator` should not be used. As the same approach was used in #237 this issue also blocks this PR.




Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/cta-observatory/pyirf/pull/240?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) Report
Patch coverage: **`100.00`**% and project coverage change: **`+0.16`** :tada:
> Comparison is base [(`642940a`)](https://app.codecov.io/gh/cta-observatory/pyirf/commit/642940afda460c8eeb5624b5c220bc9d8753db45?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) 94.33% compared to head [(`5b9fe86`)](https://app.codecov.io/gh/cta-observatory/pyirf/pull/240?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) 94.50%.

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main     #240      +/-   ##
==========================================
+ Coverage   94.33%   94.50%   +0.16%     
==========================================
  Files          58       60       +2     
  Lines        2738     2821      +83     
==========================================
+ Hits         2583     2666      +83     
  Misses        155      155              
```


| [Impacted Files](https://app.codecov.io/gh/cta-observatory/pyirf/pull/240?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) | Coverage Δ | |
|---|---|---|
| [...yirf/interpolation/nearest\_simplex\_extrapolator.py](https://app.codecov.io/gh/cta-observatory/pyirf/pull/240?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-cHlpcmYvaW50ZXJwb2xhdGlvbi9uZWFyZXN0X3NpbXBsZXhfZXh0cmFwb2xhdG9yLnB5) | `100.00% <100.00%> (ø)` | |
| [.../interpolation/tests/test\_griddata\_interpolator.py](https://app.codecov.io/gh/cta-observatory/pyirf/pull/240?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-cHlpcmYvaW50ZXJwb2xhdGlvbi90ZXN0cy90ZXN0X2dyaWRkYXRhX2ludGVycG9sYXRvci5weQ==) | `100.00% <100.00%> (ø)` | |
| [pyirf/interpolation/tests/test\_utils.py](https://app.codecov.io/gh/cta-observatory/pyirf/pull/240?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-cHlpcmYvaW50ZXJwb2xhdGlvbi90ZXN0cy90ZXN0X3V0aWxzLnB5) | `100.00% <100.00%> (ø)` | |
| [pyirf/interpolation/utils.py](https://app.codecov.io/gh/cta-observatory/pyirf/pull/240?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory#diff-cHlpcmYvaW50ZXJwb2xhdGlvbi91dGlscy5weQ==) | `100.00% <100.00%> (ø)` | |


</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/cta-observatory/pyirf/pull/240?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory).   
:loudspeaker: Do you have feedback about the report comment? [Let us know in this issue](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory).


Patch:
diff --git a/docs/changes/240.bugfix.rst b/docs/changes/240.bugfix.rst
new file mode 100644
index 0000000..e16b615
--- /dev/null
+++ b/docs/changes/240.bugfix.rst
@@ -0,0 +1,1 @@
+Fixes a bug, where the nearest simplex in NearestSimplexExtrapolator was not computed correct on a 2D grid.
diff --git a/pyirf/interpolation/nearest_simplex_extrapolator.py b/pyirf/interpolation/nearest_simplex_extrapolator.py
index ec5e473..70bc236 100644
--- a/pyirf/interpolation/nearest_simplex_extrapolator.py
+++ b/pyirf/interpolation/nearest_simplex_extrapolator.py
@@ -10,6 +10,7 @@ from .moment_morph_interpolator import (
     barycentric_2D_interpolation_coefficients,
     linesegment_1D_interpolation_coefficients,
 )
+from .utils import find_nearest_simplex
 
 __all__ = ["ParametrizedNearestSimplexExtrapolator"]
 
@@ -103,8 +104,9 @@ class ParametrizedNearestSimplexExtrapolator(ParametrizedExtrapolator):
 
             extrapolant = self._extrapolate1D(segment_inds, target_point)
         elif self.grid_dim == 2:
-            dists = self.triangulation.plane_distance(target_point)
-            nearest_simplex_ind = np.argmax(dists)
+            nearest_simplex_ind = find_nearest_simplex(
+                self.triangulation, target_point.squeeze()
+            )
             simplex_indices = self.triangulation.simplices[nearest_simplex_ind]
 
             extrapolant = self._extrapolate2D(simplex_indices, target_point)
diff --git a/pyirf/interpolation/utils.py b/pyirf/interpolation/utils.py
new file mode 100644
index 0000000..aa54ad1
--- /dev/null
+++ b/pyirf/interpolation/utils.py
@@ -0,0 +1,167 @@
+import numpy as np
+
+
+def plumb_point_dist(line, target):
+    """
+    Compute minimal distance between target and line under the constraint, that it has
+    to lay between the points building the line and not on the extension of it.
+
+    Parameters
+    ----------
+    line: np.ndarray, shape=(2, M)
+        Array of two points spanning a line segment. Might be in two or three dims M.
+    target: np.ndarray, shape=(M)
+        Target point, of which the minimal distance to line segement is needed
+
+    Returns
+    -------
+    d_min: float
+        Minimal distance to line segement between points in line
+    """
+    A = line[0]
+    B = line[1]
+    P = target
+
+    # Costruct the footpoint/plumb point of the target projected onto
+    # both lines OA + r1*AB (F1) and OB + r1*BA (F2), for details see
+    # https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Vector_formulation
+    F1 = A + np.dot(P - A, B - A) * (B - A) / np.dot(B - A, B - A)
+    F2 = B + np.dot(P - B, A - B) * (A - B) / np.dot(A - B, A - B)
+
+    # Find, at which parameter value r1/r2 the plumb point lies on line F1/F2
+    if B[0] - A[0] == 0:
+        r1 = (F1[1] - A[1]) / (B[1] - A[1])
+    else:
+        r1 = (F1[0] - A[0]) / (B[0] - A[0])
+
+    if A[0] - B[0] == 0:
+        r2 = (F2[1] - B[1]) / (A[1] - B[1])
+    else:
+        r2 = (F2[0] - B[0]) / (A[0] - B[0])
+
+    # If |r1| + |r2| == 1, the plomb point lies between A and B, thus the
+    # distance is the seareched one. In cases where the plumb point is A or B
+    # use the second method consistently, as there might be +/- eps differences
+    # due to different methods of computation.
+    if np.isclose(np.abs(r1) + np.abs(r2), 1) and not (
+        np.isclose(r1, 0) or np.isclose(r2, 0)
+    ):
+        # Compute distance of plumb point to line, for details see
+        # https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Another_vector_formulation
+        return np.linalg.norm(np.cross(P - A, B - A)) / np.linalg.norm(B - A)
+    # If not, the nearest point A <= x <= B is one A and B, thus the searched distance
+    # is the one to this nearest point
+    else:
+        return np.min(np.array([np.linalg.norm(A - P), np.linalg.norm(B - P)]))
+
+
+def point_facet_angle(line, target):
+    """
+    Compute cos(angle) between target and a line segment"
+
+    Parameters
+    ----------
+    line: np.ndarray, shape=(2, M)
+        Array of two points spanning a line segment. Might be in two or three dims M.
+    target: np.ndarray, shape=(M)
+        Target point, of which the angle is needed.
+
+    Returns
+    -------
+    cos_angle: float
+        Cosine of angle at target in the triangle ABT with line-points A and B and target T.
+    """
+    PB = line[1] - target
+    PA = line[0] - target
+
+    # For details see https://en.wikipedia.org/wiki/Angle#Dot_product_and_generalisations
+    return np.dot(PB, PA) / (np.linalg.norm(PA) * np.linalg.norm(PB))
+
+
+def find_nearest_facet(qhull_points, target):
+    """
+    Search nearest facet by looking for the closest point on a facet. If this point is an edge point
+    (which by definition lays on two facets) use the one with the lowest cos of angle
+    (maximising angle) as fallback
+
+    Parameters
+    ----------
+    qhull_points: np.ndarray, shape=(n_facets, 2, 2)
+        Array containting both points building a facet for all facets in a
+        grids convex hull.
+    target: np.ndarray, shape=(2)
+        Target point, of which the nearest facet on the convex hull is wanted.
+
+    Returns
+    -------
+    nearest_facet_ind: int
+        Index of the nearest facet in qhull_points.
+    """
+    plumbs = np.array(
+        [
+            (plumb_point_dist(line, target), point_facet_angle(line, target))
+            for line in qhull_points
+        ],
+        dtype=[("plumb_dist", "<f16"), ("facet_point_angle", "<f16")],
+    )
+
+    return np.argsort(plumbs, order=["plumb_dist", "facet_point_angle"])[0]
+
+
+def find_simplex_to_facet(simplices_points, facet_points):
+    """
+    Find simplex of triangulation corresponding to one facet of the convex hull
+    utilizing that each facet is only part of one triangulation simplex
+
+    Parameters
+    ----------
+    simplices_points: np.ndarray, shape=(n_simplices, 3, 2)
+        Array containting all three points building a triangluation simplex for all
+        simplices in a grid's triangulation.
+    facet_points: np.ndarray, shape=(2, 2)
+        Points of which the containing simplex is wanted.
+
+    Returns
+    -------
+    simplex_ind: int
+        Index of the simplex corresponding to a facet spanned by facet_points.
+    """
+    # A qhull facet can only be part of one simplex, both end-points have to be part
+    # of the points spanning the simplex
+    lookup = np.logical_or(
+        np.all(simplices_points == facet_points[0], axis=-1),
+        np.all(simplices_points == facet_points[1], axis=-1),
+    )
+
+    # There is only one simplex, where both points are part of, thus sum(lookup) = 2
+    return np.argmax(np.sum(lookup, axis=-1))
+
+
+def find_nearest_simplex(triangulation, target):
+    """
+    Find nearest simplex to target by finding the nearest convex hull facet and
+    the corresponding simplex
+
+    Parameters
+    ----------
+    triangulation: scipy.spatial.Delaunay object
+        Delaunay triangulation of an input grid, shapes of original grid has to match
+        target's one.
+    target: np.ndarray, shape=(2)
+        Target point, of which the nearest facet on the convex hull is wanted.
+        Has to lay outside of the grid for a usable result, which is not checked here.
+
+    Returns
+    -------
+    simplex_ind: int
+        Index of the nearest simplex to a target point.
+    """
+    qhull = triangulation.convex_hull
+    qhull_points = triangulation.points[qhull]
+
+    nearest_facet_ind = find_nearest_facet(qhull_points, target)
+    facet_points = qhull_points[nearest_facet_ind]
+
+    simplices_points = triangulation.points[triangulation.simplices]
+
+    return find_simplex_to_facet(simplices_points, facet_points)


Test patch:
diff --git a/pyirf/interpolation/tests/test_utils.py b/pyirf/interpolation/tests/test_utils.py
new file mode 100644
index 0000000..60fe5f7
--- /dev/null
+++ b/pyirf/interpolation/tests/test_utils.py
@@ -0,0 +1,149 @@
+import numpy as np
+import pytest
+from scipy.spatial import Delaunay
+
+
+@pytest.fixture
+def non_rect_grid():
+    grid = np.array([[0, 0], [10, 20], [30, 20], [20, 0], [40, 0]])
+
+    return Delaunay(grid)
+
+
+def test_plumb_point_distance():
+    """Test line-segment to point distance computation"""
+    from pyirf.interpolation.utils import plumb_point_dist
+
+    # Test vertical line
+    line = np.array([[0, 0], [0, 1]])
+
+    # Plumb point between end-points
+    assert plumb_point_dist(line, np.array([-1, 0.5])) == 1
+    assert plumb_point_dist(line, np.array([-0.7, 0.25])) == 0.7
+
+    # Plumb point on one end-point
+    assert plumb_point_dist(line, np.array([0, 2.1])) == 1.1
+    assert plumb_point_dist(line, np.array([0, -1])) == 1
+
+    # Plumb point not between end-points, nearest point is (0, 0)
+    assert plumb_point_dist(line, np.array([-1, -1])) == np.sqrt(2)
+    # Nearest point is (0, 1)
+    assert plumb_point_dist(line, np.array([3, 3])) == np.sqrt(13)
+
+    # Test horzontal line
+    line = np.array([[0, 0], [1, 0]])
+
+    # Plumb point between end-points
+    assert plumb_point_dist(line, np.array([0.5, 0.5])) == 0.5
+
+    # Plumb point in extention of line, nearest point is (1, 0)
+    assert plumb_point_dist(line, np.array([2, 0])) == 1
+
+    # Plumb point on end point
+    assert plumb_point_dist(line, np.array([1, 1])) == 1
+
+    # Nearest point is (0, 0)
+    assert plumb_point_dist(line, np.array([-1, -1])) == np.sqrt(2)
+
+    # Test arbitrary line
+    line = np.array([[1, 1], [-1, -1]])
+    # isclose needed here, as there is a small numerical deviation
+    # of +/- eps in this case. Plumb point between end-points
+    assert np.isclose(plumb_point_dist(line, np.array([-1, 1])), np.sqrt(2))
+
+    # Nearest point is (-1, -1)
+    assert plumb_point_dist(line, np.array([-2, -3])) == np.sqrt(5)
+
+
+def test_point_facet_angle():
+    """Test angle computation in triangle, function should return cos(angle)"""
+    from pyirf.interpolation.utils import point_facet_angle
+
+    line = np.array([[0, 0], [0, 1]])
+
+    assert np.isclose(
+        point_facet_angle(line, np.array([1, 0])), np.cos(45 * np.pi / 180)
+    )
+    assert np.isclose(
+        point_facet_angle(line, np.array([-1, 0])), np.cos(45 * np.pi / 180)
+    )
+    # these points build a same-side triangle
+    assert np.isclose(
+        point_facet_angle(line, np.array([np.sqrt(3) / 2, 0.5])),
+        np.cos(60 * np.pi / 180),
+    )
+
+
+def test_find_nearest_facet_rect_grid():
+    """Test nearest facet finding on rectanguar grid"""
+    from pyirf.interpolation.utils import find_nearest_facet
+
+    rect_grid = Delaunay(
+        np.array([[0, 0], [0, 20], [20, 20], [20, 0], [40, 0], [40, 20]])
+    )
+    qhull_points = rect_grid.points[rect_grid.convex_hull]
+
+    nearest_facet_ind = find_nearest_facet(qhull_points, np.array([10, -5]))
+    assert np.logical_or(
+        np.array_equal(qhull_points[nearest_facet_ind], np.array([[0, 0], [20, 0]])),
+        np.array_equal(qhull_points[nearest_facet_ind], np.array([[20, 0], [0, 0]])),
+    )
+
+    nearest_facet_ind = find_nearest_facet(qhull_points, np.array([45, 15]))
+    assert np.logical_or(
+        np.array_equal(qhull_points[nearest_facet_ind], np.array([[40, 0], [40, 20]])),
+        np.array_equal(qhull_points[nearest_facet_ind], np.array([[40, 20], [40, 0]])),
+    )
+
+    nearest_facet_ind = find_nearest_facet(qhull_points, np.array([-10, -1]))
+    assert np.logical_or(
+        np.array_equal(qhull_points[nearest_facet_ind], np.array([[0, 0], [0, 20]])),
+        np.array_equal(qhull_points[nearest_facet_ind], np.array([[0, 20], [0, 0]])),
+    )
+
+
+def test_find_nearest_facet_non_rect_grid(non_rect_grid):
+    """Test nearest facet finding on a non rectanguar grid to catch some more cases"""
+    from pyirf.interpolation.utils import find_nearest_facet
+
+    qhull_points = non_rect_grid.points[non_rect_grid.convex_hull]
+
+    nearest_facet_ind = find_nearest_facet(qhull_points, np.array([5, 20]))
+    assert np.logical_or(
+        np.array_equal(qhull_points[nearest_facet_ind], np.array([[0, 0], [10, 20]])),
+        np.array_equal(qhull_points[nearest_facet_ind], np.array([[10, 20], [0, 0]])),
+    )
+
+    nearest_facet_ind = find_nearest_facet(qhull_points, np.array([35, 30]))
+    assert np.logical_or(
+        np.array_equal(qhull_points[nearest_facet_ind], np.array([[10, 20], [30, 20]])),
+        np.array_equal(qhull_points[nearest_facet_ind], np.array([[30, 20], [10, 20]])),
+    )
+
+
+def test_find_simplex_to_facet(non_rect_grid):
+    """
+    Test facet-to-simplex finding on non rectangular grid, as the triangulation
+    is clear in this case as it is build from left to right and not ambiguous.
+    For the rectangular grid used above two triangulations exist.
+    """
+    from pyirf.interpolation.utils import find_simplex_to_facet
+
+    simplices_points = non_rect_grid.points[non_rect_grid.simplices]
+
+    assert find_simplex_to_facet(simplices_points, np.array([[0, 0], [0, 20]])) == 0
+    assert find_simplex_to_facet(simplices_points, np.array([[10, 20], [30, 20]])) == 1
+    assert find_simplex_to_facet(simplices_points, np.array([[30, 20], [40, 0]])) == 2
+
+
+def test_find_nearest_simplex(non_rect_grid):
+    """
+    Test whole nearest simplex finding on non rectangular grid, as the triangulation
+    is clear in this case as it is build from left to right and not ambiguous.
+    For the rectangular grid used above two triangulations exist.
+    """
+    from pyirf.interpolation.utils import find_nearest_simplex
+
+    assert find_nearest_simplex(non_rect_grid, np.array([-10, -10])) == 0
+    assert find_nearest_simplex(non_rect_grid, np.array([10, 30])) == 1
+    assert find_nearest_simplex(non_rect_grid, np.array([20.00000000001, -10])) == 2