@@ -14,17 +14,6 @@ from .triangulation import Triangulation

 import random

-def find_initial_simplex(pts, ndim):

-    origin = pts[0]

-    vecs = pts[1:] - origin

-    if np.linalg.matrix_rank(vecs) < ndim:

-        return None  # all points are coplanar

-

-    tri = scipy.spatial.Delaunay(pts)

-    simplex = tri.simplices[0]  # take a random simplex from tri

-    return simplex

-

-

 def volume(simplex, ys=None):

     # Notice the parameter ys is there so you can use this volume method as

     # as loss function

@@ -198,19 +187,15 @@ class LearnerND(BaseLearner):

         if self._tri is not None:

             return self._tri

-        if len(self.data) < 2:

-            return None

-

-        initial_simplex = find_initial_simplex(self.points, self.ndim)

-        if initial_simplex is None:

+        try:

+            self._tri = Triangulation(self.points)

+            return self._tri

+        except ValueError:

+            # A ValueError is raised if we do not have enough points or

+            # the provided points are coplanar, so we need more points to create

+            # a valid triangulation

             return None

-        pts = self.points[initial_simplex]

-        self._tri = Triangulation([tuple(p) for p in pts])

-        to_add = [p for i, p in enumerate(self.points)

-                  if i not in initial_simplex]

-        for p in to_add:

-            self._tri.add_point(p, transform=self._transform)

         # XXX: also compute losses of initial simplex

     @property

@@ -225,19 +210,20 @@ class LearnerND(BaseLearner):

         point = tuple(point)

         if point in self.data:

-            return

+            return  # we already know about the point

         if value is None:

             return self._tell_pending(point)

         self._pending.discard(point)

+        tri = self.tri

         self.data[point] = value

-        if self.tri is not None:

+        if tri is not None:

             simplex = self._pending_to_simplex.get(point)

             if simplex is not None and not self._simplex_exists(simplex):

                 simplex = None

-            to_delete, to_add = self._tri.add_point(

+            to_delete, to_add = tri.add_point(

                 point, simplex, transform=self._transform)

             self.update_losses(to_delete, to_add)

@@ -313,14 +299,16 @@ class LearnerND(BaseLearner):

             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

@@ -350,7 +338,8 @@ class LearnerND(BaseLearner):

         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)

@@ -415,8 +404,8 @@ class LearnerND(BaseLearner):

             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]

@@ -3,6 +3,7 @@ from itertools import combinations, chain

 import numpy as np

 import math

+import scipy.spatial

 def fast_norm(v):

@@ -142,7 +143,7 @@ class Triangulation:

     Parameters

     ----------

     coords : 2d array-like of floats

-        Coordinates of vertices of the first simplex.

+        Coordinates of vertices.

     Attributes

     ----------

@@ -155,14 +156,23 @@ class Triangulation:

         index of the list.

     hull : set of int

         Exterior vertices

+

+    Raises

+    ------

+    ValueError

+        if the list of coordinates is incorrect or the points do not form one 

+        or more simplices in the 

     """

     def __init__(self, coords):

         if not is_iterable_and_sized(coords):

-            raise ValueError("Please provide a 2-dimensional list of points")

+            raise TypeError("Please provide a 2-dimensional list of points")

         coords = list(coords)

         if not all(is_iterable_and_sized(coord) for coord in coords):

-            raise ValueError("Please provide a 2-dimensional list of points")

+            raise TypeError("Please provide a 2-dimensional list of points")

+        if len(coords) == 0:

+            raise ValueError("Please provide at least one simplex") 

+            # raise now because otherwise the next line will raise a less

         dim = len(coords[0])

         if any(len(coord) != dim for coord in coords):

@@ -171,12 +181,12 @@ class Triangulation:

         if dim == 1:

             raise ValueError("Triangulation class only supports dim >= 2")

-        if len(coords) != dim + 1:

-            raise ValueError("Can only add one simplex on initialization")

+        if len(coords) < dim + 1:

+            raise ValueError("Please provide at least one simplex")

         coords = list(map(tuple, coords))

         vectors = np.subtract(coords[1:], coords[0])

-        if np.linalg.det(vectors) == 0:

+        if np.linalg.matrix_rank(vectors) < dim:

             raise ValueError("Initial simplex has zero volumes "

                              "(the points are linearly dependent)")

@@ -184,7 +194,12 @@ class Triangulation:

         self.simplices = set()

         # initialise empty set for each vertex

         self.vertex_to_simplices = [set() for _ in coords]

-        self.add_simplex(range(len(self.vertices)))

+

+        # find a Delaunay triangulation to start with, then we will throw it 

+        # away and continue with our own algorithm

+        initial_tri = scipy.spatial.Delaunay(coords)

+        for simplex in initial_tri.simplices:

+            self.add_simplex(simplex)

     def delete_simplex(self, simplex):

         simplex = tuple(sorted(simplex))

@@ -75,7 +75,7 @@ def test_triangulation_raises_exception_for_1d_list():

     # We could support 1d, but we don't for now, because it is not relevant

     # so a user has to be aware

     pts = [0, 1]

-    with pytest.raises(ValueError):

+    with pytest.raises(TypeError):

         Triangulation(pts)

@@ -292,3 +292,40 @@ def test_triangulation_is_deterministic(dim):

     t1 = _make_triangulation(points)

     t2 = _make_triangulation(points)

     assert t1.simplices == t2.simplices

+

+

+@with_dimension

+def test_initialisation_raises_when_not_enough_points(dim):

+    deficient_simplex = _make_standard_simplex(dim)[:-1]

+    

+    with pytest.raises(ValueError):

+        Triangulation(deficient_simplex)

+

+

+@with_dimension

+def test_initialisation_raises_when_points_coplanar(dim):

+    zero_volume_simplex = _make_standard_simplex(dim)[:-1]

+    

+    new_point1 = np.average(zero_volume_simplex, axis=0)

+    new_point2 = np.sum(zero_volume_simplex, axis=0)

+    zero_volume_simplex = np.vstack((zero_volume_simplex, 

+                                     new_point1, new_point2))

+    

+    with pytest.raises(ValueError):

+        Triangulation(zero_volume_simplex)

+

+

+@with_dimension

+def test_initialisation_accepts_more_than_one_simplex(dim):

+    points = _make_standard_simplex(dim)

+    new_point = [1.1] * dim  # Point oposing the origin but outside circumsphere

+    points = np.vstack((points, new_point))

+

+    tri = Triangulation(points)

+

+    simplex1 = tuple(range(dim+1))

+    simplex2 = tuple(range(1, dim+2))

+

+    _check_triangulation_is_valid(tri)    

+    

+    assert tri.simplices == {simplex1, simplex2}