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")
+        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,175 @@
+import itertools
+
+import numpy as np
+import scipy.linalg
+import scipy.spatial
+
+import hypothesis.strategies as st
+from adaptive.learner.new_learnerND import ConvexHull, Interval
+from hypothesis.extra import numpy as hynp
+
+
+def reflections(ndim):
+    return map(np.diag, itertools.product([1, -1], repeat=ndim))
+
+
+reals = st.floats(min_value=-100, max_value=100, allow_nan=False, allow_infinity=False)
+positive_reals = st.floats(
+    min_value=1e-3, max_value=100, allow_nan=False, allow_infinity=False
+)
+
+
+@st.composite
+def point(draw, ndim):
+    return draw(reals if ndim == 1 else st.tuples(*[reals] * ndim))
+
+
+def unique_vectors(xs):
+    xs = np.asarray(xs)
+    if len(xs.shape) == 1:
+        xs = xs[:, None]
+    c = np.max(np.linalg.norm(xs, axis=1))
+    if c == 0:
+        return False
+    d = scipy.spatial.distance_matrix(xs, xs)
+    d = np.extract(1 - np.identity(d.shape[0]), d)
+    return not np.any(d < 1e-3 / c)
+
+
+@st.composite
+def point_inside_simplex(draw, simplex):
+    simplex = draw(simplex)
+    simplex = np.asarray(simplex)
+    dim = simplex.shape[1]
+    # Set the numpy random seed
+    draw(st.random_module())
+    # 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):
+    kwargs = dict(
+        allow_nan=False, allow_infinity=False, exclude_min=True, exclude_max=True
+    )
+    if isinstance(domain, Interval):
+        a, b = domain.bounds
+        eps = (b - a) * 1e-2
+        x = st.floats(min_value=(a + eps), max_value=(b - eps), **kwargs)
+    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]
+        )
+        x = point_inside_simplex(simplex)
+
+    xs = st.tuples(*[x] * n).filter(unique_vectors)
+    return draw(xs)
+
+
+@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):
+    kwargs = dict(allow_nan=False, allow_infinity=False)
+    if isinstance(domain, Interval):
+        a, b = domain.bounds
+        length = b - a
+        before_domain = st.floats(a - 10 * length, a, exclude_max=True, **kwargs)
+        after_domain = st.floats(b, b + 10 * length, exclude_min=True, **kwargs)
+        x = before_domain | after_domain
+    else:
+        assert isinstance(domain, ConvexHull)
+        hull = domain.bounds
+        # Generate point between bounding box and bounding box * 10
+        points = hull.points[hull.vertices]
+        x = st.tuples(
+            *[
+                (
+                    st.floats(a - 10 * (b - a), a, exclude_max=True, **kwargs)
+                    | st.floats(b, b + 10 * (b - a), exclude_min=True, **kwargs)
+                )
+                for a, b in zip(points.min(axis=0), points.max(axis=0))
+            ]
+        )
+
+    xs = st.tuples(*[x] * n).filter(unique_vectors)
+    return draw(xs)
+
+
+@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
+
+    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])
+
+    f = st.floats(1e-3, 1 - 1e-3, allow_nan=False, allow_infinity=False)
+    xb = draw(st.tuples(*[f] * 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):
+    if ndim == 1:
+        limits = draw(st.tuples(reals, reals).map(sorted).filter(lambda x: x[0] < x[1]))
+        domain = Interval(*limits)
+    else:
+        points = draw(
+            hynp.arrays(np.float, (10, ndim), elements=reals, unique=True).filter(
+                unique_vectors
+            )
+        )
+        domain = ConvexHull(points)
+    return domain
+
+
+@st.composite
+def make_hypercube_domain(draw, ndim, fill=True):
+    if ndim == 1:
+        limit = draw(positive_reals)
+        subdomain = Interval(-limit, limit)
+    else:
+        x = draw(positive_reals)
+        point = np.full(ndim, x)
+        boundary_points = [r @ point for r in reflections(ndim)]
+        subdomain = ConvexHull(boundary_points)
+    return subdomain