From 4d03a752b46c1dc9925d5a0de2b6251fd2eedbdc Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Tue, 17 Jul 2018 20:52:23 +0200
Subject: [PATCH 01/32] make loss function based on simplex volume in higher
 dimensional space taking into account the output dimension of the function

---
 adaptive/learner/learnerND.py     | 48 +++++++++++++++++++++++++++++--
 adaptive/learner/triangulation.py |  2 +-
 2 files changed, 47 insertions(+), 3 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index ce6ec128d..9409d4b58 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -1,5 +1,5 @@
 # -*- coding: utf-8 -*-
-from collections import OrderedDict
+from collections import OrderedDict, Iterable
 import itertools
 import heapq
 
@@ -12,6 +12,7 @@
 
 from .triangulation import Triangulation
 import random
+from math import factorial
 
 
 def find_initial_simplex(pts, ndim):
@@ -56,8 +57,51 @@ def std_loss(simplex, ys):
     return r.flat * np.power(vol, 1. / dim) + vol
 
 
+def simplex_volume(vertices) -> float:
+    """
+    Return the volume of the simplex with given vertices.
+
+    If vertices are given they must be in a arraylike with shape (N+1, N):
+    the position vectors of the N+1 vertices in N dimensions. If the sides
+    are given, they must be the compressed pairwise distance matrix as
+    returned from scipy.spatial.distance.pdist.
+
+    Raises a ValueError if the vertices do not form a simplex (for example,
+    because they are coplanar, colinear or coincident).
+
+    Warning: this algorithm has not been tested for numerical stability.
+    """
+
+    # Implements http://mathworld.wolfram.com/Cayley-MengerDeterminant.html
+    # Modified from https://codereview.stackexchange.com/questions/77593/calculating-the-volume-of-a-tetrahedron
+
+    # β_ij = |v_i - v_k|²
+    vertices = np.array(vertices, dtype=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 = np.concatenate((np.ones(num_verts), 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
+    vol_square = np.linalg.det(sq_dists_mat) / coeff
+
+    if vol_square <= 0:
+        raise ValueError('Provided vertices do not form a tetrahedron')
+
+    return np.sqrt(vol_square)
+
+
 def default_loss(simplex, ys):
-    return std_loss(simplex, ys)
+    # return std_loss(simplex, ys)
+    if isinstance(ys[0], Iterable):
+        pts = [(*s, *ys[i]) for i, s in enumerate(simplex)]
+    else:
+        pts = [(*s, ys[i]) for i, s in enumerate(simplex)]
+    return simplex_volume(pts)
 
 
 def choose_point_in_simplex(simplex, transform=None):
diff --git a/adaptive/learner/triangulation.py b/adaptive/learner/triangulation.py
index b1c669231..0a24f29d4 100644
--- a/adaptive/learner/triangulation.py
+++ b/adaptive/learner/triangulation.py
@@ -1,4 +1,4 @@
-from collections import defaultdict, Counter, Sized, Iterable
+from collections import Counter, Sized, Iterable
 from itertools import combinations, chain
 
 import numpy as np

From d9a45e1bcc51508d5879d65a0ca97717f4730275 Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Tue, 17 Jul 2018 21:24:43 +0200
Subject: [PATCH 02/32] make point_in_simplex and circumcircle global methods

this way you can import them into other files
---
 adaptive/learner/triangulation.py | 71 +++++++++++++++++--------------
 1 file changed, 40 insertions(+), 31 deletions(-)

diff --git a/adaptive/learner/triangulation.py b/adaptive/learner/triangulation.py
index 0a24f29d4..df7df1cf5 100644
--- a/adaptive/learner/triangulation.py
+++ b/adaptive/learner/triangulation.py
@@ -29,6 +29,17 @@ def fast_2d_point_in_simplex(point, simplex, eps=1e-8):
     return (t >= -eps) and (s + t <= 1 + eps)
 
 
+def point_in_simplex(point, simplex, eps=1e-8):
+    if len(point) == 2:
+        return fast_2d_point_in_simplex(point, simplex, eps)
+
+    x0 = np.array(simplex[0], dtype=float)
+    vectors = np.array(simplex[1:], dtype=float) - x0
+    alpha = np.linalg.solve(vectors.T, point - x0)
+
+    return all(alpha > -eps) and sum(alpha) < 1 + eps
+
+
 def fast_2d_circumcircle(points):
     """Compute the center and radius of the circumscribed circle of a triangle
 
@@ -106,6 +117,32 @@ def fast_3d_circumcircle(points):
     return tuple(center), radius
 
 
+def circumsphere(pts):
+    dim = len(pts) - 1
+    if dim == 2:
+        return fast_2d_circumcircle(pts)
+    if dim == 3:
+        return fast_3d_circumcircle(pts)
+
+    # Modified method from http://mathworld.wolfram.com/Circumsphere.html
+    mat = [[np.sum(np.square(pt)), *pt, 1] for pt in pts]
+
+    center = []
+    for i in range(1, len(pts)):
+        r = np.delete(mat, i, 1)
+        factor = (-1) ** (i + 1)
+        center.append(factor * np.linalg.det(r))
+
+    a = np.linalg.det(np.delete(mat, 0, 1))
+    center = [x / (2 * a) for x in center]
+
+    x0 = pts[0]
+    vec = np.subtract(center, x0)
+    radius = fast_norm(vec)
+
+    return tuple(center), radius
+
+
 def orientation(face, origin):
     """Compute the orientation of the face with respect to a point, origin.
 
@@ -227,15 +264,8 @@ def get_reduced_simplex(self, point, simplex, eps=1e-8) -> list:
         return [simplex[i] for i in result]
 
     def point_in_simplex(self, point, simplex, eps=1e-8):
-        if self.dim == 2:
-            vertices = self.get_vertices(simplex)
-            return fast_2d_point_in_simplex(point, vertices, eps)
-        elif self.dim == 3:
-            # XXX: Better to write a separate function for this
-            return len(self.get_reduced_simplex(point, simplex, eps)) > 0
-        else:
-            # XXX: Better to write a separate function for this
-            return len(self.get_reduced_simplex(point, simplex, eps)) > 0
+        vertices = self.get_vertices(simplex)
+        return point_in_simplex(point, vertices, eps)
 
     def locate_point(self, point):
         """Find to which simplex the point belongs.
@@ -330,28 +360,7 @@ def circumscribed_circle(self, simplex, transform):
             The center and radius of the circumscribed circle
         """
         pts = np.dot(self.get_vertices(simplex), transform)
-        if self.dim == 2:
-            return fast_2d_circumcircle(pts)
-        if self.dim == 3:
-            return fast_3d_circumcircle(pts)
-
-        # Modified method from http://mathworld.wolfram.com/Circumsphere.html
-        mat = [[np.sum(np.square(pt)), *pt, 1] for pt in pts]
-
-        center = []
-        for i in range(1, len(simplex)):
-            r = np.delete(mat, i, 1)
-            factor = (-1) ** (i+1)
-            center.append(factor * np.linalg.det(r))
-
-        a = np.linalg.det(np.delete(mat, 0, 1))
-        center = [x / (2*a) for x in center]
-
-        x0 = self.vertices[next(iter(simplex))]
-        vec = np.subtract(center, x0)
-        radius = fast_norm(vec)
-
-        return tuple(center), radius
+        return circumsphere(pts)
 
     def point_in_cicumcircle(self, pt_index, simplex, transform):
         # return self.fast_point_in_circumcircle(pt_index, simplex, transform)

From 5f1eb3b3e9f16dae28eef9320cf9c628f20901c5 Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Tue, 17 Jul 2018 21:25:36 +0200
Subject: [PATCH 03/32] change choose_point_in_simplex such that it sometimes
 takes the center

---
 adaptive/learner/learnerND.py | 19 +++++++++++++------
 1 file changed, 13 insertions(+), 6 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index 9409d4b58..aece633b5 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -10,7 +10,7 @@
 from ..notebook_integration import ensure_holoviews
 from .base_learner import BaseLearner
 
-from .triangulation import Triangulation
+from .triangulation import Triangulation, point_in_simplex, circumsphere
 import random
 from math import factorial
 
@@ -126,12 +126,19 @@ def choose_point_in_simplex(simplex, transform=None):
     if transform is not None:
         simplex = np.dot(simplex, transform)
 
-    distances = scipy.spatial.distance.pdist(simplex)
-    distance_matrix = scipy.spatial.distance.squareform(distances)
-    i, j = np.unravel_index(np.argmax(distance_matrix), distance_matrix.shape)
+    # choose center iff the shape of the simplex is nice,
+    # otherwise the longest edge
+    center, radius = circumsphere(simplex)
+    if point_in_simplex(center, simplex):
+        point = np.average(simplex, axis=0)
+    else:
+        distances = scipy.spatial.distance.pdist(simplex)
+        distance_matrix = scipy.spatial.distance.squareform(distances)
+        i, j = np.unravel_index(np.argmax(distance_matrix), distance_matrix.shape)
+
+        point = (simplex[i, :] + simplex[j, :]) / 2
 
-    point = (simplex[i, :] + simplex[j, :]) / 2
-    return np.linalg.solve(transform, point)
+    return np.linalg.solve(transform, point)  # undo the transform
 
 
 class LearnerND(BaseLearner):

From b90a25d382442fdfc6547c8e65441aaf0f745ae0 Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Fri, 20 Jul 2018 09:59:54 +0200
Subject: [PATCH 04/32] pep8

---
 adaptive/learner/learnerND.py | 20 ++++++++++++--------
 1 file changed, 12 insertions(+), 8 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index aece633b5..c5ed9eb3e 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -40,7 +40,7 @@ def volume(simplex, ys=None):
 def orientation(simplex):
     matrix = np.subtract(simplex[:-1], simplex[-1])
     # See https://www.jstor.org/stable/2315353
-    sign, logdet = np.linalg.slogdet(matrix)
+    sign, _logdet = np.linalg.slogdet(matrix)
     return sign
 
 
@@ -128,13 +128,14 @@ def choose_point_in_simplex(simplex, transform=None):
 
     # choose center iff the shape of the simplex is nice,
     # otherwise the longest edge
-    center, radius = circumsphere(simplex)
+    center, _radius = circumsphere(simplex)
     if point_in_simplex(center, simplex):
         point = np.average(simplex, axis=0)
     else:
         distances = scipy.spatial.distance.pdist(simplex)
         distance_matrix = scipy.spatial.distance.squareform(distances)
-        i, j = np.unravel_index(np.argmax(distance_matrix), distance_matrix.shape)
+        i, j = np.unravel_index(np.argmax(distance_matrix), 
+                                distance_matrix.shape)
 
         point = (simplex[i, :] + simplex[j, :]) / 2
 
@@ -364,14 +365,16 @@ def _ask(self, n=1):
             if len(losses):
                 loss, simplex = heapq.heappop(losses)
 
-                assert self._simplex_exists(simplex), "all simplices in the heap should exist"
+                assert self._simplex_exists(simplex), \
+                    "all simplices in the heap should exist"
 
                 if simplex in self._subtriangulations:
                     subtri = self._subtriangulations[simplex]
                     loss_density = loss / self.tri.volume(simplex)
                     for pend_simplex in subtri.simplices:
                         pend_loss = subtri.volume(pend_simplex) * loss_density
-                        heapq.heappush(pending_losses, (pend_loss, simplex, pend_simplex))
+                        heapq.heappush(
+                            pending_losses, (pend_loss, simplex, pend_simplex))
                     continue
             else:
                 loss = 0
@@ -401,7 +404,8 @@ def _ask(self, n=1):
         return new_points[0], new_loss_improvements[0]
 
     def update_losses(self, to_delete: set, to_add: set):
-        pending_points_unbound = set()  # XXX: add the points outside the triangulation to this as well
+        # XXX: add the points outside the triangulation to this as well
+        pending_points_unbound = set()
 
         for simplex in to_delete:
             self._losses.pop(simplex, None)
@@ -466,8 +470,8 @@ def plot(self, n=None, tri_alpha=0):
             raise NotImplementedError('holoviews currently does not support',
                                       '3D surface plots in bokeh.')
         if len(self.bounds) != 2:
-           raise NotImplementedError("Only 2D plots are implemented: You can "
-                                     "plot a 2D slice with 'plot_slice'.")
+            raise NotImplementedError("Only 2D plots are implemented: You can "
+                                      "plot a 2D slice with 'plot_slice'.")
         x, y = self.bounds
         lbrt = x[0], y[0], x[1], y[1]
 

From 9ae74260a888ed8237e3beeb6c316b2487bf44f0 Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Fri, 20 Jul 2018 10:35:42 +0200
Subject: [PATCH 05/32] update docstring

---
 adaptive/learner/learnerND.py | 42 +++++++++++++++++++++++------------
 1 file changed, 28 insertions(+), 14 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index c5ed9eb3e..b2593c14d 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -58,16 +58,25 @@ def std_loss(simplex, ys):
 
 
 def simplex_volume(vertices) -> float:
-    """
-    Return the volume of the simplex with given vertices.
+    """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 
+    surface area of a triangle in a 3d space.
 
-    If vertices are given they must be in a arraylike with shape (N+1, N):
-    the position vectors of the N+1 vertices in N dimensions. If the sides
-    are given, they must be the compressed pairwise distance matrix as
-    returned from scipy.spatial.distance.pdist.
+    Parameters
+    ----------
+    vertices: arraylike (2 dimensional)
+        array of points
 
-    Raises a ValueError if the vertices do not form a simplex (for example,
-    because they are coplanar, colinear or coincident).
+    Returns
+    -------
+    volume: int
+        the volume of the simplex with given vertices.
+
+    Raises
+    ------
+    ValueError:
+        if the vertices do not form a simplex (for example,
+        because they are coplanar, colinear or coincident).
 
     Warning: this algorithm has not been tested for numerical stability.
     """
@@ -107,7 +116,9 @@ def default_loss(simplex, ys):
 def choose_point_in_simplex(simplex, transform=None):
     """Choose a new point in inside a simplex.
 
-    Pick the center of the longest edge of this simplex
+    Pick the center of the simplex if the shape is nice (that is, the 
+    circumcenter lies within the simplex). Otherwise take the middle of the 
+    longest edge.
 
     Parameters
     ----------
@@ -118,7 +129,7 @@ def choose_point_in_simplex(simplex, transform=None):
 
     Returns
     -------
-    point : numpy array
+    point : numpy array of length N
         The coordinates of the suggested new point.
     """
 
@@ -126,20 +137,23 @@ def choose_point_in_simplex(simplex, transform=None):
     if transform is not None:
         simplex = np.dot(simplex, transform)
 
-    # choose center iff the shape of the simplex is nice,
-    # otherwise the longest edge
+    # choose center if and only if the shape of the simplex is nice,
+    # otherwise: the longest edge
     center, _radius = circumsphere(simplex)
     if point_in_simplex(center, simplex):
         point = np.average(simplex, axis=0)
     else:
         distances = scipy.spatial.distance.pdist(simplex)
         distance_matrix = scipy.spatial.distance.squareform(distances)
-        i, j = np.unravel_index(np.argmax(distance_matrix), 
+        i, j = np.unravel_index(np.argmax(distance_matrix),
                                 distance_matrix.shape)
 
         point = (simplex[i, :] + simplex[j, :]) / 2
 
-    return np.linalg.solve(transform, point)  # undo the transform
+    if transform is not None:
+        point = np.linalg.solve(transform, point)  # undo the transform
+        
+    return point
 
 
 class LearnerND(BaseLearner):

From d73443c7926e92358ce270731b6414e20c0f826e Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Fri, 13 Jul 2018 15:57:55 +0200
Subject: [PATCH 06/32] remove n from _ask signature and split off
 _ask_bound_point

---
 adaptive/learner/learnerND.py | 29 +++++++++++++----------------
 1 file changed, 13 insertions(+), 16 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index ce6ec128d..fe97390d7 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -275,31 +275,28 @@ def ask(self, n=1):
         xs, losses = zip(*(self._ask() for _ in range(n)))
         return list(xs), list(losses)
 
-    def _ask(self, n=1):
+    def _ask_bound_point(self):
+        bounds_to_do = [p for p in self._bounds_points
+                        if p not in self.data and p not in self._pending]
+        new_point = bounds_to_do[0]
+        self._tell_pending(new_point)
+        return new_point, np.inf
+
+    def _ask(self):
         # Complexity: O(N log N)
         # XXX: adapt this function
         # XXX: allow cases where n > 1
-        assert n == 1
 
-        new_points = []
-        new_loss_improvements = []
         if not self.bounds_are_done:
-            bounds_to_do = [p for p in self._bounds_points
-                            if p not in self.data and p not in self._pending]
-            new_points = bounds_to_do[:n]
-            new_loss_improvements = [np.inf] * len(new_points)
-            n = n - len(new_points)
-            for p in new_points:
-                self._tell_pending(p)
-
-        if n == 0:
-            return new_points[0], new_loss_improvements[0]
+            return self._ask_bound_point()
 
         losses = [(-v, k) for k, v in self.losses().items()]
         heapq.heapify(losses)
-
         pending_losses = []  # also a heap
 
+        new_points = []
+        new_loss_improvements = []
+
         if len(losses) == 0:
             # pick a random point inside the bounds
             a = np.diff(self.bounds).flat
@@ -309,7 +306,7 @@ def _ask(self, n=1):
             p = tuple(p)
             return p, np.inf
 
-        while len(new_points) < n:
+        while len(new_points) < 1:
             if len(losses):
                 loss, simplex = heapq.heappop(losses)
 

From 23efe579338b350a99ba7b911ee1f07ba6dc68ed Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Fri, 13 Jul 2018 16:02:19 +0200
Subject: [PATCH 07/32] split _ask_point_without_known_simplices into separate
 function

---
 adaptive/learner/learnerND.py | 26 ++++++++++++++------------
 1 file changed, 14 insertions(+), 12 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index fe97390d7..67f5cc1cd 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -282,29 +282,31 @@ def _ask_bound_point(self):
         self._tell_pending(new_point)
         return new_point, np.inf
 
+    def _ask_point_without_known_simplices(self):
+        # pick a random point inside the bounds
+        # XXX: change this into picking a point based on volume loss
+        a = np.diff(self.bounds).flat
+        b = np.array(self.bounds)[:, 0]
+        r = np.array([self._random.random() for _ in range(self.ndim)])
+        p = r * a + b
+        p = tuple(p)
+        return p, np.inf
+
     def _ask(self):
         # Complexity: O(N log N)
-        # XXX: adapt this function
-        # XXX: allow cases where n > 1
-
         if not self.bounds_are_done:
             return self._ask_bound_point()
 
-        losses = [(-v, k) for k, v in self.losses().items()]
-        heapq.heapify(losses)
+        losses = [(-v, k) for k, v in self.losses().items()]  # O(N)
+        heapq.heapify(losses)  # O(N log N)
         pending_losses = []  # also a heap
 
         new_points = []
         new_loss_improvements = []
 
         if len(losses) == 0:
-            # pick a random point inside the bounds
-            a = np.diff(self.bounds).flat
-            b = np.array(self.bounds)[:, 0]
-            r = np.array([self._random.random() for _ in range(self.ndim)])
-            p = r * a + b
-            p = tuple(p)
-            return p, np.inf
+            # we have no known simplices
+            return self._ask_point_without_known_simplices()
 
         while len(new_points) < 1:
             if len(losses):

From ebc6dc0ecc773f9f6e1c7078c29e0527740d60d8 Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Mon, 16 Jul 2018 12:46:28 +0200
Subject: [PATCH 08/32] use SortedList instead of heap

---
 adaptive/learner/learnerND.py | 28 +++++++++++++++++-----------
 1 file changed, 17 insertions(+), 11 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index 67f5cc1cd..414a300d5 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -12,6 +12,7 @@
 
 from .triangulation import Triangulation
 import random
+from sortedcontainers import SortedList
 
 
 def find_initial_simplex(pts, ndim):
@@ -283,6 +284,7 @@ def _ask_bound_point(self):
         return new_point, np.inf
 
     def _ask_point_without_known_simplices(self):
+        assert self.bounds_are_done
         # pick a random point inside the bounds
         # XXX: change this into picking a point based on volume loss
         a = np.diff(self.bounds).flat
@@ -292,22 +294,15 @@ def _ask_point_without_known_simplices(self):
         p = tuple(p)
         return p, np.inf
 
-    def _ask(self):
-        # Complexity: O(N log N)
-        if not self.bounds_are_done:
-            return self._ask_bound_point()
-
-        losses = [(-v, k) for k, v in self.losses().items()]  # O(N)
-        heapq.heapify(losses)  # O(N log N)
+    def _ask_best_point(self):
+        assert self.bounds_are_done
+        assert self.tri is not None
+        losses = SortedList([(v, k) for k, v in self.losses().items()])  # O(N log N)
         pending_losses = []  # also a heap
 
         new_points = []
         new_loss_improvements = []
 
-        if len(losses) == 0:
-            # we have no known simplices
-            return self._ask_point_without_known_simplices()
-
         while len(new_points) < 1:
             if len(losses):
                 loss, simplex = heapq.heappop(losses)
@@ -348,6 +343,16 @@ def _ask(self):
 
         return new_points[0], new_loss_improvements[0]
 
+    def _ask(self):
+        if not self.bounds_are_done:
+            return self._ask_bound_point()  # O(1)
+
+        if self.tri is None:
+            # we have no known simplices
+            return self._ask_point_without_known_simplices()  # O(1)
+
+        return self._ask_best_point()  # O(N)
+
     def update_losses(self, to_delete: set, to_add: set):
         pending_points_unbound = set()  # XXX: add the points outside the triangulation to this as well
 
@@ -384,6 +389,7 @@ def losses(self):
         losses : dict
             the key is a simplex, the value is the loss of this simplex
         """
+        # XXX could be a property
         if self.tri is None:
             return dict()
 

From 65e04acc38d0a05acfc300a9be55769b8a7aaa4f Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Mon, 16 Jul 2018 20:56:21 +0200
Subject: [PATCH 09/32] lose indentation

---
 adaptive/learner/learnerND.py | 31 +++++++++++++++++++++----------
 1 file changed, 21 insertions(+), 10 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index 414a300d5..bc4798ecf 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -167,11 +167,20 @@ def __init__(self, func, bounds, loss_per_simplex=None):
         # triangulation of the pending points inside a specific simplex
         self._subtriangulations = dict()  # simplex -> triangulation
 
+        # scale to unit
         self._transform = np.linalg.inv(np.diag(np.diff(bounds).flat))
 
         # create a private random number generator with fixed seed
         self._random = random.Random(1)
 
+        # All real triangles
+        # list of tuple (loss, simplex)
+        self._losses_real = SortedList()
+
+        # all real triangles that have not been subdivided and the pending triangles
+        # list of tuple (loss, real simplex, sub_simplex or None)
+        self._losses_combined = SortedList()
+
     @property
     def npoints(self):
         return len(self.data)
@@ -263,14 +272,16 @@ def _tell_pending(self, point, simplex=None):
         # Neighbours also includes the simplex itself
 
         for simpl in neighbours:
-            if self.tri.point_in_simplex(point, simpl):
-                subtri = self._subtriangulations
-                if simpl not in subtri:
-                    vertices = self.tri.get_vertices(simpl)
-                    tr = subtri[simpl] = Triangulation(vertices)
-                    tr.add_point(point, next(iter(tr.simplices)))
-                else:
-                    subtri[simpl].add_point(point)
+            if not self.tri.point_in_simplex(point, simpl):
+                continue  # point not in simplex
+            
+            subtri = self._subtriangulations
+            if simpl not in subtri:
+                vertices = self.tri.get_vertices(simpl)
+                tr = subtri[simpl] = Triangulation(vertices)
+                tr.add_point(point, next(iter(tr.simplices)))
+            else:
+                subtri[simpl].add_point(point)
 
     def ask(self, n=1):
         xs, losses = zip(*(self._ask() for _ in range(n)))
@@ -305,7 +316,7 @@ def _ask_best_point(self):
 
         while len(new_points) < 1:
             if len(losses):
-                loss, simplex = heapq.heappop(losses)
+                loss, simplex = losses.pop()
 
                 assert self._simplex_exists(simplex), "all simplices in the heap should exist"
 
@@ -314,7 +325,7 @@ def _ask_best_point(self):
                     loss_density = loss / self.tri.volume(simplex)
                     for pend_simplex in subtri.simplices:
                         pend_loss = subtri.volume(pend_simplex) * loss_density
-                        heapq.heappush(pending_losses, (pend_loss, simplex, pend_simplex))
+                        heapq.heappush(pending_losses, (-pend_loss, simplex, pend_simplex))
                     continue
             else:
                 loss = 0

From d8f508aeaca4f6802164e2e6bc4374da763ea04b Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Tue, 17 Jul 2018 14:01:29 +0200
Subject: [PATCH 10/32] make learnerND.ask O(log N)

---
 adaptive/learner/learnerND.py | 138 +++++++++++++++++++---------------
 1 file changed, 78 insertions(+), 60 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index bc4798ecf..103b65bd5 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -1,7 +1,6 @@
 # -*- coding: utf-8 -*-
 from collections import OrderedDict
 import itertools
-import heapq
 
 import numpy as np
 from scipy import interpolate
@@ -175,11 +174,12 @@ def __init__(self, func, bounds, loss_per_simplex=None):
 
         # All real triangles
         # list of tuple (loss, simplex)
-        self._losses_real = SortedList()
+        self._losses_list = SortedList()
 
         # all real triangles that have not been subdivided and the pending triangles
         # list of tuple (loss, real simplex, sub_simplex or None)
-        self._losses_combined = SortedList()
+        self._losses_combined_list = SortedList()
+        self._sub_losses = dict()  # map the other way around
 
     @property
     def npoints(self):
@@ -203,6 +203,11 @@ def ip(self):
         # XXX: take our own triangulation into account when generating the ip
         return interpolate.LinearNDInterpolator(self.points, self.values)
 
+    def _get_loss(self, simplex, subsimplex=None):
+        if subsimplex is None:
+            return self._losses[simplex]
+        return self._sub_losses[(simplex, subsimplex)]
+
     @property
     def tri(self):
         if self._tri is not None:
@@ -265,6 +270,8 @@ def _tell_pending(self, point, simplex=None):
         simplex = tuple(simplex or self.tri.locate_point(point))
         if not simplex:
             return
+            # Simplex is None if pending point is outside the triangulation,
+            # then you do not have subtriangles
 
         simplex = tuple(simplex)
         simplices = [self.tri.vertex_to_simplices[i] for i in simplex]
@@ -274,14 +281,30 @@ def _tell_pending(self, point, simplex=None):
         for simpl in neighbours:
             if not self.tri.point_in_simplex(point, simpl):
                 continue  # point not in simplex
-            
-            subtri = self._subtriangulations
-            if simpl not in subtri:
+
+            if simpl not in self._subtriangulations:
                 vertices = self.tri.get_vertices(simpl)
-                tr = subtri[simpl] = Triangulation(vertices)
-                tr.add_point(point, next(iter(tr.simplices)))
+                tr = self._subtriangulations[simpl] = Triangulation(vertices)
+                todel, toadd = tr.add_point(point, next(iter(tr.simplices)))
             else:
-                subtri[simpl].add_point(point)
+                todel, toadd = self._subtriangulations[simpl].add_point(point)
+
+            for subsimplex in todel:
+                if subsimplex == tuple(range(self.ndim + 1)):
+                    continue
+                subloss = self._sub_losses.pop((simpl, subsimplex))
+                self._losses_combined_list.discard((subloss, simpl, subsimplex))
+
+            loss = self._losses[simpl]
+            self._losses_combined_list.discard((loss, simpl, None))
+
+            loss = self._losses[simpl]
+            loss_density = loss / self.tri.volume(simpl)
+            subtriangulation = self._subtriangulations[simpl]
+            for subsimplex in toadd:
+                subloss = subtriangulation.volume(subsimplex) * loss_density
+                self._losses_combined_list.add((subloss, simpl, subsimplex))
+                self._sub_losses[(simpl, subsimplex)] = subloss
 
     def ask(self, n=1):
         xs, losses = zip(*(self._ask() for _ in range(n)))
@@ -308,51 +331,25 @@ def _ask_point_without_known_simplices(self):
     def _ask_best_point(self):
         assert self.bounds_are_done
         assert self.tri is not None
-        losses = SortedList([(v, k) for k, v in self.losses().items()])  # O(N log N)
-        pending_losses = []  # also a heap
-
-        new_points = []
-        new_loss_improvements = []
-
-        while len(new_points) < 1:
-            if len(losses):
-                loss, simplex = losses.pop()
-
-                assert self._simplex_exists(simplex), "all simplices in the heap should exist"
-
-                if simplex in self._subtriangulations:
-                    subtri = self._subtriangulations[simplex]
-                    loss_density = loss / self.tri.volume(simplex)
-                    for pend_simplex in subtri.simplices:
-                        pend_loss = subtri.volume(pend_simplex) * loss_density
-                        heapq.heappush(pending_losses, (-pend_loss, simplex, pend_simplex))
-                    continue
-            else:
-                loss = 0
-                simplex = ()
-
-            points = np.array(self.tri.get_vertices(simplex))
-            loss = abs(loss)
-            if len(pending_losses):
-                pend_loss, real_simp, pend_simp = pending_losses[0]
-                pend_loss = abs(pend_loss)
 
-                if pend_loss > loss:
-                    subtri = self._subtriangulations[real_simp]
-                    points = np.array(subtri.get_vertices(pend_simp))
-                    simplex = real_simp
-                    loss = pend_loss
+        # _losses_combined is a SortedList. pop() will give us the highest loss
+        loss, simplex, subsimplex = self._losses_combined_list.pop()  # O(log N)
 
-            point_new = tuple(choose_point_in_simplex(
-                points, transform=self._transform))
-            self._pending_to_simplex[point_new] = simplex
+        if subsimplex is None:
+            # We found a real simplex and want to subdivide it
+            points = self.tri.get_vertices(simplex)
+        else:
+            # We found a pending simplex and want to subdivide it
+            subtri = self._subtriangulations[simplex]
+            points = subtri.get_vertices(subsimplex)
 
-            new_points.append(point_new)
-            new_loss_improvements.append(-loss)
+        point_new = tuple(choose_point_in_simplex(points,
+                                                  transform=self._transform))
 
-            self._tell_pending(point_new, simplex)
+        self._pending_to_simplex[point_new] = simplex
+        self._tell_pending(point_new, simplex)  # O(??)
 
-        return new_points[0], new_loss_improvements[0]
+        return point_new, loss
 
     def _ask(self):
         if not self.bounds_are_done:
@@ -362,16 +359,23 @@ def _ask(self):
             # we have no known simplices
             return self._ask_point_without_known_simplices()  # O(1)
 
-        return self._ask_best_point()  # O(N)
+        return self._ask_best_point()  # O(??)
 
     def update_losses(self, to_delete: set, to_add: set):
         pending_points_unbound = set()  # XXX: add the points outside the triangulation to this as well
 
         for simplex in to_delete:
-            self._losses.pop(simplex, None)
+            loss = self._losses.pop(simplex, None)
+            self._losses_list.discard((loss, simplex))
+            self._losses_combined_list.discard((loss, simplex, None))
             subtri = self._subtriangulations.pop(simplex, None)
-            if subtri is not None:
-                pending_points_unbound.update(subtri.vertices)
+            if subtri is None:
+                continue
+
+            pending_points_unbound.update(subtri.vertices)
+            for subsimplex in subtri.simplices:
+                subloss = self._sub_losses[(simplex, subsimplex)]
+                self._losses_combined_list.discard((subloss, simplex, subsimplex))
 
         pending_points_unbound = set(p for p in pending_points_unbound
                                      if p not in self.data)
@@ -381,16 +385,30 @@ def update_losses(self, to_delete: set, to_add: set):
             values = [self.data[tuple(v)] for v in vertices]
             loss = self.loss_per_simplex(vertices, values)
             self._losses[simplex] = float(loss)
+            self._losses_list.add((loss, simplex))
 
             for p in pending_points_unbound:
                 # try to insert it
-                if self.tri.point_in_simplex(p, simplex):
-                    subtri = self._subtriangulations
-                    if simplex not in subtri:
-                        vertices = self.tri.get_vertices(simplex)
-                        subtri[simplex] = Triangulation(vertices)
-                    subtri[simplex].add_point(p)
-                    self._pending_to_simplex[p] = simplex
+                if not self.tri.point_in_simplex(p, simplex):
+                    continue
+
+                if simplex not in self._subtriangulations:
+                    vertices = self.tri.get_vertices(simplex)
+                    self._subtriangulations[simplex] = Triangulation(vertices)
+
+                self._subtriangulations[simplex].add_point(p)
+                self._pending_to_simplex[p] = simplex
+
+            if simplex not in self._subtriangulations:
+                self._losses_combined_list.add((loss, simplex, None))
+                continue
+
+            loss_density = loss / self.tri.volume(simplex)
+            subtriangulation = self._subtriangulations[simplex]
+            for subsimplex in subtriangulation.simplices:
+                subloss = subtriangulation.volume(subsimplex) * loss_density
+                self._losses_combined_list.add((subloss, simplex, subsimplex))
+                self._sub_losses[(simplex, subsimplex)] = subloss
 
     def losses(self):
         """Get the losses of each simplex in the current triangulation, as dict

From fd047a06b2e4fa7de761a90a8cc7d5d13ab33ccb Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Tue, 17 Jul 2018 14:34:40 +0200
Subject: [PATCH 11/32] improve triangulation performance in 3d

---
 adaptive/learner/triangulation.py | 15 +++++++++++----
 1 file changed, 11 insertions(+), 4 deletions(-)

diff --git a/adaptive/learner/triangulation.py b/adaptive/learner/triangulation.py
index b1c669231..59b7f34b8 100644
--- a/adaptive/learner/triangulation.py
+++ b/adaptive/learner/triangulation.py
@@ -9,6 +9,8 @@ def fast_norm(v):
     # notice this method can be even more optimised
     if len(v) == 2:
         return math.sqrt(v[0] * v[0] + v[1] * v[1])
+    if len(v) == 3:
+        return math.sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2])
     return math.sqrt(np.dot(v, v))
 
 
@@ -81,9 +83,12 @@ def fast_3d_circumcircle(points):
     points = np.array(points)
     pts = points[1:] - points[0]
 
-    l1, l2, l3 = [np.dot(p, p) for p in pts]  # length squared
     (x1, y1, z1), (x2, y2, z2), (x3, y3, z3) = pts
 
+    l1 = x1 * x1 + y1 * y1 + z1 * z1
+    l2 = x2 * x2 + y2 * y2 + z2 * z2
+    l3 = x3 * x3 + y3 * y3 + z3 * z3
+
     # Compute some determinants:
     dx = (+ l1 * (y2 * z3 - z2 * y3)
           - l2 * (y1 * z3 - z1 * y3)
@@ -99,11 +104,13 @@ def fast_3d_circumcircle(points):
           + x3 * (y1 * z2 - z1 * y2))
     a = 2 * aa
 
-    center = [dx / a, -dy / a, dz / a]
+    center = (dx / a, -dy / a, dz / a)
     radius = fast_norm(center)
-    center = np.add(center, points[0])
+    center = (center[0] + points[0][0],
+              center[1] + points[0][1],
+              center[2] + points[0][2])
 
-    return tuple(center), radius
+    return center, radius
 
 
 def orientation(face, origin):

From fce140aa53a789dfdf0797b04accb2f244dfc798 Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Tue, 17 Jul 2018 15:09:18 +0200
Subject: [PATCH 12/32] make learner pass the tests

---
 adaptive/learner/learnerND.py     | 14 +++++++++-----
 adaptive/learner/triangulation.py |  2 +-
 2 files changed, 10 insertions(+), 6 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index 103b65bd5..04302c052 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -155,7 +155,7 @@ def __init__(self, func, bounds, loss_per_simplex=None):
         self.data = OrderedDict()
         self._pending = set()
 
-        self._bounds_points = list(itertools.product(*bounds))
+        self._bounds_points = list(map(tuple, itertools.product(*bounds)))
 
         self.function = func
         self._tri = None
@@ -318,7 +318,7 @@ def _ask_bound_point(self):
         return new_point, np.inf
 
     def _ask_point_without_known_simplices(self):
-        assert self.bounds_are_done
+        assert not self._bounds_available
         # pick a random point inside the bounds
         # XXX: change this into picking a point based on volume loss
         a = np.diff(self.bounds).flat
@@ -329,7 +329,7 @@ def _ask_point_without_known_simplices(self):
         return p, np.inf
 
     def _ask_best_point(self):
-        assert self.bounds_are_done
+        assert not self._bounds_available
         assert self.tri is not None
 
         # _losses_combined is a SortedList. pop() will give us the highest loss
@@ -351,8 +351,12 @@ def _ask_best_point(self):
 
         return point_new, loss
 
+    @property
+    def _bounds_available(self):
+        return any((p not in self._pending and p not in self.data) for p in self._bounds_points)
+
     def _ask(self):
-        if not self.bounds_are_done:
+        if self._bounds_available:
             return self._ask_bound_point()  # O(1)
 
         if self.tri is None:
@@ -383,7 +387,7 @@ def update_losses(self, to_delete: set, to_add: set):
         for simplex in to_add:
             vertices = self.tri.get_vertices(simplex)
             values = [self.data[tuple(v)] for v in vertices]
-            loss = self.loss_per_simplex(vertices, values)
+            loss = float(self.loss_per_simplex(vertices, values))
             self._losses[simplex] = float(loss)
             self._losses_list.add((loss, simplex))
 
diff --git a/adaptive/learner/triangulation.py b/adaptive/learner/triangulation.py
index 59b7f34b8..f2eadc6cd 100644
--- a/adaptive/learner/triangulation.py
+++ b/adaptive/learner/triangulation.py
@@ -487,7 +487,7 @@ def volume(self, simplex):
         prefactor = np.math.factorial(self.dim)
         vertices = np.array(self.get_vertices(simplex))
         vectors = vertices[1:] - vertices[0]
-        return abs(np.linalg.det(vectors)) / prefactor
+        return float(abs(np.linalg.det(vectors)) / prefactor)
 
     def volumes(self):
         return [self.volume(sim) for sim in self.simplices]

From bded62bf47a927937d93dca567f2f51c8a07d546 Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Fri, 20 Jul 2018 17:18:25 +0200
Subject: [PATCH 13/32] use heaps instead of sortedlist

---
 adaptive/learner/learnerND.py | 71 +++++++++++++++--------------------
 1 file changed, 30 insertions(+), 41 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index 04302c052..8844486fe 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -11,7 +11,7 @@
 
 from .triangulation import Triangulation
 import random
-from sortedcontainers import SortedList
+import heapq
 
 
 def find_initial_simplex(pts, ndim):
@@ -172,14 +172,9 @@ def __init__(self, func, bounds, loss_per_simplex=None):
         # create a private random number generator with fixed seed
         self._random = random.Random(1)
 
-        # All real triangles
-        # list of tuple (loss, simplex)
-        self._losses_list = SortedList()
-
-        # all real triangles that have not been subdivided and the pending triangles
-        # list of tuple (loss, real simplex, sub_simplex or None)
-        self._losses_combined_list = SortedList()
-        self._sub_losses = dict()  # map the other way around
+        # 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
 
     @property
     def npoints(self):
@@ -203,11 +198,6 @@ def ip(self):
         # XXX: take our own triangulation into account when generating the ip
         return interpolate.LinearNDInterpolator(self.points, self.values)
 
-    def _get_loss(self, simplex, subsimplex=None):
-        if subsimplex is None:
-            return self._losses[simplex]
-        return self._sub_losses[(simplex, subsimplex)]
-
     @property
     def tri(self):
         if self._tri is not None:
@@ -285,26 +275,19 @@ def _tell_pending(self, point, simplex=None):
             if simpl not in self._subtriangulations:
                 vertices = self.tri.get_vertices(simpl)
                 tr = self._subtriangulations[simpl] = Triangulation(vertices)
-                todel, toadd = tr.add_point(point, next(iter(tr.simplices)))
+                _, toadd = tr.add_point(point, next(iter(tr.simplices)))
             else:
-                todel, toadd = self._subtriangulations[simpl].add_point(point)
-
-            for subsimplex in todel:
-                if subsimplex == tuple(range(self.ndim + 1)):
-                    continue
-                subloss = self._sub_losses.pop((simpl, subsimplex))
-                self._losses_combined_list.discard((subloss, simpl, subsimplex))
+                _, toadd = self._subtriangulations[simpl].add_point(point)
 
             loss = self._losses[simpl]
-            self._losses_combined_list.discard((loss, simpl, None))
 
-            loss = self._losses[simpl]
             loss_density = loss / self.tri.volume(simpl)
             subtriangulation = self._subtriangulations[simpl]
             for subsimplex in toadd:
                 subloss = subtriangulation.volume(subsimplex) * loss_density
-                self._losses_combined_list.add((subloss, simpl, subsimplex))
-                self._sub_losses[(simpl, subsimplex)] = subloss
+                heapq.heappush(self._losses_combined, 
+                               (-subloss, simpl, subsimplex))
+                
 
     def ask(self, n=1):
         xs, losses = zip(*(self._ask() for _ in range(n)))
@@ -328,12 +311,26 @@ def _ask_point_without_known_simplices(self):
         p = tuple(p)
         return p, np.inf
 
+    def _pop_highest_existing_simplex(self):
+        # find the simplex with the highest loss, we do need to check that the
+        # simplex hasn't been deleted yet
+        while True:
+            loss, simplex, subsimplex = heapq.heappop(self._losses_combined) 
+            if (subsimplex is None and 
+                    simplex in self.tri.simplices and 
+                    simplex not in self._subtriangulations):
+                return abs(loss), simplex, subsimplex
+            if (simplex in self._subtriangulations and
+                    simplex in self.tri.simplices and 
+                    subsimplex in self._subtriangulations[simplex].simplices):
+                return abs(loss), simplex, subsimplex
+
+
     def _ask_best_point(self):
         assert not self._bounds_available
         assert self.tri is not None
 
-        # _losses_combined is a SortedList. pop() will give us the highest loss
-        loss, simplex, subsimplex = self._losses_combined_list.pop()  # O(log N)
+        loss, simplex, subsimplex = self._pop_highest_existing_simplex()
 
         if subsimplex is None:
             # We found a real simplex and want to subdivide it
@@ -370,16 +367,9 @@ def update_losses(self, to_delete: set, to_add: set):
 
         for simplex in to_delete:
             loss = self._losses.pop(simplex, None)
-            self._losses_list.discard((loss, simplex))
-            self._losses_combined_list.discard((loss, simplex, None))
             subtri = self._subtriangulations.pop(simplex, None)
-            if subtri is None:
-                continue
-
-            pending_points_unbound.update(subtri.vertices)
-            for subsimplex in subtri.simplices:
-                subloss = self._sub_losses[(simplex, subsimplex)]
-                self._losses_combined_list.discard((subloss, simplex, subsimplex))
+            if subtri is not None:
+                pending_points_unbound.update(subtri.vertices)
 
         pending_points_unbound = set(p for p in pending_points_unbound
                                      if p not in self.data)
@@ -389,7 +379,6 @@ def update_losses(self, to_delete: set, to_add: set):
             values = [self.data[tuple(v)] for v in vertices]
             loss = float(self.loss_per_simplex(vertices, values))
             self._losses[simplex] = float(loss)
-            self._losses_list.add((loss, simplex))
 
             for p in pending_points_unbound:
                 # try to insert it
@@ -404,15 +393,15 @@ def update_losses(self, to_delete: set, to_add: set):
                 self._pending_to_simplex[p] = simplex
 
             if simplex not in self._subtriangulations:
-                self._losses_combined_list.add((loss, simplex, None))
+                heapq.heappush(self._losses_combined, (-loss, simplex, None))
                 continue
 
             loss_density = loss / self.tri.volume(simplex)
             subtriangulation = self._subtriangulations[simplex]
             for subsimplex in subtriangulation.simplices:
                 subloss = subtriangulation.volume(subsimplex) * loss_density
-                self._losses_combined_list.add((subloss, simplex, subsimplex))
-                self._sub_losses[(simplex, subsimplex)] = subloss
+                heapq.heappush(self._losses_combined, 
+                               (-subloss, simplex, subsimplex))
 
     def losses(self):
         """Get the losses of each simplex in the current triangulation, as dict

From 79ddeff09e8d17346ee13ec4527cbb60a559f364 Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Fri, 20 Jul 2018 17:34:43 +0200
Subject: [PATCH 14/32] add tell_pending to _ask_point_without_known_simplices

---
 adaptive/learner/learnerND.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index 8844486fe..a05eaee8d 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -309,6 +309,8 @@ def _ask_point_without_known_simplices(self):
         r = np.array([self._random.random() for _ in range(self.ndim)])
         p = r * a + b
         p = tuple(p)
+
+        self._tell_pending(p)
         return p, np.inf
 
     def _pop_highest_existing_simplex(self):
@@ -360,7 +362,7 @@ def _ask(self):
             # we have no known simplices
             return self._ask_point_without_known_simplices()  # O(1)
 
-        return self._ask_best_point()  # O(??)
+        return self._ask_best_point()  # O(log N)
 
     def update_losses(self, to_delete: set, to_add: set):
         pending_points_unbound = set()  # XXX: add the points outside the triangulation to this as well

From d136d19effdfbf830e40aabd29cd5ae5b031899b Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Fri, 20 Jul 2018 18:04:49 +0200
Subject: [PATCH 15/32] make bowyer_watson faster by looking at less
 circumcircles

it prunes some circumcircles faster
it results in:
~20% faster in 2d
~40% faster in 3d
---
 adaptive/learner/triangulation.py | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/adaptive/learner/triangulation.py b/adaptive/learner/triangulation.py
index f2eadc6cd..fb0a5663f 100644
--- a/adaptive/learner/triangulation.py
+++ b/adaptive/learner/triangulation.py
@@ -411,15 +411,18 @@ def bowyer_watson(self, pt_index, containing_simplex=None, transform=None):
 
             if self.point_in_cicumcircle(pt_index, simplex, transform):
                 self.delete_simplex(simplex)
-                todo_points = set(simplex) - done_points
+                todo_points = set(simplex)
                 done_points.update(simplex)
 
-                if len(todo_points):
-                    neighbours = set.union(*[self.vertex_to_simplices[p]
-                                             for p in todo_points])
-                    queue.update(neighbours - done_simplices)
-
                 bad_triangles.add(simplex)
+                if len(todo_points) == 0:
+                    continue
+                neighbours = set.union(*[self.vertex_to_simplices[p]
+                                            for p in todo_points])
+                neighbours = neighbours - done_simplices
+                neighbours = set(simpl for simpl in neighbours if len(set(simpl) & done_points) >= self.dim)
+                queue.update(neighbours)
+
 
         faces = list(self.faces(simplices=bad_triangles))
 

From d1c51924f97f4c45dd6c5b48c8aaf4a8f68d2c57 Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Mon, 23 Jul 2018 12:05:53 +0200
Subject: [PATCH 16/32] move simplex_volume function to triangulation file

---
 adaptive/learner/learnerND.py     | 56 +++----------------------------
 adaptive/learner/triangulation.py | 49 +++++++++++++++++++++++++++
 2 files changed, 53 insertions(+), 52 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index b2593c14d..71ebd1598 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -10,9 +10,9 @@
 from ..notebook_integration import ensure_holoviews
 from .base_learner import BaseLearner
 
-from .triangulation import Triangulation, point_in_simplex, circumsphere
+from .triangulation import Triangulation, point_in_simplex, \
+                           circumsphere, simplex_volume_in_embedding
 import random
-from math import factorial
 
 
 def find_initial_simplex(pts, ndim):
@@ -57,60 +57,13 @@ def std_loss(simplex, ys):
     return r.flat * np.power(vol, 1. / dim) + vol
 
 
-def simplex_volume(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 
-    surface area of a triangle in a 3d space.
-
-    Parameters
-    ----------
-    vertices: arraylike (2 dimensional)
-        array of points
-
-    Returns
-    -------
-    volume: int
-        the volume of the simplex with given vertices.
-
-    Raises
-    ------
-    ValueError:
-        if the vertices do not form a simplex (for example,
-        because they are coplanar, colinear or coincident).
-
-    Warning: this algorithm has not been tested for numerical stability.
-    """
-
-    # Implements http://mathworld.wolfram.com/Cayley-MengerDeterminant.html
-    # Modified from https://codereview.stackexchange.com/questions/77593/calculating-the-volume-of-a-tetrahedron
-
-    # β_ij = |v_i - v_k|²
-    vertices = np.array(vertices, dtype=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 = np.concatenate((np.ones(num_verts), 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
-    vol_square = np.linalg.det(sq_dists_mat) / coeff
-
-    if vol_square <= 0:
-        raise ValueError('Provided vertices do not form a tetrahedron')
-
-    return np.sqrt(vol_square)
-
-
 def default_loss(simplex, ys):
     # return std_loss(simplex, ys)
     if isinstance(ys[0], Iterable):
         pts = [(*s, *ys[i]) for i, s in enumerate(simplex)]
     else:
         pts = [(*s, ys[i]) for i, s in enumerate(simplex)]
-    return simplex_volume(pts)
+    return simplex_volume_in_embedding(pts)
 
 
 def choose_point_in_simplex(simplex, transform=None):
@@ -133,12 +86,11 @@ def choose_point_in_simplex(simplex, transform=None):
         The coordinates of the suggested new point.
     """
 
-    # XXX: find a better selection algorithm
     if transform is not None:
         simplex = np.dot(simplex, transform)
 
     # choose center if and only if the shape of the simplex is nice,
-    # otherwise: the longest edge
+    # otherwise: the center the longest edge
     center, _radius = circumsphere(simplex)
     if point_in_simplex(center, simplex):
         point = np.average(simplex, axis=0)
diff --git a/adaptive/learner/triangulation.py b/adaptive/learner/triangulation.py
index df7df1cf5..d70ac1795 100644
--- a/adaptive/learner/triangulation.py
+++ b/adaptive/learner/triangulation.py
@@ -3,6 +3,8 @@
 
 import numpy as np
 import math
+import scipy.spatial
+from math import factorial
 
 
 def fast_norm(v):
@@ -173,6 +175,53 @@ def is_iterable_and_sized(obj):
     return isinstance(obj, Iterable) and isinstance(obj, Sized)
 
 
+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 
+    surface area of a triangle in a 3d space.
+
+    Parameters
+    ----------
+    vertices: arraylike (2 dimensional)
+        array of points
+
+    Returns
+    -------
+    volume: int
+        the volume of the simplex with given vertices.
+
+    Raises
+    ------
+    ValueError:
+        if the vertices do not form a simplex (for example,
+        because they are coplanar, colinear or coincident).
+
+    Warning: this algorithm has not been tested for numerical stability.
+    """
+
+    # Implements http://mathworld.wolfram.com/Cayley-MengerDeterminant.html
+    # Modified from https://codereview.stackexchange.com/questions/77593/calculating-the-volume-of-a-tetrahedron
+
+    # β_ij = |v_i - v_k|²
+    vertices = np.array(vertices, dtype=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 = np.concatenate((np.ones(num_verts), 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
+    vol_square = np.linalg.det(sq_dists_mat) / coeff
+
+    if vol_square <= 0:
+        raise ValueError('Provided vertices do not form a tetrahedron')
+
+    return np.sqrt(vol_square)
+
+
 class Triangulation:
     """A triangulation object.
 

From 01c53f4959d886cd9c55ff098e2e1d690bd58e1a Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Mon, 23 Jul 2018 12:13:39 +0200
Subject: [PATCH 17/32] merge lists with zip

---
 adaptive/learner/learnerND.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index 71ebd1598..f7deebc58 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -60,9 +60,9 @@ def std_loss(simplex, ys):
 def default_loss(simplex, ys):
     # return std_loss(simplex, ys)
     if isinstance(ys[0], Iterable):
-        pts = [(*s, *ys[i]) for i, s in enumerate(simplex)]
+        pts = [(*x, *y) for x, y in zip(simplex, ys)]
     else:
-        pts = [(*s, ys[i]) for i, s in enumerate(simplex)]
+        pts = [(*x, y) for x, y in zip(simplex, ys)]
     return simplex_volume_in_embedding(pts)
 
 

From c1d59173a9b1a4ba8da0e1e58268ee63610bfb07 Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Mon, 23 Jul 2018 12:14:36 +0200
Subject: [PATCH 18/32] tetrahedron -> simplex

---
 adaptive/learner/triangulation.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/adaptive/learner/triangulation.py b/adaptive/learner/triangulation.py
index d70ac1795..38a799771 100644
--- a/adaptive/learner/triangulation.py
+++ b/adaptive/learner/triangulation.py
@@ -217,7 +217,7 @@ def simplex_volume_in_embedding(vertices) -> float:
     vol_square = np.linalg.det(sq_dists_mat) / coeff
 
     if vol_square <= 0:
-        raise ValueError('Provided vertices do not form a tetrahedron')
+        raise ValueError('Provided vertices do not form a simplex')
 
     return np.sqrt(vol_square)
 

From 7921d1de1fb6055e96bd50c924fe81270b299ef0 Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Mon, 23 Jul 2018 14:12:44 +0200
Subject: [PATCH 19/32] add quick volume computation for dim == 2

---
 adaptive/learner/triangulation.py | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/adaptive/learner/triangulation.py b/adaptive/learner/triangulation.py
index 38a799771..8237d9eaf 100644
--- a/adaptive/learner/triangulation.py
+++ b/adaptive/learner/triangulation.py
@@ -202,8 +202,16 @@ def simplex_volume_in_embedding(vertices) -> float:
     # Implements http://mathworld.wolfram.com/Cayley-MengerDeterminant.html
     # Modified from https://codereview.stackexchange.com/questions/77593/calculating-the-volume-of-a-tetrahedron
 
-    # β_ij = |v_i - v_k|²
+    
     vertices = np.array(vertices, dtype=float)
+    dim = len(vertices[0])
+    if dim == 2:
+        # Heron's formula
+        a, b, c = scipy.spatial.distance.pdist(vertices, metric='euclidean')
+        s = 0.5 * (a + b + c)
+        return math.sqrt(s*(s-a)*(s-b)*(s-c))
+
+    # β_ij = |v_i - v_k|²
     sq_dists = scipy.spatial.distance.pdist(vertices, metric='sqeuclidean')
 
     # Add border while compressed

From 74601c6c790fec16c26982fa788ba9514ad89d96 Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Mon, 23 Jul 2018 21:23:05 +0200
Subject: [PATCH 20/32] try to add anisotropicity, however, this fails...

---
 adaptive/learner/learnerND.py | 69 +++++++++++++++++++++++++++++++----
 1 file changed, 62 insertions(+), 7 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index f7deebc58..a2c1d2462 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -14,6 +14,9 @@
                            circumsphere, simplex_volume_in_embedding
 import random
 
+import rtree
+import math
+
 
 def find_initial_simplex(pts, ndim):
     origin = pts[0]
@@ -165,7 +168,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
@@ -184,11 +187,15 @@ def __init__(self, func, bounds, loss_per_simplex=None):
         # triangulation of the pending points inside a specific simplex
         self._subtriangulations = dict()  # simplex -> triangulation
 
-        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)
 
+        props = rtree.index.Property(dimension=self.ndim)
+        self._point_tree = rtree.index.Index(properties=props)
+        self.anisotripic = anisotropic
+        
     @property
     def npoints(self):
         return len(self.data)
@@ -228,7 +235,7 @@ def tri(self):
                   if i not in initial_simplex]
 
         for p in to_add:
-            self._tri.add_point(p, transform=self._transform)
+            self._tri.add_point(p, transform=self._scale_matrix)
         # XXX: also compute losses of initial simplex
 
     @property
@@ -241,7 +248,7 @@ def points(self):
 
     def tell(self, point, value):
         point = tuple(point)
-
+        print(point)
         if point in self.data:
             return
 
@@ -249,16 +256,64 @@ def tell(self, point, value):
             return self._tell_pending(point)
 
         self._pending.discard(point)
+        # Insert the point into our Rtree
+        scaled_point = tuple(np.dot(self._scale_matrix, point))
+        # print(scaled_point)
+        self._point_tree.insert(self.npoints, scaled_point)  
         self.data[point] = value
 
         if self.tri is not None:
             simplex = self._pending_to_simplex.get(point)
             if simplex is not None and not self._simplex_exists(simplex):
                 simplex = None
+            transform = self.get_local_transform_matrix(point)
             to_delete, to_add = self._tri.add_point(
-                point, simplex, transform=self._transform)
+                point, simplex, transform=transform)
             self.update_losses(to_delete, to_add)
 
+    def get_local_transform_matrix(self, point):
+        scale = self._scale_matrix
+        if self.tri is None or not self.anisotripic:
+            return scale
+        
+        # Get some points in the neighbourhood
+        scaled_point = tuple(np.dot(scale, point))
+        indices = self._point_tree.nearest(scaled_point, 10)
+        
+        points = [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
+        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)
+        *gradient, _constant = fit
+        
+        # we do not need the constant
+        # 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) / 10
+
+        scale_along_gradient = projection_matrix * factor + identity
+        m = np.dot(scale_along_gradient, scale)
+        print("m:", m)
+        return m
+
+
     def _simplex_exists(self, simplex):
         simplex = tuple(sorted(simplex))
         return simplex in self.tri.simplices
@@ -359,7 +414,7 @@ def _ask(self, n=1):
                     loss = pend_loss
 
             point_new = tuple(choose_point_in_simplex(
-                points, transform=self._transform))
+                points, transform=self._scale_matrix))
             self._pending_to_simplex[point_new] = simplex
 
             new_points.append(point_new)
@@ -502,7 +557,7 @@ def plot_slice(self, cut_mapping, n=None):
             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))
 

From b1fc0bd3f0f40123007c38563b820c22450099c5 Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Wed, 25 Jul 2018 13:27:07 +0200
Subject: [PATCH 21/32] make anisotropicity work

---
 adaptive/learner/learnerND.py     | 22 +++++++++++-----------
 adaptive/learner/triangulation.py | 18 ++++++++++++++----
 2 files changed, 25 insertions(+), 15 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index 56aec1c9d..4111432cf 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -240,6 +240,7 @@ def tri(self):
         for p in to_add:
             self._tri.add_point(p, transform=self._scale_matrix)
         # XXX: also compute losses of initial simplex
+        self.update_losses(set(), self._tri.simplices)
 
     @property
     def values(self):
@@ -251,7 +252,6 @@ def points(self):
 
     def tell(self, point, value):
         point = tuple(point)
-        print(point)
         if point in self.data:
             return
 
@@ -276,6 +276,7 @@ def tell(self, point, value):
 
     def get_local_transform_matrix(self, point):
         scale = self._scale_matrix
+
         if self.tri is None or not self.anisotripic:
             return scale
         
@@ -283,7 +284,7 @@ def get_local_transform_matrix(self, point):
         scaled_point = tuple(np.dot(scale, point))
         indices = self._point_tree.nearest(scaled_point, 10)
         
-        points = [self.points[i] for i in indices]
+        points = np.array([self.points[i] for i in indices])
         values = [self.data[tuple(p)] for p in points]
 
         if isinstance(values[0], Iterable):
@@ -291,12 +292,13 @@ def get_local_transform_matrix(self, point):
 
         # 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)
         *gradient, _constant = fit
-        
+
         # we do not need the constant
         # gradient is a vector of the amount of the slope in each direction
         magnitude = np.linalg.norm(gradient)
@@ -309,13 +311,12 @@ def get_local_transform_matrix(self, point):
         projection_matrix = (gradient.T @ gradient)
         identity = np.eye(self.ndim)
 
-        factor = math.sqrt(magnitude ** 2 + 1) / 10
+        factor = math.sqrt(magnitude ** 2 + 1) - 1
+        factor = min(factor, 2)
 
-        scale_along_gradient = projection_matrix * factor + identity
+        scale_along_gradient = projection_matrix * factor * 1 + identity
         m = np.dot(scale_along_gradient, scale)
-        print("m:", m)
-        return m
-
+        return m.T
 
     def _simplex_exists(self, simplex):
         simplex = tuple(sorted(simplex))
@@ -358,7 +359,6 @@ def _tell_pending(self, point, simplex=None):
                 subloss = subtriangulation.volume(subsimplex) * loss_density
                 heapq.heappush(self._losses_combined, 
                                (-subloss, simpl, subsimplex))
-                
 
     def ask(self, n=1):
         xs, losses = zip(*(self._ask() for _ in range(n)))
@@ -412,8 +412,8 @@ def _ask_best_point(self):
             subtri = self._subtriangulations[simplex]
             points = subtri.get_vertices(subsimplex)
 
-        point_new = tuple(choose_point_in_simplex(points,
-                                                  transform=self._scale_matrix))
+        transform = self.get_local_transform_matrix(points[0])
+        point_new = tuple(choose_point_in_simplex(points, transform=transform))
 
         self._pending_to_simplex[point_new] = simplex
         self._tell_pending(point_new, simplex)  # O(??)
diff --git a/adaptive/learner/triangulation.py b/adaptive/learner/triangulation.py
index 64fd528f6..3d0073a05 100644
--- a/adaptive/learner/triangulation.py
+++ b/adaptive/learner/triangulation.py
@@ -485,11 +485,10 @@ def bowyer_watson(self, pt_index, containing_simplex=None, transform=None):
                     continue
                 neighbours = set.union(*[self.vertex_to_simplices[p]
                                             for p in todo_points])
-                neighbours = neighbours - done_simplices
-                neighbours = set(simpl for simpl in neighbours if len(set(simpl) & done_points) >= self.dim)
+                neighbours = neighbours - done_simplices - queue
+                neighbours = set(simpl for simpl in neighbours if len(set(simpl) & set(simplex)) >= self.dim)
                 queue.update(neighbours)
 
-
         faces = list(self.faces(simplices=bad_triangles))
 
         multiplicities = Counter(face for face in faces)
@@ -502,7 +501,18 @@ def bowyer_watson(self, pt_index, containing_simplex=None, transform=None):
                 self.add_simplex((*face, pt_index))
 
         new_triangles = self.vertex_to_simplices[pt_index]
-        return bad_triangles - new_triangles, new_triangles - bad_triangles
+        deleted_simplices = bad_triangles - new_triangles
+        new_simplices = new_triangles - bad_triangles
+
+        if not np.isclose(sum([self.volume(s) for s in deleted_simplices]),
+                          sum([self.volume(s) for s in new_simplices])):
+            print("vertices = ", self.vertices, "\n\n\n")
+            print("simplices = ", self.simplices, "\n\n\n")
+            print("new_s = ", new_simplices)
+            print("old_s = ", deleted_simplices)
+            print("holes_faces = ", hole_faces)
+            raise AssertionError("debug me")
+        return deleted_simplices, new_simplices
 
     def add_point(self, point, simplex=None, transform=None):
         """Add a new vertex and create simplices as appropriate.

From 971043096daddc077074f7d8d48b5f6a4f2d4192 Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Mon, 20 Aug 2018 11:09:41 +0200
Subject: [PATCH 22/32] out-comment the maximum shape factor of two

---
 adaptive/learner/learnerND.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index 4111432cf..6b26b1099 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -312,7 +312,7 @@ def get_local_transform_matrix(self, point):
         identity = np.eye(self.ndim)
 
         factor = math.sqrt(magnitude ** 2 + 1) - 1
-        factor = min(factor, 2)
+        # factor = min(factor, 2)
 
         scale_along_gradient = projection_matrix * factor * 1 + identity
         m = np.dot(scale_along_gradient, scale)

From 1df1216aed0c8ed0c2ad8f3eff0b693b7cef3b7a Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Mon, 20 Aug 2018 13:49:57 +0200
Subject: [PATCH 23/32] remove dependence on RTree

---
 adaptive/learner/learnerND.py     | 31 ++++++++++---------------------
 adaptive/learner/triangulation.py |  8 --------
 2 files changed, 10 insertions(+), 29 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index 6b26b1099..06a050e9e 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -13,8 +13,6 @@
                            circumsphere, simplex_volume_in_embedding
 import random
 import heapq
-
-import rtree
 import math
 
 
@@ -190,9 +188,6 @@ def __init__(self, func, bounds, loss_per_simplex=None, anisotropic=False):
 
         # create a private random number generator with fixed seed
         self._random = random.Random(1)
-
-        props = rtree.index.Property(dimension=self.ndim)
-        self._point_tree = rtree.index.Index(properties=props)
         self.anisotripic = anisotropic
         
         # all real triangles that have not been subdivided and the pending
@@ -259,31 +254,26 @@ def tell(self, point, value):
             return self._tell_pending(point)
 
         self._pending.discard(point)
-        # Insert the point into our Rtree
-        scaled_point = tuple(np.dot(self._scale_matrix, point))
-        # print(scaled_point)
-        self._point_tree.insert(self.npoints, scaled_point)  
         self.data[point] = value
 
         if self.tri is not None:
             simplex = self._pending_to_simplex.get(point)
             if simplex is not None and not self._simplex_exists(simplex):
-                simplex = None
-            transform = self.get_local_transform_matrix(point)
+                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, point):
+    def get_local_transform_matrix(self, simplex):
         scale = self._scale_matrix
 
-        if self.tri is None or not self.anisotripic:
+        if simplex is None or self.tri is None or not self.anisotripic:
             return scale
         
-        # Get some points in the neighbourhood
-        scaled_point = tuple(np.dot(scale, point))
-        indices = self._point_tree.nearest(scaled_point, 10)
-        
+        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]
 
@@ -298,8 +288,8 @@ def get_local_transform_matrix(self, point):
         B = np.array(values, dtype=float)
         fit, *_ = np.linalg.lstsq(A, B, rcond=None)
         *gradient, _constant = fit
+        # we do not need the constant, only the gradient
 
-        # we do not need the constant
         # gradient is a vector of the amount of the slope in each direction
         magnitude = np.linalg.norm(gradient)
 
@@ -312,9 +302,8 @@ def get_local_transform_matrix(self, point):
         identity = np.eye(self.ndim)
 
         factor = math.sqrt(magnitude ** 2 + 1) - 1
-        # factor = min(factor, 2)
 
-        scale_along_gradient = projection_matrix * factor * 1 + identity
+        scale_along_gradient = projection_matrix * factor + identity
         m = np.dot(scale_along_gradient, scale)
         return m.T
 
@@ -412,7 +401,7 @@ def _ask_best_point(self):
             subtri = self._subtriangulations[simplex]
             points = subtri.get_vertices(subsimplex)
 
-        transform = self.get_local_transform_matrix(points[0])
+        transform = self.get_local_transform_matrix(simplex)
         point_new = tuple(choose_point_in_simplex(points, transform=transform))
 
         self._pending_to_simplex[point_new] = simplex
diff --git a/adaptive/learner/triangulation.py b/adaptive/learner/triangulation.py
index 3d0073a05..6966905de 100644
--- a/adaptive/learner/triangulation.py
+++ b/adaptive/learner/triangulation.py
@@ -504,14 +504,6 @@ def bowyer_watson(self, pt_index, containing_simplex=None, transform=None):
         deleted_simplices = bad_triangles - new_triangles
         new_simplices = new_triangles - bad_triangles
 
-        if not np.isclose(sum([self.volume(s) for s in deleted_simplices]),
-                          sum([self.volume(s) for s in new_simplices])):
-            print("vertices = ", self.vertices, "\n\n\n")
-            print("simplices = ", self.simplices, "\n\n\n")
-            print("new_s = ", new_simplices)
-            print("old_s = ", deleted_simplices)
-            print("holes_faces = ", hole_faces)
-            raise AssertionError("debug me")
         return deleted_simplices, new_simplices
 
     def add_point(self, point, simplex=None, transform=None):

From 672247c3d91c2026e9bc05de25b42ee1142906af Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Thu, 27 Sep 2018 20:32:19 +0200
Subject: [PATCH 24/32] enable plotting of custom triangulations

---
 adaptive/learner/learnerND.py     |  8 ++++--
 adaptive/learner/triangulation.py | 47 +++++++++++++++++++++++++++++++
 2 files changed, 52 insertions(+), 3 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index fd717f0e8..ad0473215 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -12,8 +12,9 @@
 
 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):
@@ -217,8 +218,9 @@ def bounds_are_done(self):
         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)
+        fd = FakeDelaunay(self.tri.vertices, self.tri.simplices)
+        values = list(map(lambda k: self.data[k], self.tri.vertices))
+        return interpolate.LinearNDInterpolator(fd, values) 
 
     @property
     def tri(self):
diff --git a/adaptive/learner/triangulation.py b/adaptive/learner/triangulation.py
index 41ebed3ee..7f57576a2 100644
--- a/adaptive/learner/triangulation.py
+++ b/adaptive/learner/triangulation.py
@@ -622,3 +622,50 @@ def hull(self):
     def convex_invariant(self, vertex):
         """Hull is convex."""
         raise NotImplementedError
+
+
+class FakeDelaunay(scipy.spatial.Delaunay):
+    def __init__(self, vertices, simplices):
+        """
+        Throw in a Triangulation object, and receive a brandnew 
+        scipy.spatial.Delaunay
+
+        Parameters
+        ----------
+        vertices : 2d arraylike of floats
+
+        simplices : 2d arraylike of integers
+        """
+
+        # just create some Delaunay triangulation object, any would do.
+        super().__init__([[0,1], [1,0], [0,0]])
+
+        self.ndim = len(vertices[0])
+        self._points = np.array(vertices, dtype='float64')
+        self.simplices = np.array(list(simplices), dtype='int32')
+        # self.neighbors = np.zero # () todo
+        self.nsimplex = len(simplices)
+        self.npoints = len(vertices)
+        self.min_bound = self._points.min(axis=0)
+        self.max_bound = self._points.max(axis=0)
+        self._transform = None # Set it to None to recompute it when needed
+
+        # now the hard part: finding the neighbours
+        self.neighbors = -np.ones_like(self.simplices) # set all to -1
+
+        faces_dict = {}
+        for n, simplex in enumerate(self.simplices):
+            simplex = tuple(simplex)
+            for i, p in enumerate(simplex):
+                face = tuple(sorted(simplex[:i] + simplex[i+1:]))
+                if face in faces_dict:
+                    other_simpl_index, other_point_index = faces_dict[face]
+                    self.neighbors[other_simpl_index, other_point_index] = n
+                    self.neighbors[n, i] = other_simpl_index
+                else:
+                    faces_dict[face] = (n, i)
+
+
+
+
+

From d777dc1eeaac8c3f9c6a0ad3c3d8f740e9033691 Mon Sep 17 00:00:00 2001
From: Bas Nijholt <basnijholt@gmail.com>
Date: Fri, 28 Sep 2018 00:02:57 +0200
Subject: [PATCH 25/32] small style change

---
 adaptive/learner/learnerND.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index ad0473215..941f856f4 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -218,9 +218,9 @@ def bounds_are_done(self):
         return all(p in self.data for p in self._bounds_points)
 
     def ip(self):
-        fd = FakeDelaunay(self.tri.vertices, self.tri.simplices)
-        values = list(map(lambda k: self.data[k], self.tri.vertices))
-        return interpolate.LinearNDInterpolator(fd, 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):
@@ -290,7 +290,7 @@ def get_local_transform_matrix(self, simplex):
         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)
+        fit = np.linalg.lstsq(A, B, rcond=None)[0]
         *gradient, _constant = fit
         # we do not need the constant, only the gradient
 

From d9d981e1453464d7d8b098e832df8e8222b5db4e Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Sun, 30 Sep 2018 16:14:04 +0200
Subject: [PATCH 26/32] change inside bounds check to have eps precision

---
 adaptive/learner/learnerND.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index 941f856f4..881e14486 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -315,8 +315,8 @@ 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)

From f5eb8ef30f7439738b4663265010252904c73204 Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Sun, 30 Sep 2018 16:14:48 +0200
Subject: [PATCH 27/32] add checks for volume consistency

---
 adaptive/learner/triangulation.py |  4 ++++
 adaptive/tests/test_learnernd.py  | 26 ++++++++++++++++++++++++--
 2 files changed, 28 insertions(+), 2 deletions(-)

diff --git a/adaptive/learner/triangulation.py b/adaptive/learner/triangulation.py
index 7f57576a2..c36ff1aa7 100644
--- a/adaptive/learner/triangulation.py
+++ b/adaptive/learner/triangulation.py
@@ -520,6 +520,10 @@ def bowyer_watson(self, pt_index, containing_simplex=None, transform=None):
         deleted_simplices = bad_triangles - new_triangles
         new_simplices = new_triangles - bad_triangles
 
+        old_vol = sum([self.volume(simplex) for simplex in deleted_simplices])
+        new_vol = sum([self.volume(simplex) for simplex in new_simplices])
+        assert np.isclose(old_vol, new_vol), f"{old_vol} !== {new_vol}"
+
         return deleted_simplices, new_simplices
 
     def add_point(self, point, simplex=None, transform=None):
diff --git a/adaptive/tests/test_learnernd.py b/adaptive/tests/test_learnernd.py
index 0f28ecab2..3ffcf9995 100644
--- a/adaptive/tests/test_learnernd.py
+++ b/adaptive/tests/test_learnernd.py
@@ -1,10 +1,17 @@
 # -*- coding: utf-8 -*-
 
 from ..learner import LearnerND
-from ..runner import replay_log
+from ..runner import replay_log, BlockingRunner
+import numpy as np
 
 
-def test_faiure_case_LearnerND():
+def sphere(xyz):
+    x, y, z = xyz
+    a = 0.4
+    return x + z**2 + np.exp(-(x**2 + y**2 + z**2 - 0.75**2)**2/a**4)
+
+
+def test_failure_case_LearnerND():
     log = [
         ('ask', 4),
         ('tell', (-1, -1, -1), 1.607873907219222e-101),
@@ -21,3 +28,18 @@ def test_faiure_case_LearnerND():
     ]
     learner = LearnerND(lambda *x: x, bounds=[(-1, 1), (-1, 1), (-1, 1)])
     replay_log(learner, log)
+
+
+def test_anisotropic_3d():
+    # there was this bug that the total volume would exceed the bounding box 
+    # volume for the anisotropic 3d learnerND
+    # learner = adaptive.LearnerND(ring, bounds=[(-1, 1), (-1, 1)], anisotropic=True)
+    learner = LearnerND(sphere, bounds=[(-1, 1), (-1, 1), (-1, 1)], anisotropic=True)
+    def goal(l):
+        if l.tri:
+            sum_of_simplex_volumes = np.sum(l.tri.volumes())
+            assert sum_of_simplex_volumes < 8.00000000001
+        return l.npoints >= 1000
+    BlockingRunner(learner, goal, ntasks=1)
+
+    assert learner.npoints >= 1000

From ff8845acf6aa8a8bc578f4d84174872380ad0e84 Mon Sep 17 00:00:00 2001
From: Jorn Hoofwijk <jornhoofwijk@gmail.com>
Date: Mon, 1 Oct 2018 11:31:41 +0200
Subject: [PATCH 28/32] add .pytest_cache to gitignore

---
 .gitignore | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.gitignore b/.gitignore
index 89fe43568..c70fa195a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -44,6 +44,7 @@ nosetests.xml
 coverage.xml
 *,cover
 .hypothesis/
+.pytest_cache/
 
 # Translations
 *.mo

From 2efbfb6df20acfed3e637dc41c9d4218a7531975 Mon Sep 17 00:00:00 2001
From: Bas Nijholt <bas@nijho.lt>
Date: Fri, 3 Apr 2026 13:08:53 -0700
Subject: [PATCH 29/32] fix: Python 3.10+ compatibility (collections.abc,
 math.factorial)

---
 adaptive/learner/learnerND.py     | 5 +++--
 adaptive/learner/triangulation.py | 5 +++--
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index 881e14486..0a10f6f00 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -1,5 +1,6 @@
 # -*- coding: utf-8 -*-
-from collections import OrderedDict, Iterable
+from collections import OrderedDict
+from collections.abc import Iterable
 import heapq
 import itertools
 import random
@@ -24,7 +25,7 @@ def volume(simplex, ys=None):
     dim = len(simplex) - 1
 
     # See https://www.jstor.org/stable/2315353
-    vol = np.abs(np.linalg.det(matrix)) / np.math.factorial(dim)
+    vol = np.abs(np.linalg.det(matrix)) / math.factorial(dim)
     return vol
 
 
diff --git a/adaptive/learner/triangulation.py b/adaptive/learner/triangulation.py
index c36ff1aa7..caec070fa 100644
--- a/adaptive/learner/triangulation.py
+++ b/adaptive/learner/triangulation.py
@@ -1,4 +1,5 @@
-from collections import Counter, Sized, Iterable
+from collections import Counter
+from collections.abc import Sized, Iterable
 from itertools import combinations, chain
 
 import numpy as np
@@ -575,7 +576,7 @@ def add_point(self, point, simplex=None, transform=None):
             return self.bowyer_watson(pt_index, actual_simplex, transform)
 
     def volume(self, simplex):
-        prefactor = np.math.factorial(self.dim)
+        prefactor = factorial(self.dim)
         vertices = np.array(self.get_vertices(simplex))
         vectors = vertices[1:] - vertices[0]
         return float(abs(np.linalg.det(vectors)) / prefactor)

From f84cab565b91425469d89d2b7a6df8efc584a150 Mon Sep 17 00:00:00 2001
From: Bas Nijholt <bas@nijho.lt>
Date: Fri, 3 Apr 2026 13:27:33 -0700
Subject: [PATCH 30/32] fix: use uniform scale for Bowyer-Watson insertion, not
 anisotropic transform

The per-simplex anisotropic transform (based on local gradient) varies
across the mesh, which can produce disconnected or non-star-shaped
cavities in the Bowyer-Watson algorithm, breaking its volume
conservation invariant.

Fix: use the uniform bounds-scaling matrix for triangulation
maintenance (Bowyer-Watson insertion). The anisotropic transform is
still used for point selection (choose_point_in_simplex in
_ask_best_point), so anisotropic sampling behavior is preserved.

Closes #74.
---
 adaptive/learner/learnerND.py | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index 0a10f6f00..da8ac2338 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -265,9 +265,14 @@ def tell(self, point, value):
             simplex = self._pending_to_simplex.get(point)
             if simplex is not None and not self._simplex_exists(simplex):
                 simplex = self.tri.locate_point(point)
-            transform = self.get_local_transform_matrix(simplex)
+            # Use uniform scale matrix for Bowyer-Watson insertion rather than
+            # the per-simplex anisotropic transform. The anisotropic transform
+            # varies locally (based on gradient), which can produce disconnected
+            # or non-star-shaped cavities in the original space, breaking the
+            # volume conservation invariant of Bowyer-Watson. Anisotropy still
+            # drives point selection via choose_point_in_simplex in _ask_best_point.
             to_delete, to_add = self._tri.add_point(
-                point, simplex, transform=transform)
+                point, simplex, transform=self._scale_matrix)
             self.update_losses(to_delete, to_add)
 
     def get_local_transform_matrix(self, simplex):

From 94214ca89d301531fb553fb17e5be35195cf924e Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Sun, 5 Apr 2026 20:43:28 +0000
Subject: [PATCH 31/32] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 docs/source/logo.md | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/docs/source/logo.md b/docs/source/logo.md
index 5cb740c2f..33939990b 100644
--- a/docs/source/logo.md
+++ b/docs/source/logo.md
@@ -4,13 +4,11 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.5
+    jupytext_version: 1.19.1
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
   name: python3
-execution:
-  timeout: 300
 ---
 
 ```{code-cell} ipython3
@@ -168,7 +166,7 @@ def animate_png(folder=None, nseconds=15):
 
 
 def save_webp(fname_webp, ims):
-    (im, *_ims) = ims
+    im, *_ims = ims
     im.save(
         fname_webp,
         save_all=True,

From 4339e65ea691d97ae585a217033118856c84e581 Mon Sep 17 00:00:00 2001
From: Bas Nijholt <bas@nijho.lt>
Date: Sun, 5 Apr 2026 17:07:46 -0700
Subject: [PATCH 32/32] Make anisotropic keyword-only in LearnerND

---
 adaptive/learner/learnerND.py         | 2 +-
 adaptive/tests/unit/test_learnernd.py | 5 +++++
 2 files changed, 6 insertions(+), 1 deletion(-)

diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py
index ac192f2a1..eeb94656c 100644
--- a/adaptive/learner/learnerND.py
+++ b/adaptive/learner/learnerND.py
@@ -308,7 +308,7 @@ class LearnerND(BaseLearner):
     children based on volume.
     """
 
-    def __init__(self, func, bounds, loss_per_simplex=None, anisotropic=False):
+    def __init__(self, func, bounds, loss_per_simplex=None, *, anisotropic=False):
         self._vdim = None
         self.loss_per_simplex = loss_per_simplex or default_loss
 
diff --git a/adaptive/tests/unit/test_learnernd.py b/adaptive/tests/unit/test_learnernd.py
index ecd00d6da..def3e5908 100644
--- a/adaptive/tests/unit/test_learnernd.py
+++ b/adaptive/tests/unit/test_learnernd.py
@@ -34,6 +34,11 @@ def test_learnerND_accepts_ConvexHull_as_input():
     assert np.allclose(learner._bbox, [(0, 2), (0, 1)])
 
 
+def test_learnerND_anisotropic_is_keyword_only():
+    with pytest.raises(TypeError):
+        LearnerND(ring_of_fire, [(-1, 1), (-1, 1)], None, True)
+
+
 def test_learnerND_raises_if_too_many_neigbors():
     @uses_nth_neighbors(2)
     def loss(*args):