added a curvature_loss function to learner1D

Author: jhoofwijk

adaptive/learner/learner1D.py

@@ -3,11 +3,14 @@
 import heapq
 import itertools
 import math
+from collections import Iterable
 
 import numpy as np
 import sortedcontainers
 
 from .base_learner import BaseLearner
+from .learnerND import volume
+from .triangulation import simplex_volume_in_embedding
 from ..notebook_integration import ensure_holoviews
 from ..utils import cache_latest
 
@@ -56,6 +59,45 @@ def default_loss(interval, scale, function_values):
     return loss
 
 
+def _loss_of_multi_interval(xs, ys):
+    N = len(xs) - 2
+    if isinstance(ys[0], Iterable):
+        pts = [(x, *y) for x, y in zip(xs, ys)]
+        vol = simplex_volume_in_embedding
+    else:
+        pts = [(x, y) for x, y in zip(xs, ys)]
+        vol = volume
+    return sum(vol(pts[i:i+3]) for i in range(N)) / N
+
+
+def triangle_loss(interval, neighbours, scale, function_values):
+    x_left, x_right = interval
+    neighbour_left, neighbour_right = neighbours
+    xs = [neighbour_left, x_left, x_right, neighbour_right]
+    # The neighbours could be None if we are at the boundary, in that case we
+    # have to filter this out
+    xs = [x for x in xs if x is not None]
+
+    if len(xs) <= 2:
+        return (x_right - x_left) / scale[0]
+    else:
+        y_scale = scale[1] or 1
+        ys_scaled = [function_values[x] / y_scale for x in xs]
+        xs_scaled = [x / scale[0] for x in xs]
+        return _loss_of_multi_interval(xs_scaled, ys_scaled)
+
+
+def get_curvature_loss(area_factor=1, euclid_factor=0.02, horizontal_factor=0.02):
+    def curvature_loss(interval, neighbours, scale, function_values):
+        triangle_loss_ = triangle_loss(interval, neighbours, scale, function_values)
+        default_loss_ = default_loss(interval, scale, function_values)
+        dx = (interval[1] - interval[0]) / scale[0]
+        return (area_factor * (triangle_loss_**0.5)
+                + euclid_factor * default_loss_
+                + horizontal_factor * dx)
+    return curvature_loss
+
+
 def linspace(x_left, x_right, n):
     """This is equivalent to
     'np.linspace(x_left, x_right, n, endpoint=False)[1:]',
@@ -116,9 +158,14 @@ class Learner1D(BaseLearner):
         to have values for both of the points in 'interval'.
     """
 
-    def __init__(self, function, bounds, loss_per_interval=None):
+    def __init__(self, function, bounds, loss_per_interval=None, loss_depends_on_neighbours=False):
         self.function = function
-        self.loss_per_interval = loss_per_interval or default_loss
+        self._loss_depends_on_neighbours = loss_depends_on_neighbours
+
+        if loss_depends_on_neighbours:
+            self.loss_per_interval = loss_per_interval or get_curvature_loss()
+        else:
+            self.loss_per_interval = loss_per_interval or default_loss
 
         # A dict storing the loss function for each interval x_n.
         self.losses = {}
@@ -176,25 +223,42 @@ def loss(self, real=True):
         losses = self.losses if real else self.losses_combined
         return max(losses.values()) if len(losses) > 0 else float('inf')
 
+    def _get_loss_in_interval(self, x_left, x_right):
+        assert x_left is not None and x_right is not None
+
+        if x_right - x_left < self._dx_eps:
+            return 0
+
+        # we need to compute the loss for this interval
+        interval = (x_left, x_right)
+        if self._loss_depends_on_neighbours:
+            neighbour_left = self.neighbors.get(x_left, (None, None))[0]
+            neighbour_right = self.neighbors.get(x_right, (None, None))[1]
+            neighbours = neighbour_left, neighbour_right
+            return self.loss_per_interval(interval, neighbours,
+                                          self._scale, self.data)
+        else:
+            return self.loss_per_interval(interval, self._scale, self.data)
+
+
     def _update_interpolated_loss_in_interval(self, x_left, x_right):
-        if x_left is not None and x_right is not None:
-            dx = x_right - x_left
-            if dx < self._dx_eps:
-                loss = 0
-            else:
-                loss = self.loss_per_interval((x_left, x_right),
-                                              self._scale, self.data)
-            self.losses[x_left, x_right] = loss
-
-            # Iterate over all interpolated intervals in between
-            # x_left and x_right and set the newly interpolated loss.
-            a, b = x_left, None
-            while b != x_right:
-                b = self.neighbors_combined[a][1]
-                self.losses_combined[a, b] = (b - a) * loss / dx
-                a = b
+        if x_left is None or x_right is None:
+            return
+
+        loss = self._get_loss_in_interval(x_left, x_right)
+        self.losses[x_left, x_right] = loss
+
+        # Iterate over all interpolated intervals in between
+        # x_left and x_right and set the newly interpolated loss.
+        a, b = x_left, None
+        dx = x_right - x_left
+        while b != x_right:
+            b = self.neighbors_combined[a][1]
+            self.losses_combined[a, b] = (b - a) * loss / dx
+            a = b
 
     def _update_losses(self, x, real=True):
+        """Update all losses that depend on x"""
         # When we add a new point x, we should update the losses
         # (x_left, x_right) are the "real" neighbors of 'x'.
         x_left, x_right = self._find_neighbors(x, self.neighbors)
@@ -212,6 +276,13 @@ def _update_losses(self, x, real=True):
             self._update_interpolated_loss_in_interval(x_left, x)
             self._update_interpolated_loss_in_interval(x, x_right)
 
+            # if the loss depends on the neighbors we should also update those losses
+            if self._loss_depends_on_neighbours:
+                neighbour_left = self.neighbors.get(x_left, (None, None))[0]
+                neighbour_right = self.neighbors.get(x_right, (None, None))[1]
+                self._update_interpolated_loss_in_interval(neighbour_left, x_left)
+                self._update_interpolated_loss_in_interval(x_right, neighbour_right)
+
             # Since 'x' is in between (x_left, x_right),
             # we get rid of the interval.
             self.losses.pop((x_left, x_right), None)
@@ -358,7 +429,7 @@ def tell_many(self, xs, ys, *, force=False):
         self.losses = {}
         for x_left, x_right in intervals:
             self.losses[x_left, x_right] = (
-                self.loss_per_interval((x_left, x_right), self._scale, self.data)
+                self._get_loss_in_interval(x_left, x_right)
                 if x_right - x_left >= self._dx_eps else 0)
 
         # List with "real" intervals that have interpolated intervals inside

adaptive/learner/triangulation.py

@@ -229,7 +229,9 @@ def simplex_volume_in_embedding(vertices) -> float:
     coeff = - (-2) ** (num_verts-1) * factorial(num_verts-1) ** 2
     vol_square = np.linalg.det(sq_dists_mat) / coeff
 
-    if vol_square <= 0:
+    if vol_square < 0:
+        if abs(vol_square) < 1e-15:
+            return 0
         raise ValueError('Provided vertices do not form a simplex')
 
     return np.sqrt(vol_square)

adaptive/tests/test_learner1d.py

@@ -4,6 +4,7 @@
 import numpy as np
 
 from ..learner import Learner1D
+from ..learner.learner1D import get_curvature_loss
 from ..runner import simple
 
 
@@ -120,9 +121,9 @@ def test_termination_on_discontinuities():
     smallest_interval = min(abs(a - b) for a, b in learner.losses.keys())
     assert smallest_interval >= np.finfo(float).eps
 
-    learner = _run_on_discontinuity(0.5E3, (-1E3, 1E3))
+    learner = _run_on_discontinuity(0.5e3, (-1e3, 1e3))
     smallest_interval = min(abs(a - b) for a, b in learner.losses.keys())
-    assert smallest_interval >= 0.5E3 * np.finfo(float).eps
+    assert smallest_interval >= 0.5e3 * np.finfo(float).eps
 
 
 def test_order_adding_points():
@@ -340,3 +341,21 @@ def _random_run(learner, learner2, scale_doubling=True):
     learner2 = Learner1D(f, bounds=(-1, 1))
     _random_run(learner, learner2, scale_doubling=True)
     test_equal(learner, learner2)
+
+
+def test_curvature_loss():
+    def f(x):
+        return np.tanh(20*x)
+
+    learner = Learner1D(f, (-1, 1), loss_per_interval=get_curvature_loss(), loss_depends_on_neighbours=True)
+    simple(learner, goal=lambda l: l.npoints > 100)
+    # assert this is reached without error
+
+
+def test_curvature_loss_vectors():
+    def f(x):
+        return np.tanh(20*x), np.tanh(20*(x-0.4))
+
+    learner = Learner1D(f, (-1, 1), loss_per_interval=get_curvature_loss(), loss_depends_on_neighbours=True)
+    simple(learner, goal=lambda l: l.npoints > 100)
+    assert learner.npoints > 100