diff --git a/.gitignore b/.gitignore index a4ba2335f..a562a9918 100644 --- a/.gitignore +++ b/.gitignore @@ -44,6 +44,7 @@ nosetests.xml coverage.xml *,cover .hypothesis/ +.pytest_cache/ # Translations *.mo diff --git a/adaptive/learner/learnerND.py b/adaptive/learner/learnerND.py index 33bbbbb07..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): + def __init__(self, func, bounds, loss_per_simplex=None, *, anisotropic=False): self._vdim = None self.loss_per_simplex = loss_per_simplex or default_loss @@ -363,6 +363,7 @@ def __init__(self, func, bounds, loss_per_simplex=None): # create a private random number generator with fixed seed self._random = random.Random(1) + self.anisotropic = anisotropic # all real triangles that have not been subdivided and the pending # triangles heap of tuples (-loss, real simplex, sub_simplex or None) @@ -379,7 +380,12 @@ def __init__(self, func, bounds, loss_per_simplex=None): def new(self) -> LearnerND: """Create a new learner with the same function and bounds.""" - return LearnerND(self.function, self.bounds, self.loss_per_simplex) + return LearnerND( + self.function, + self.bounds, + self.loss_per_simplex, + anisotropic=self.anisotropic, + ) @property def npoints(self): @@ -529,7 +535,6 @@ def points(self): def tell(self, point, value): point = tuple(point) - if point in self.data: return # we already know about the point @@ -548,9 +553,54 @@ def tell(self, point, value): simplex = self._pending_to_simplex.get(point) if simplex is not None and not self._simplex_exists(simplex): simplex = None + # Keep Bowyer-Watson in the globally scaled coordinates. The + # anisotropic per-simplex transform is only used for point + # selection, not for cavity construction. to_delete, to_add = tri.add_point(point, simplex, transform=self._transform) self._update_losses(to_delete, to_add) + def get_local_transform_matrix(self, simplex): + scale = self._transform + + if simplex is None or self.tri is None or not self.anisotropic: + return scale + + neighbors = set.union(*[self.tri.vertex_to_simplices[i] for i in simplex]) + indices = set.union(set(), *neighbors) + + points = np.array([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 + points = points @ scale + 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)[0] + *gradient, _constant = fit + # we do not need the constant, only the gradient + + # 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) - 1 + + scale_along_gradient = projection_matrix * factor + identity + m = np.dot(scale_along_gradient, scale) + return m.T + def _simplex_exists(self, simplex): simplex = tuple(sorted(simplex)) return simplex in self.tri.simplices @@ -679,22 +729,29 @@ def _pop_highest_existing_simplex(self): def _ask_best_point(self): assert self.tri is not None - loss, simplex, subsimplex = self._pop_highest_existing_simplex() + while True: + loss, simplex, subsimplex = self._pop_highest_existing_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) + 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) - point_new = tuple(choose_point_in_simplex(points, transform=self._transform)) + transform = self.get_local_transform_matrix(simplex) + point_new = tuple(choose_point_in_simplex(points, transform=transform)) + + # The priority queue may still contain stale simplices from before a + # pending split. Skip those rather than returning duplicate work. + if point_new in self.data or point_new in self.pending_points: + continue - self._pending_to_simplex[point_new] = simplex - self.tell_pending(point_new, simplex=simplex) # O(??) + self._pending_to_simplex[point_new] = simplex + self.tell_pending(point_new, simplex=simplex) # O(??) - return point_new, loss + return point_new, loss @property def _bounds_available(self): diff --git a/adaptive/learner/triangulation.py b/adaptive/learner/triangulation.py index 26a5ebc2a..8cf377248 100644 --- a/adaptive/learner/triangulation.py +++ b/adaptive/learner/triangulation.py @@ -3,6 +3,7 @@ from itertools import chain, combinations from math import factorial, sqrt +import numpy as np import scipy.spatial from numpy import abs as np_abs from numpy import ( @@ -252,6 +253,8 @@ def simplex_volume_in_embedding(vertices) -> float: 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 @@ -566,7 +569,14 @@ def bowyer_watson(self, pt_index, containing_simplex=None, transform=None): self.add_simplex(simplex) 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 + + 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 _simplex_is_almost_flat(self, simplex): return self._relative_volume(simplex) < 1e-8 @@ -713,3 +723,45 @@ 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, _point 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) diff --git a/adaptive/tests/test_learnernd.py b/adaptive/tests/test_learnernd.py index 0884b7eeb..5cfe3feb2 100644 --- a/adaptive/tests/test_learnernd.py +++ b/adaptive/tests/test_learnernd.py @@ -2,12 +2,18 @@ import scipy.spatial from adaptive.learner import LearnerND -from adaptive.runner import replay_log, simple +from adaptive.runner import BlockingRunner, replay_log, simple from .test_learners import generate_random_parametrization, ring_of_fire -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), @@ -26,6 +32,24 @@ def test_faiure_case_LearnerND(): replay_log(learner, log) +def test_anisotropic_3d(): + # There was a bug where the total simplex volume would exceed the bounding + # box volume for the anisotropic 3D learner. + # 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(learner): + if learner.tri: + sum_of_simplex_volumes = np.sum(learner.tri.volumes()) + assert sum_of_simplex_volumes < 8.00000000001 + return learner.npoints >= 1000 + + BlockingRunner(learner, goal, ntasks=1) + + assert learner.npoints >= 1000 + + def test_interior_vs_bbox_gives_same_result(): f = generate_random_parametrization(ring_of_fire) 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): 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,