Resolve "(LearnerND) add iso-surface plot feature"

Closes #112 (closed)

adaptive/learner/learnerND.py

@@ -463,6 +463,10 @@ class LearnerND(BaseLearner):
         self._subtriangulations = dict()
         self._pending_to_simplex = dict()
 
+    ##########################
+    # Plotting related stuff #
+    ##########################
+
     def plot(self, n=None, tri_alpha=0):
         """Plot the function we want to learn, only works in 2D.
 
@@ -639,6 +643,103 @@ class LearnerND(BaseLearner):
             )
 
         return vertices, faces
+    
+    def _get_isoline(self, level=0.0):
+        # Very similar to _get_isosurface, maybe merge the two functions
+        if self.ndim != 2 or self.vdim != 1:
+            raise Exception('Isoline plotting is only supported'
+                            ' for a 2D input and 1D output')
+
+        vertices = []  # index -> (x,y,z)
+        lines = []  # tuple of indices of the corner points
+        from_line_to_vertex = {}  # the interpolated vertex (index) between two known points
+
+        def _get_vertex_index(a, b):
+            if (a, b) in from_line_to_vertex:
+                return from_line_to_vertex[(a, b)]
+
+            # Otherwise compute it and cache the result.
+            vertex_a = self.tri.vertices[a]
+            vertex_b = self.tri.vertices[b]
+            value_a = self.data[vertex_a]
+            value_b = self.data[vertex_b]
+
+            da = abs(value_a - level)
+            db = abs(value_b - level)
+            dab = da + db
+
+            new_pt = (db / dab * np.array(vertex_a)
+                      + da / dab * np.array(vertex_b))
+
+            new_index = len(vertices)
+            vertices.append(new_pt)
+            from_line_to_vertex[(a, b)] = new_index
+            return new_index
+
+        for simplex in self.tri.simplices:
+            line = []
+            for a, b in itertools.combinations(simplex, 2):
+                va = self.data[self.tri.vertices[a]]
+                vb = self.data[self.tri.vertices[b]]
+                if min(va, vb) < level <= max(va, vb):
+                    vi = _get_vertex_index(a, b)
+                    should_add = True
+                    for pi in line:
+                        if np.allclose(vertices[vi], vertices[pi]):
+                            should_add = False
+                    if should_add:
+                        line.append(vi)
+
+            if len(line) == 2:
+                lines.append(line)
+
+        if len(lines) == 0:
+            r_min = min(self.data[v] for v in self.tri.vertices)
+            r_max = max(self.data[v] for v in self.tri.vertices)
+
+            raise ValueError(
+                f"Could not draw isosurface for level={level}, as"
+                " this value is not inside the function range. Please choose"
+                f" a level strictly inside interval ({r_min}, {r_max})"
+            )
+
+        return vertices, lines
+
+    def plot_isoline(self, level=0.0, n=None):
+        """Plot the isoline at a specific level of the function we want to
+        learn, only works in 2D.
+
+        Parameters
+        ----------
+        level : float
+            the value of the function at which you would like to see the isoline
+        n : int
+            the number of boxes in the interpolation grid along each axis
+        """
+        hv = ensure_holoviews()
+
+        vertices, lines = self._get_isoline(level)
+        paths = [[vertices[i], vertices[j]] for i,j in lines]
+        contour = hv.Path(paths)
+
+        x, y = self.bounds
+        lbrt = x[0], y[0], x[1], y[1]
+
+        if n is None:
+            # Calculate how many grid points are needed.
+            # factor from A=√3/4 * a² (equilateral triangle)
+            n = int(0.658 / np.sqrt(np.min(self.tri.volumes())))
+
+        xs = ys = np.linspace(0, 1, n)
+        xs = xs * (x[1] - x[0]) + x[0]
+        ys = ys * (y[1] - y[0]) + y[0]
+        z = self.ip()(xs[:, None], ys[None, :]).squeeze()
+
+        im = hv.Image(np.rot90(z), bounds=lbrt)
+
+        im_opts = dict(cmap='viridis')
+        contour_opts = dict(color='black')
+        return im.opts(style=im_opts) * contour.opts(style=contour_opts)
 
     def plot_isosurface(self, level=0.0, hull_opacity=0.2):
         """Plots the linearly interpolated isosurface of the function,