refactor: clarify 1D triangulation and interpolation

Author: basnijholt

adaptive/learner/learnerND.py

@@ -493,21 +493,39 @@ def load_dataframe(  # type: ignore[override]
     def bounds_are_done(self):
         return all(p in self.data for p in self._bounds_points)
 
+    def _sorted_line_data(self):
+        coordinates = self.points[:, 0]
+        sorted_indices = np.argsort(coordinates)
+        return coordinates[sorted_indices], self.values[sorted_indices]
+
     def _ip(self):
-        """A `scipy.interpolate.LinearNDInterpolator` instance
-        containing the learner's data."""
+        """A SciPy interpolator containing the learner's data."""
         # XXX: take our own triangulation into account when generating the _ip
-        if self.ndim == 1:
-            points = self.points.ravel()
-            sorted_idx = np.argsort(points)
-            return interpolate.interp1d(
-                points[sorted_idx],
-                self.values[sorted_idx],
-                axis=0,
-                bounds_error=False,
-                fill_value=np.nan,
-            )
-        return interpolate.LinearNDInterpolator(self.points, self.values)
+        if self.ndim != 1:
+            return interpolate.LinearNDInterpolator(self.points, self.values)
+
+        coordinates, values = self._sorted_line_data()
+        return interpolate.interp1d(
+            coordinates,
+            values,
+            axis=0,
+            bounds_error=False,
+            fill_value=np.nan,
+        )
+
+    def _plot_1d(self, n=None):
+        hv = ensure_holoviews()
+        if len(self.data) < 2:
+            return hv.Path([]) * hv.Scatter([]).opts(size=5)
+
+        (x_bounds,) = self._bbox
+        n = n or 201
+        xs = np.linspace(*x_bounds, n)
+        ys = self._ip()(xs)
+        scatter_points = [
+            (point[0], value) for point, value in sorted(self.data.items())
+        ]
+        return hv.Path((xs, ys)) * hv.Scatter(scatter_points).opts(size=5)
 
     @property
     def tri(self):
@@ -891,34 +909,23 @@ def remove_unfinished(self):
     ##########################
 
     def plot(self, n=None, tri_alpha=0):
-        """Plot the function we want to learn, only works in 2D.
+        """Plot the function we want to learn in 1D or 2D.
 
         Parameters
         ----------
         n : int
-            the number of boxes in the interpolation grid along each axis
+            The number of interpolation points per axis.
         tri_alpha : float (0 to 1)
             Opacity of triangulation lines
         """
-        hv = ensure_holoviews()
         if self.vdim > 1:
             raise NotImplementedError(
                 "holoviews currently does not support", "3D surface plots in bokeh."
             )
         if self.ndim == 1:
-            if len(self.data) >= 2:
-                (x,) = self._bbox
-                n = n or 201
-                xs = np.linspace(x[0], x[1], n)
-                ys = self._ip()(xs)
-                path = hv.Path((xs, ys))
-                points = [(x[0], y) for x, y in sorted(self.data.items())]
-                scatter = hv.Scatter(points)
-            else:
-                path = hv.Path([])
-                scatter = hv.Scatter([])
-            return path * scatter.opts(size=5)
+            return self._plot_1d(n)
 
+        hv = ensure_holoviews()
         if self.ndim != 2:
             raise NotImplementedError(
                 "Only 1D and 2D plots are implemented: You can "

adaptive/learner/triangulation.py

@@ -231,6 +231,14 @@ def is_iterable_and_sized(obj):
     return isinstance(obj, Iterable) and isinstance(obj, Sized)
 
 
+def _flat_simplices(coords):
+    """Yield the intervals of a 1D triangulation in sorted coordinate order."""
+    sorted_indices = sorted(range(len(coords)), key=coords.__getitem__)
+    for left, right in zip(sorted_indices, sorted_indices[1:]):
+        if coords[left] != coords[right]:
+            yield left, right
+
+
 def simplex_volume_in_embedding(vertices) -> float:
     """Calculate the volume of a simplex in a higher dimensional embedding.
     That is: dim > len(vertices) - 1. For example if you would like to know the
@@ -257,12 +265,13 @@ def simplex_volume_in_embedding(vertices) -> float:
     # Modified from https://codereview.stackexchange.com/questions/77593/calculating-the-volume-of-a-tetrahedron
 
     vertices = asarray(vertices, dtype=float)
-    if len(vertices) == 2:
-        # 1-simplex (line segment): volume is the Euclidean distance
+    num_vertices = len(vertices)
+    if num_vertices == 2:
+        # A 1-simplex is just a line segment.
         return float(norm(vertices[1] - vertices[0]))
 
-    dim = len(vertices[0])
-    if dim == 2:
+    embedding_dim = len(vertices[0])
+    if embedding_dim == 2:
         # Heron's formula
         a, b, c = scipy.spatial.distance.pdist(vertices, metric="euclidean")
         s = 0.5 * (a + b + c)
@@ -272,13 +281,12 @@ def simplex_volume_in_embedding(vertices) -> float:
     sq_dists = scipy.spatial.distance.pdist(vertices, metric="sqeuclidean")
 
     # Add border while compressed
-    num_verts = scipy.spatial.distance.num_obs_y(sq_dists)
-    bordered = concatenate((ones(num_verts), sq_dists))
+    bordered = concatenate((ones(num_vertices), sq_dists))
 
     # Make matrix and find volume
     sq_dists_mat = scipy.spatial.distance.squareform(bordered)
 
-    coeff = -((-2) ** (num_verts - 1)) * factorial(num_verts - 1) ** 2
+    coeff = -((-2) ** (num_vertices - 1)) * factorial(num_vertices - 1) ** 2
     vol_square = fast_det(sq_dists_mat) / coeff
 
     if vol_square < 0:
@@ -346,19 +354,14 @@ def __init__(self, coords):
         self.vertex_to_simplices = [set() for _ in coords]
 
         if dim == 1:
-            # For 1D, sort points and create intervals as simplices,
-            # skipping adjacent duplicates to avoid degenerate zero-volume simplices
-            sorted_indices = sorted(range(len(coords)), key=lambda i: coords[i])
-            for i in range(len(sorted_indices) - 1):
-                if coords[sorted_indices[i]] == coords[sorted_indices[i + 1]]:
-                    continue
-                self.add_simplex((sorted_indices[i], sorted_indices[i + 1]))
+            simplices = _flat_simplices(coords)
         else:
-            # find a Delaunay triangulation to start with, then we will throw it
-            # away and continue with our own algorithm
-            initial_tri = scipy.spatial.Delaunay(coords)
-            for simplex in initial_tri.simplices:
-                self.add_simplex(simplex)
+            # Find a Delaunay triangulation to start with, then we will throw it
+            # away and continue with our own algorithm.
+            simplices = scipy.spatial.Delaunay(coords).simplices
+
+        for simplex in simplices:
+            self.add_simplex(simplex)
 
     def delete_simplex(self, simplex):
         simplex = tuple(sorted(simplex))