WIP: unify Learner1D, Learner2D and LearnerND

adaptive/learner/triangulation.py

@@ -266,7 +266,7 @@ class Triangulation:
  or more simplices in the
  """
 
- def __init__(self, coords):
+ def __init__(self, coords, *, _check_vertices=True):
  if not is_iterable_and_sized(coords):
  raise TypeError("Please provide a 2-dimensional list of points")
  coords = list(coords)
@@ -287,23 +287,28 @@ def __init__(self, coords):
  raise ValueError("Please provide at least one simplex")
 
  coords = list(map(tuple, coords))
- vectors = np.subtract(coords[1:], coords[0])
- if np.linalg.matrix_rank(vectors) < dim:
- raise ValueError(
- "Initial simplex has zero volumes "
- "(the points are linearly dependent)"
- )
+ if _check_vertices:
+ vectors = np.subtract(coords[1:], coords[0])
+ if np.linalg.matrix_rank(vectors) < dim:
+ raise ValueError(
+ "Initial simplex has zero volumes "
+ "(the points are linearly dependent)"
+ )
 
  self.vertices = list(coords)
  self.simplices = set()
  # initialise empty set for each vertex
  self.vertex_to_simplices = [set() for _ in coords]
 
- # 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)
+ if len(coords) == dim + 1:
+ # There is just a single simplex
+ self.add_simplex(tuple(range(dim + 1)))
+ else:
+ # 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))

adaptive/priority_queue.py

@@ -0,0 +1,115 @@
+from sortedcontainers import SortedDict, SortedList
+
+__all__ = ["Empty", "Queue"]
+
+
+class Empty(KeyError):
+ pass
+
+
+class Queue:
+ """Priority queue supporting update and removal at arbitrary position.
+
+ Parameters
+ ----------
+ entries : iterable of (item, priority)
+ The initial data in the queue. Providing this is faster than
+ calling 'insert' a bunch of times.
+ """
+
+ def __init__(self, entries=()):
+ self._queue = SortedDict(
+ ((priority, -n), item) for n, (item, priority) in enumerate(entries)
+ )
+ # 'self._queue' cannot be keyed only on priority, as there may be several
+ # items that have the same priority. To keep unique elements the key
+ # will be '(priority, self._n)', where 'self._n' is decremented whenever
+ # we add a new element. 'self._n' is negative so that elements with equal
+ # priority are sorted by insertion order.
+ self._n = -len(self._queue)
+ # To efficiently support updating and removing items if their priority
+ # is unknown we have to keep the reverse map of 'self._queue'. Because
+ # items may not be hashable we cannot use a SortedDict, so we use a
+ # SortedList storing '(item, key)'.
+ self._items = SortedList(((v, k) for k, v in self._queue.items()))
+
+ def __len__(self):
+ return len(self._queue)
+
+ def items(self):
+ "Return an iterator over the items in the queue in priority order."
+ return reversed(self._queue.values())
+
+ def peek(self):
+ """Return the item and priority at the front of the queue.
+
+ Raises
+ ------
+ Empty : if the queue is empty
+ """
+ self._check_nonempty()
+ ((priority, _), item) = self._queue.peekitem()
+ return item, priority
+
+ def pop(self):
+ """Remove and return the item and priority at the front of the queue.
+
+ Raises
+ ------
+ Empty : if the queue is empty
+ """
+ self._check_nonempty()
+ (key, item) = self._queue.popitem()
+ i = self._items.index((item, key))
+ del self._items[i]
+ priority, _ = key
+ return item, priority
+
+ def insert(self, item, priority):
+ "Insert 'item' into the queue with the given priority."
+ key = (priority, self._n)
+ self._items.add((item, key))
+ self._queue[key] = item
+ self._n -= 1
+
+ def _check_nonempty(self):
+ if not self._queue:
+ raise Empty()
+
+ def _find_first(self, item):
+ self._check_nonempty()
+ i = self._items.bisect_left((item, ()))
+ try:
+ should_be, key = self._items[i]
+ except IndexError:
+ raise KeyError("item is not in queue") from None
+ if item != should_be:
+ raise KeyError("item is not in queue")
+ return i, key
+
+ def remove(self, item):
+ """Remove the 'item' from the queue.
+
+ Raises
+ ------
+ KeyError : if 'item' is not in the queue.
+ """
+ i, key = self._find_first(item)
+ del self._queue[key]
+ del self._items[i]
+
+ def update(self, item, priority):
+ """Update 'item' in the queue to have the given priority.
+
+ Raises
+ ------
+ KeyError : if 'item' is not in the queue.
+ """
+ i, key = self._find_first(item)
+ _, n = key
+ new_key = (priority, n)
+
+ del self._queue[key]
+ del self._items[i]
+ self._queue[new_key] = item
+ self._items.add((item, new_key))

adaptive/tests/domain_utils.py

@@ -0,0 +1,151 @@
+import itertools
+
+import numpy as np
+
+import hypothesis.strategies as st
+from adaptive.learner.new_learnerND import ConvexHull, Interval
+
+# This module contains utilities for producing domains and points inside and outside of them.
+# Because we typically do not want to test very degenerate cases (e.g. points that are almost
+# coincident, very large or very small) we prefer generating points in the interval [0, 1)
+# using numpy.random, rather than drawing from Hypothesis' "floats" strategy.
+
+
+# Return an iterator that yields matrices reflecting in the cartesian
+# coordinate axes in 'ndim' dimensions.
+def reflections(ndim):
+ return map(np.diag, itertools.product([1, -1], repeat=ndim))
+
+
+def point_inside_simplex(simplex):
+ simplex = np.asarray(simplex)
+ dim = simplex.shape[1]
+ # Generate a point in the unit simplex.
+ # https://cs.stackexchange.com/questions/3227/uniform-sampling-from-a-simplex
+ # We avoid using Hypothesis to generate the points as it typically chooses
+ # very annoying points, which we want to avoid testing for now.
+ xb = np.random.rand(dim)
+ xb = np.array(sorted(xb))
+ xb[1:] = xb[1:] - xb[:-1]
+ # Transform into the simplex we need
+ v0, vecs = simplex[0], simplex[1:] - simplex[0]
+ x = tuple(v0 + (vecs.T @ xb))
+ return x
+
+
+@st.composite
+def points_inside(draw, domain, n):
+ # Set the numpy random seed
+ draw(st.random_module())
+ if isinstance(domain, Interval):
+ a, b = domain.bounds
+ return a + (b - a) * np.random.rand(n)
+ else:
+ assert isinstance(domain, ConvexHull)
+ tri = domain.triangulation
+ simplices = list(tri.simplices)
+ simplex = st.sampled_from(simplices).map(
+ lambda simplex: [tri.vertices[s] for s in simplex]
+ )
+ # "point_inside_simplex" uses the numpy RNG, and we set the seed above.
+ # Together this means we're almost guaranteed not to get coinciding points.
+ # Note that we draw from the 'simplex' strategy on each iteration, so we
+ # distribute the points between the different simplices in the domain.
+ return [tuple(point_inside_simplex(draw(simplex))) for _ in range(n)]
+
+
+@st.composite
+def point_inside(draw, domain):
+ return draw(points_inside(domain, 1))[0]
+
+
+@st.composite
+def a_few_points_inside(draw, domain):
+ n = draw(st.integers(3, 20))
+ return draw(points_inside(domain, n))
+
+
+@st.composite
+def points_outside(draw, domain, n):
+ # set numpy random seed
+ draw(st.random_module())
+
+ if isinstance(domain, Interval):
+ a, b = domain.bounds
+ ndim = 1
+ else:
+ assert isinstance(domain, ConvexHull)
+ hull = domain.bounds
+ points = hull.points[hull.vertices]
+ ndim = points.shape[1]
+ a, b = points.min(axis=0)[None, :], points.max(axis=0)[None, :]
+
+ # Generate a point outside the bounding box of the domain.
+ center = (a + b) / 2
+ border = (b - a) / 2
+ r = border + 10 * border * np.random.rand(n, ndim)
+ quadrant = np.sign(np.random.rand(n, ndim) - 0.5)
+ assert not np.any(quadrant == 0)
+ return center + quadrant * r
+
+
+@st.composite
+def point_outside(draw, domain):
+ return draw(points_outside(domain, 1))[0]
+
+
+@st.composite
+def point_on_shared_face(draw, domain, dim):
+ # Return a point that is shared by at least 2 subdomains
+ assert isinstance(domain, ConvexHull)
+ assert 0 < dim < domain.ndim
+
+ # Set the numpy random seed
+ draw(st.random_module())
+
+ tri = domain.triangulation
+
+ for face in tri.faces(dim + 1):
+ containing_subdomains = tri.containing(face)
+ if len(containing_subdomains) > 1:
+ break
+
+ vertices = np.array([tri.vertices[i] for i in face])
+
+ xb = np.random.rand(dim)
+
+ x = tuple(vertices[0] + xb @ (vertices[1:] - vertices[0]))
+
+ assert all(tri.point_in_simplex(x, s) for s in containing_subdomains)
+
+ return x
+
+
+@st.composite
+def make_random_domain(draw, ndim, fill=True):
+ # Set the numpy random seed
+ draw(st.random_module())
+
+ if ndim == 1:
+ a, b = sorted(np.random.rand(2) - 0.5)
+ domain = Interval(a, b)
+ else:
+ # Generate points in a hypercube around the origin
+ points = np.random.rand(10, ndim) - 0.5
+ domain = ConvexHull(points)
+ return domain
+
+
+@st.composite
+def make_hypercube_domain(draw, ndim, fill=True):
+ # Set the numpy random seed
+ draw(st.random_module())
+ limit = np.random.rand()
+
+ if ndim == 1:
+ subdomain = Interval(-limit, limit)
+ else:
+ point = np.full(ndim, limit)
+ boundary_points = [r @ point for r in reflections(ndim)]
+ subdomain = ConvexHull(boundary_points)
+ return subdomain