WIP: (feature) add anisotropic meshing to LearnerND

Closes #74

Depends on !86 (merged) and #80 and therefore it has the corresponding branches included as well

Still has a few to-do's:

• let LearnerND.ip() make use of our triangulation rather than building a new one
• make it work in arbitrary dimensions
• verify that it is beneficial
• let the user configure the parameters (maximum stretch factor and number of points to take into account) Use one simplex and it's neighbours
• make test pass
  • add rtree as install requirement No more RTree anymore :)
  • raise exception if anistropic=True and rtree not installed, pass if anisotropic=False
• refactor code to be human-readable
• let's make it fast :) #80

As it seems it doesn't work that well with the ring, since this ring has a relative low average gradient and very high second derivative. So maybe this second derivative might be a more useful property to determine the

Sneak peek:

adaptive/learner/learnerND.py

@@ -12,8 +12,9 @@ from .base_learner import BaseLearner
 
 from ..notebook_integration import ensure_holoviews
 from .triangulation import (Triangulation, point_in_simplex,
-                            circumsphere, simplex_volume_in_embedding)
+                            circumsphere, simplex_volume_in_embedding, FakeDelaunay)
 from ..utils import restore
+import math
 
 
 def volume(simplex, ys=None):
@@ -155,7 +156,7 @@ class LearnerND(BaseLearner):
     children based on volume.
     """
 
-    def __init__(self, func, bounds, loss_per_simplex=None):
+    def __init__(self, func, bounds, loss_per_simplex=None, anisotropic=False):
         self.ndim = len(bounds)
         self._vdim = None
         self.loss_per_simplex = loss_per_simplex or default_loss
@@ -175,10 +176,15 @@ class LearnerND(BaseLearner):
         self._subtriangulations = dict()  # simplex -> triangulation
 
         # scale to unit
-        self._transform = np.linalg.inv(np.diag(np.diff(bounds).flat))
+        self._scale_matrix = np.linalg.inv(np.diag(np.diff(bounds).flat))
 
         # create a private random number generator with fixed seed
         self._random = random.Random(1)
+        self.anisotripic = anisotropic
+        
+        # all real triangles that have not been subdivided and the pending
+        # triangles heap of tuples (-loss, real simplex, sub_simplex or None)
+        self._losses_combined = [] # heap
 
         # all real triangles that have not been subdivided and the pending 
         # triangles heap of tuples (-loss, real simplex, sub_simplex or None)
@@ -212,8 +218,9 @@ class LearnerND(BaseLearner):
         return all(p in self.data for p in self._bounds_points)
 
     def ip(self):
-        # XXX: take our own triangulation into account when generating the ip
-        return interpolate.LinearNDInterpolator(self.points, self.values)
+        tri = FakeDelaunay(self.tri.vertices, self.tri.simplices)
+        values = [self.data[k] for k in self.tri.vertices]
+        return interpolate.LinearNDInterpolator(tri, values)
 
     @property
     def tri(self):
@@ -240,7 +247,6 @@ class LearnerND(BaseLearner):
 
     def tell(self, point, value):
         point = tuple(point)
-
         if point in self.data:
             return  # we already know about the point
 
@@ -257,17 +263,60 @@ class LearnerND(BaseLearner):
         if tri is not None:
             simplex = self._pending_to_simplex.get(point)
             if simplex is not None and not self._simplex_exists(simplex):
-                simplex = None
-            to_delete, to_add = tri.add_point(
-                point, simplex, transform=self._transform)
+                simplex = self.tri.locate_point(point)
+            transform = self.get_local_transform_matrix(simplex)
+            to_delete, to_add = self._tri.add_point(
+                point, simplex, transform=transform)
             self.update_losses(to_delete, to_add)
 
+    def get_local_transform_matrix(self, simplex):
+        scale = self._scale_matrix
+
+        if simplex is None or self.tri is None or not self.anisotripic:
+            return scale
+        
+        neighbours = set.union(*[self.tri.vertex_to_simplices[i] for i in simplex])
+        indices = set.union(set(), *neighbours)
+
+        points = np.array([self.points[i] for i in indices])
+        values = [self.data[tuple(p)] for p in points]
+
+        if isinstance(values[0], Iterable):
+            raise ValueError("Anisotropic learner only works with 1D output")
+
+        # Do a linear least square fit
+        # A x = B, find x
+        points = np.dot(scale, points.T).T
+        ones = np.ones((len(points), 1))
+        A = np.hstack((points, ones))
+        B = np.array(values, dtype=float)
+        fit = np.linalg.lstsq(A, B, rcond=None)[0]
+        *gradient, _constant = fit
+        # we do not need the constant, only the gradient
+
+        # gradient is a vector of the amount of the slope in each direction
+        magnitude = np.linalg.norm(gradient)
+
+        if np.isclose(magnitude, 0):
+            return scale  # there is no noticable gradient
+
+        # Make it a 2d numpy array and normalise it.
+        gradient = np.array([gradient], dtype=float) / magnitude
+        projection_matrix = (gradient.T @ gradient)
+        identity = np.eye(self.ndim)
+
+        factor = math.sqrt(magnitude ** 2 + 1) - 1
+
+        scale_along_gradient = projection_matrix * factor + identity
+        m = np.dot(scale_along_gradient, scale)
+        return m.T
+
     def _simplex_exists(self, simplex):
         simplex = tuple(sorted(simplex))
         return simplex in self.tri.simplices
 
-    def inside_bounds(self, point):
-        return all(mn <= p <= mx for p, (mn, mx) in zip(point, self.bounds))
+    def inside_bounds(self, point, eps = 1e-8):
+        return all((mn - eps) <= p <= (mx + eps) for p, (mn, mx) in zip(point, self.bounds))
             
     def tell_pending(self, point, *, simplex=None):
         point = tuple(point)
@@ -384,8 +433,8 @@ class LearnerND(BaseLearner):
             subtri = self._subtriangulations[simplex]
             points = subtri.get_vertices(subsimplex)
 
-        point_new = tuple(choose_point_in_simplex(points,
-                                                  transform=self._transform))
+        transform = self.get_local_transform_matrix(simplex)
+        point_new = tuple(choose_point_in_simplex(points, transform=transform))
 
         self._pending_to_simplex[point_new] = simplex
         self.tell_pending(point_new, simplex=simplex)  # O(??)
@@ -543,7 +592,7 @@ class LearnerND(BaseLearner):
             p = hv.Path((x, y))
 
             # Plot with 5% margins such that the boundary points are visible
-            margin = 0.05 / self._transform[ind, ind]
+            margin = 0.05 / self._scale_matrix[ind, ind]
             plot_bounds = (x.min() - margin, x.max() + margin)
             return p.redim(x=dict(range=plot_bounds))